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