Implementation of the GM/Redi/Griffies schemes in the MIT OGCM

A.Adcroft and C.Hill, MIT

July 2000

1 Overview

There are two parts to the Redi/GM parameterization of geostrophic eddies. The first aims to mix tracer properties along isentropes (neutral surfaces) by means of a diffusion operator oriented along the local isentropic surface (Redi). The second part, adiabatically re-arranges tracers through an advective flux where the advecting flow is a function of slope of the isentropic surfaces (GM).

The first GCM implementation of the Redi scheme was by Cox 1987 in the GFDL ocean circulation model. The original approach failed to distinguish between isopycnals and surfaces of locally referenced potential density (now called neutral surfaces) which are proper isentropes for the ocean. As will be discussed later, it also appears that the Cox implementation is susceptible to a computational mode. Due to this mode, the Cox scheme requires a background lateral diffusion to be present to conserve the integrity of the model fields.

The GM parameterization was then added to the GFDL code in the form of a non-divergent bolus velocity. The method defines two streamfunctions expressed in terms of the isoneutral slopes subject to the boundary condition of zero value on upper and lower boundaries. The horizontal bolus velocities are then the vertical derivative of these functions. Here in lies a problem highlighted by Griffies et al., 1997: the bolus velocities involve multiple derivatives on the potential density field, which can consequently give rise to noise. Griffies et al. point out that the GM bolus fluxes can be identically written as a skew flux which involves fewer differential operators. Further, combining the skew flux formulation and Redi scheme, substantial cancellations take place to the point that the horizontal fluxes are unmodified from the lateral diffusion parameterization.

2 Redi scheme: Isopycnal diffusion

The Redi scheme diffuses tracers along isopycnals and introduces a term in the tendency (rhs) of such a tracer (here τ) of the form:

∇ · κ_{ρ}K_{Redi}∇τ  (1)

where κ_{ρ} is the along isopycnal diffusivity and K_{Redi} is a rank 2 tensor that projects the gradient of τ onto the isopycnal surface. The unapproximated projection tensor is:

K_{Redi} = ( \begin{matrix}1 + S_{x} & S_{x}S_{y} & S_{x} \\ S_{x}S_{y} & 1 + S_{y} & S_{y} \\ S_{x} & S_{y} & |S|^{2} \end{matrix} )  (2)

Here, S_{x} = −∂_{x}σ/∂_{z}σ and S_{y} = −∂_{y}σ/∂_{z}σ are the components of the isoneutral slope.

The first point to note is that a typical slope in the ocean interior is small, say of the order 10^{−4}. A maximum slope might be of order 10^{−2} and only exceeds such in unstratified regions where the slope is ill defined. It is therefore justifiable, and customary, to make the small slope approximation, |S| << 1. The Redi projection tensor then becomes:

K_{Redi} = ( \begin{matrix}1 & 0 & S_{x} \\ 0 & 1 & S_{y} \\ S_{x} & S_{y} & |S|^{2} \end{matrix} )  (3)

3 GM parameterization

The GM parameterization aims to parameterise the “advective” or “transport” effect of geostrophic eddies by means of a “bolus” velocity, u^{∗}. The divergence of this advective flux is added to the tracer tendency equation (on the rhs):

−∇ · τu^{∗}  (4)

The bolus velocity is defined as:

u^{∗} = −∂_{z}F_{x}  (5)
v^{∗} = −∂_{z}F_{y}  (6)
w^{∗} = ∂_{x}F_{x} + ∂_{y}F_{y}  (7)

where F_{x} and F_{y} are stream-functions with boundary conditions F_{x} = F_{y} = 0 on upper and lower boundaries. The virtue of casting the bolus velocity in terms of these stream-functions is that they are automatically non-divergent (∂_{x}u^{∗} + ∂_{y}v^{∗} = −∂_{xz}F_{x} − ∂_{yz}F_{y} = −∂_{z}w^{∗}). F_{x} and F_{y} are specified in terms of the isoneutral slopes S_{x} and S_{y}:

F_{x} = κ_{GM}S_{x}  (8)
F_{y} = κ_{GM}S_{y}  (9)

This is the form of the GM parameterization as applied by Donabasaglu, 1997, in MOM versions 1 and 2.

3.1 Griffies Skew Flux

Griffies notes that the discretisation of bolus velocities involves multiple layers of differencing and interpolation that potentially lead to noisy fields and computational modes. He pointed out that the bolus flux can be re-written in terms of a non-divergent flux and a skew-flux:

