GISS Dynamic ocean model
This model is a fully dynamic, non-Boussinesq, mass-conserving
free surface ocean model (Russell et al, 1995, 2000, Liu et al 2002,
2003). The dynamics is based on a modified Arakawa scheme on the
C-grid, with a linear upstream scheme for advecting tracers. Vertical
mixing uses the KPP scheme of Large et al (1996). Momentum mixing is
modelled as a spatically varying laplacian as in Wajsowicz (1993). The
effects of mesoscale eddies and isopycnal diffusion are parameterised
as in Gent and McWilliams (1996), but using variable coefficients
(Visbeck et al, 1997), and coded as in Griffies (1998).
The model contains up to 12 variable depth subgrid scale straits
which contect ocean grid cells, which would not be connected at the
resolution used. In particular, the Straits of Gibraltar, Hormuz, and
Nares straits are so modelled. All ocean components are fluxed through
these straits as a function of the end to end pressure gradients,
balanced against a drag proportional to the straits 'width' which
serves as a tuning parameter to get reasonable fluxes.
Prognostic Variables
MO(I,J,L) (kg/m²) = salt water mass per unit area
G0MO(I,J,L) (J) = potential enthapy in grid cell
GXMO,GYMO,GZMO (J) = linear gradients of potential enthalpy
GXXMO,GYYMO,GZZMO,GXYMO,GYZMO,GZXMO (J) = second order moments of potential enthalpy
S0MO(I,J,L) (kg) = salt mass in grid cell
SXMO,SYMO,SZMO (J) = linear gradients of salt
SXXMO,SYYMO,SZZMO,SXYMO,SYZMO,SZXMO (J) = second order moments of salt
UO (I,J,L) (m/s) = eastward C-grid velocity centered at eastern edge of grid cell
VO (I,J,L) (m/s) = northward C-grid velocity centered at northern edge of grid cell
UOD (I,J,L) (m/s) = eastward D-grid velocity centered at northern edge of grid cell
VOD (I,J,L) (m/s) = northward D-grid velocity centered at eastern edge of grid cell
MUST (L,N) (kg/s) = mass flux of salt water through strait
G0MST (L,N) (J) = potential enthalpy in strait cell
GXMST,GZMST (J) = linear gradients of potential enthapy
S0MST (L,N) (kg) = salt mass in strait cell
SXMST,SZMST (kg) = linear gradients of salt
Cross sectional area of salinity as a function of longitude in an ocean grid cell is:
S (kg/kg) = S0MO / M + 2 * SXMO * m / M² + 6 * SXXMO * m² / M³ - .5 * SXXMO / M
where M (kg) = MO * oDXYP and - .5 M < m < .5 M
Ocean Functions and the Equation of State
The following variables are defined for the ocean model:
C (J/kg°C) = specific heat capacity at constant pressure
g (m/s²) = Earth's vertical gravitational acceleration
G (J/kg) = potential specific enthalpy
H (J/kg) = specific enthalpy
m (kg/m²) = vertical coordinate = mass above minus mean atmospheric mass
P (Pa) = pressure above mean atmospheric pressure
Po (Pa) = global mean atmospheric pressure at sea level (101325 Pa)
S (kg/kg) = salinity
T (°C) = in-situ temperature
α (m³/kg) = specific volume
β (m³/kg) = potential specific volume
Γ (°C/Pa) = adiabatic lapse rate
Θ (°C) = potential temperature
G and S can be calculated at any location in the ocean from m and
the prognostic variables. The purpose of this section is to derive P,
H, T, α and β from G, S and m or G, S and P. Remember that P = m·g.
α(T,S,P), C(T,S,P) and Γ(T,S,P) are provided by Fofonoff and
Millard [1983], but also see Fofonoff [1985]. H(T,S,P) is derived as
follows:
∂H/∂P = α(T,S,P) (Eq. 1)
∂T/∂P = Γ(T,S,P) (Eq. 2)
are integrated simultaneously from P to Po at constant entropy and
salinity to obtain H(T,S,P) - H(Θ,S,Po). Po is the reference pressure
for potential temperature so that Θ = T at Po. Next evaluate H(Θ,S,Po)
- H(Θ,0,Po) at constant entropy and pressure which is provided by
Millero and Leung [1976]. Finally, the differential equation
∂H/∂T = C(T,0,Po) (Eq. 3)
is integrated from Θ to 0°C at constant salinity and pressure to obtain
H(Θ,0,Po) - H(0,0,Po). The reference value for specific enthalpy is
set by H(0,0,Po) = 0. Thus H(T,S,P) is determined by the above formulas.
Θ(T,S,P) and T(Θ,S,P) are derived by integrating (Eq. 2) from P to
Po. The G of a parcel is defined to be the H of the parcel at Po were
it raised adiabatically. Thus, G(T,S,P) = H(Θ(T,S,P),S,Po) = G(Θ,S).
Θ(G,S) is derived by inverting the first argument of G(Θ,S) which is
possible because G and Θ are in a one-to-one relationship at constant S and Po.
T, H, α and β can now be derived as functions of G, S and P.
T(G,S,P) = T(G(Θ,S),S,P) = T(Θ,S,P)
H(G,S,P) = H(G(Θ,S),S,P) = H(Θ,S,P) = H(Θ(T,S,P),S,P) = H(T,S,P)
α(G,S,P) = α(G(Θ,S),S,P) = α(Θ,S,P) = α(Θ(T,S,P),S,P) =
α(T,S,P)
β(G,S) = α(G,S,Po)
Evaluation of α(G,S,P), called the equation of state, is
performed by linear interpolation within a look-up table in the ocean
model's computer program. The increments in the table are 4000 (J/kg)
for G, .001 for S, and 2·10^5 (Pa) for P.
Ocean Pressure Gradient Force
The pressure gradient force accelerates the velocity between two
mass grid cells at the same layer. The force in a layer has two terms:
the gradient of the mean pressure multiplied by the average mean volume
of the two cells and the gradient of the mean geopotential multiplied
by the average mass of the two cells. Other gradients of
potential enthalpy and salt are not used in calculating the above mean
quantities, but the linear vertical gradients are used.
The mean vertical coordinate, m1, of a grid cell is the average of
the top edge m and the bottom edge m. Δm1 is the difference of m at
the two edges. The mean pressure of a grid cell is P(m1) = m1·g, where g
(m/s²) is the uniform effective gravity. To calculate the mean
specific volume of a grid cell we assume that α(G(m),S(m),P(m)) = α(m)
fits a quadratic polynomial in m. The mean value of the quadratic α(m)
from m1-Δm1/2 to m1+Δm1/2 is equal to [α(m1-x) + α(m1+x)]/2 where x =
Δm1/√12. The mean volume of a grid cell is equal to the
mass multiplied by the mean specific volume which is calculated using
quadratic precision using two evaluations of α(G,S,P) per grid cell at
the vertical coordinates m1-x and m1+x. This shows how the first term
of the pressure gradient term is modeled.
The second term is modeled as follows. We again assume that α(m)
fits a quadratic in each grid cell. The height, h (meters), at the
column bottom is specified by the bottom topography. Integrating
upwards, the change in h from the bottom layer edge to any level in a
layer is equal to the integral of α(m) dm. The change in h to the top
layer edge uses the two evaluations of α mentioned in the previous
paragraph and is calculated with quadratic precision. The top edge
value of h of a layer becomes the bottom edge value of h for the layer
above. The mean geopotential, Φ (m²/s²), of the layer is the mass
weighted g·h throughout the layer which is a double integral of α(m) and
is equal to
∫m1+Δm1/2
Φ = 1/Δm1 ∫ g·h(m) dm =
∫m1-Δm1/2
∫m1+Δm1/2 [ ∫m1+Δm1/2 ]
= 1/Δm1 ∫ g [ h(m1+Δm1/2) + ∫ α(n) dn ] dm =
∫m1-Δm1/2 [ ∫m ]
= g·h(m1+Δm1/2) + g[α(m1-x)(3-√3) + α(m1+x)(3+√3)]m1/12
The mean Φ is calculated with quadratic precision using the same
two evaluations of α(G,S,P) as used in evaluating the first term.
Given an arbitrary quadratic polynomial in an interval
[-Δm1/2,Δm1/2], the abcissas of the intersection points between the
polynomial and the least square fit line that fits the polynomial in
the interval are the points -x and x mentioned above. These points do
not depend on the coefficients of the polynomial. Conversely, there
are no two other points such that the mean value of an arbitrary
quadratic polynomial can be derived from evaluating the polynomial at
the two points. This mathematical result is even more elegant because
the double integral above can be calculated from an evaluation of
the polynomial at the points -x and x.
2024/03/28/03:30:20