u^{∗}τ = ( \begin{matrix}−∂_{z}(κ_{GM}S_{x})τ \\ −∂_{z}(κ_{GM}S_{y})τ \\ (∂_{x}κ_{GM}S_{x} + ∂_{y}κ_{GM}S_{y})τ \end{matrix} )

= ( \begin{matrix}−∂_{z}(κ_{GM}S_{x}τ) \\ −∂_{z}(κ_{GM}S_{y}τ) \\ ∂_{x}(κ_{GM}S_{x}τ) + ∂_{y}(κ_{GM}S_{y})τ) \end{matrix} ) + ( \begin{matrix}κ_{GM}S_{x}∂_{z}τ \\ κ_{GM}S_{y}∂_{z}τ \\ −κ_{GM}S_{x}∂_{x}τ − κ_{GM}S_{y})∂_{y}τ \end{matrix} )

The first vector is non-divergent and thus has no effect on the tracer field and can be dropped. The remaining flux can be written:

u^{∗}τ = −κ_{GM}K_{GM}∇τ  (10)

where

K_{GM} = ( \begin{matrix}0 & 0 & −S_{x} \\ 0 & 0 & −S_{y} \\ S_{x} & S_{y} & 0 \end{matrix} )  (11)

is an anti-symmetric tensor.

This formulation of the GM parameterization involves fewer derivatives than the original and also involves only terms that already appear in the Redi mixing scheme. Indeed, a somewhat fortunate cancellation becomes apparent when we use the GM parameterization in conjunction with the Redi isoneutral mixing scheme:

κ_{ρ}K_{Redi}∇τ − u^{∗}τ = (κ_{ρ}K_{Redi} + κ_{GM}K_{GM})∇τ  (12)

In the instance that κ_{GM} = κ_{ρ} then

κ_{ρ}K_{Redi} + κ_{GM}K_{GM} = κ_{ρ} ( \begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 2S_{x} & 2S_{y} & |S|^{2} \end{matrix} )  (13)

which differs from the variable laplacian diffusion tensor by only two non-zero elements in the z-row.

4 Variable κ_{GM}

Visbeck et al., 1996, suggest making the eddy coefficient, κ_{GM}, a function of the Eady growth rate, |f|/\sqrt{Ri}. The formula involves a non-dimensional constant, α, and a length-scale L:

κ_{GM} = αL^{2} \overline{\frac{|f|}{\sqrt{Ri}}}^{z}

where the Eady growth rate has been depth averaged (indicated by the overline). A local Richardson number is defined Ri = N^{2}/(∂u/∂z)^{2} which, when combined with thermal wind gives:

\frac{1}{Ri} = \frac{(\frac{∂u}{∂z})^{2}}{N^{2}} = \frac{(\frac{g}{fρ_{o}} |∇σ|)^{2}}{N^{2}} = \frac{M^{4}}{|f|^{2}N^{2}}

where M^{2} is defined M^{2} = \frac{g}{ρ_{o}} |∇σ|. Substituting into the formula for κ_{GM} gives:

κ_{GM} = αL^{2} \overline{\frac{M^{2}}{N}}^{z} = αL^{2} \overline{\frac{M^{2}}{N^{2}} N}^{z} = αL^{2}\overline{|S|N}^{z}

5 Tapering and stability

Experience with the GFDL model showed that the GM scheme has to be matched to the convective parameterization. This was originally expressed in connection with the introduction of the KPP boundary layer scheme (Large et al., 97) but infact, as subsequent experience with the MIT model has found, is necessary for any convective parameterization.

Deep convection sites and the mixed layer are indicated by homogenized, unstable or nearly unstable stratification. The slopes become in such regions are infinite, very large with a sign reversal or simply very large. From a numerical point of view, large slopes lead to large variations in the tensor elements (implying large bolus flow) and can be numerically unstable. This was first reognized by Cox, 1987, who implemented “slopes clipping” in the isopycnal mixing tensor. Here, the slope magnitude is simply restricted by an upper limit:

|∇σ| = \sqrt{σ_{x}^{2} + σ_{y}^{2}}  (14)
S_{lim} = −\frac{|∇σ|}{S_{max}}  where S_{max} is a parameter  (15)
σ_{z}^{⋆} = min(σ_{z}, S_{lim})  (16)
[s_{x}, s_{y}] = −\frac{[σ_{x}, σ_{y}]}{σ_{z}^{⋆}}  (17)

Notice that this algorithm assumes stable stratification through the “min” function. In the case where the fluid is well stratified (σ_{z} < S_{lim}) then the slopes evaluate to:

[s_{x}, s_{y}] = −\frac{[σ_{x}, σ_{y}]}{σ_{z}}  (18)

while in the limited regions (σ_{z} > S_{lim}) the slopes become:

[s_{x}, s_{y}] = \frac{[σ_{x}, σ_{y}]}{|∇σ|/S_{max}}  (19)

so that the slope magnitude is limited \sqrt{s_{x}^{2} + s_{y}^{2}} = S_{max}.

Even using slope clipping, it is normally the case that the vertical diffusion term (with coefficient κ_{ρ}K_{33} = κ_{ρ}S_{max}^{2}) is large and must be time-stepped using an implicit procedure (see section on discretisation and code later). Fig. ?? shows the mixed layer depth resulting from a) using the GM scheme with clipping and b) no GM scheme (horizontal diffusion). The classic result of dramatically reducd mixed layers is evident. Indeed, the deep convection sites to just one or two points each and are much shallower than we might prefer.

Figure 1: Mixed layer depth using GM parameterization with a) Cox slope clipping and for comparison b) using horizontal constant diffusion. ??

Figure 2: Effective slope as a function of “true” slope using a) Cox slope clipping, b) GKW91 limiting, c) DM95 limiting and d) LDD97 limiting.

This it turns out is due to the over zealous restratification due to the bolus transport parameterization.

Gerdes et al., 1991 ....

Danabasoglu and McWilliams, 1995 ....

Large et al., 1997 ....

Issue: should GM and Redi be tapered? cf convective paper

Issue: is adiabtic important in these regions? changes discretization

6 Discretisation and code

The Gent-McWilliams-Redi parameterization is implemented through the package “gmredi”. There are two necessary calls to “gmredi” routines other than initialization; 1) to calculate the slope tensor as a function of the current model state (gmredi calc tensor) and 2) evaluation of the lateral and vertical fluxes due to gradients along isopycnals or bolus transport (gmredi xtransport, gmredi ytransport and gmredi-rtransport).

Each element of the tensor is discretised to be adiabatic and so that there would be no flux if the gmredi operator is applied to buoyancy. To acheive this we have to consider both these constraints for each row of the tensor,

Figure 3: Effective slope magnitude at 100 m depth evaluated using a) Cox slope clipping, b) GKW91 limiting, c) DM95 limiting and d) LDD97 limiting. each row corresponding to a ’u’, ’v’ or ’w’ point on the model grid.

This is the old documentation.....

The code that implements the Redi/GM/Griffies schemes involves an original core routine inc tracer() that is used to calculate the tendency in the tracers (namely, salt and potential temperature) and a new routine RediTensor() that calculates the tensor components and κ_{GM}.

6.1 subroutine RediTensor()

subroutine RediTensor(Temp,Salt,Kredigm,K31,K32,K33, nIter,DumpFlag)
                      |---in--| |-------out-------|
! Input
real Temp(Nx,Ny,Nz)     ! Potential temperature
real Salt(Nx,Ny,Nz)     ! Salinity
! Output
real Kredigm(Nx,Ny,Nz)  ! Redi/GM eddy coefficient
real K31(Nx,Ny,Nz)      ! Redi/GM (3,1) tensor component
real K32(Nx,Ny,Nz)      ! Redi/GM (3,2) tensor component
real K33(Nx,Ny,Nz)      ! Redi/GM (3,3) tensor component
! Auxiliary input
integer nIter           ! interation/time-step number
logical DumpFlag        ! flag to indicate routine should ``dump''

The subroutine RediTensor() is called from model() with input arguments T and S. It returns the 3D-arrays Kredigm, K31, K32 and K33 which represent κ_{GM} (at T/S points) and the three components of the bottom row in the Redi/GM tensor; 2S_{x}, 2S_{y} and |S|^{2} respectively, all at W points.

The discretisations and algorithm within RediTensor() are as follows. The routine first calculates the locally reference potential density σ_{θ} from T and S and calculates the potential density gradients in subroutine grad-Sigma():

Array Grid-point Definition
SigX  U  σ_{x} = \frac{1}{Δx}δ_{x}σ|_{z(k)}
SigY  V  σ_{y} = \frac{1}{Δy}δ_{y}σ|_{z(k)}
SigZ  W  σ_{z} = \frac{1}{Δz}[σ|_{z(k)}(k − 1/2) − σ|_{z(k)}(k + 1/2)]

Note that σ_{z} is the static stability because the potential densities are referenced to the same reference level (W -level).

The next step calculates the three tensor components K13, K23 and K33 in subroutine KtensorWface(). First, the lateral gradients σ_{x} and σ_{y} are interpolated to the W points and stored in intermediate variables:

Sx = \overline{\overline{σ_{x}}^{x}}^{z}
Sy = \overline{\overline{σ_{y}}^{y}}^{z}

Next, the magnitude of ∇_{z}σ is stored in an intermediate variable:

                       Sxy2 = \sqrt{Sx^{2} + Sy^{2}}

The stratification (σ_{z}) is “checked” such that the slope vector has magnitude less than or equal to Smax and stored in an intermediate variable:

                    Sz = max(σ_{z}, −Sxy2/Smax)

This guarantees stability and at the same time retains the lateral orientation of the slope vector. The tensor components are then calculated:

                   K13 = −2Sx/Sz
                   K23 = −2Sx/Sz
                   K33 = (Sx/Sz)^{2} + (Sy/Sz)^{2}

Finally, Kredigm (κ_{GM}) is calculated in subroutine GMRediCoefficient(). First, all the gradients are interpolated to the T/S points and stored in intermediate variables:

Sx = \overline{σ_{x}}^{x}
Sy = \overline{σ_{y}}^{y}
Sz = \overline{σ_{z}}^{z}

Again, a nominal stratification is found by “check” the magnitude of the slope vector but here is converted to a Brunt-Vasala frequency:

                  M2 = \sqrt{Sx^{2} + Sy^{2}}

N2 = −\frac{g}{ρ_{o}} max(Sz, −M2/Smax

The magnitude of the slope is then |S| = M2/N2. The Eady growth rate is defined as |f|/\sqrt{(}Ri) = |S|N and is calculated as:

                       FrRi = \frac{M2}{N2}(−\frac{g}{ρ}Sz)

The Eady growth rate is then averaged over the upper layers (about 1100m) and κ_{GM} specified from this 2D-variable:

Kredigm = 0.02 ∗ (200d3 ∗ ∗2) ∗ FrRi

6.2 subroutine inc tracer()

inctracer() is called from model() and has four new arguments:

subroutine inc_tracer( ...,Kredigm,K31,K32,K33, ... )
real Kredigm(Nx,Ny,Nz)  ! Eddy coefficient
real K31(Nx,Ny,Nz)      ! (3,1) tensor coefficient
real K32(Nx,Ny,Nz)      ! (3,2) tensor coefficient
real K33(Nx,Ny,Nz)      ! (3,3) tensor coefficient

Within the routine, the lateral fluxes, fluxWest and fluxSouth, in the Redi/GM/Griffies scheme are very similar to the conventional horizontal diffusion terms except that the diffusion coefficient is a function of space and must be interpolated from the T/S points:

                fluxWest(τ) = . . . + \overline{Kredigm}^{x}∂_{x}τ
               fluxSouth(τ) = . . . + \overline{Kredigm}^{y}∂_{y}τ

The Redi/GM/Griffies scheme adds three terms to the vertical flux (fluxUpper) in the tracer equation. It is discretise simply:

fluxUpper(τ) = … + \overline{Kredigm}^{z} (K13\overline{∂_{x}τ}^{xz} + K23\overline{∂_{y}τ}^{yz} + K33∂_{z}τ )

On boundaries, fluxUpper is set to zero.

7 Outstanding issues?

A second innovation discussed by Griffies concerns the nonlinearity in the equation of state and its implication for the evaluation of the isoneutral slopes. He argues that ∇σ should be written ∇σ = ρ_{o}(β∇S − α∇θ) so that, for example, the zonal slope becomes:

S_{x} = \frac{α∂_{x}θ − β∂_{x}S}{β∂_{z}S − α∂_{z}θ}

where α = \frac{1}{ρ_{o}} \frac{∂ρ}{∂θ} and β = \frac{1}{ρ_{o}} \frac{∂ρ}{∂S} are locally evaluated. He also advocates a very complicated interpolation scheme which he calls the “triads” scheme. Both of these innovations are designed to overcome computational modes and at the same time make the scheme more adiabatic.

Compared to the current implementation out-lined above, the Griffies innovations requires two non-linear equations to replace the equation of state, twice as many differences because of the explicit expressions with θ and S, and a much wider stencil (ie. much more computation) for the “triads” interpolation.

We have not implemented these aspects of Griffies schemes on the basis that the additional computational costs and code complexity, so far, outweigh the immediate benefits.