s_phys_pkgs/text/gmredi.tex

\subsection{GMREDI: Gent-McWilliams/Redi SGS Eddy Parameterization}
\label{sec:pkg:gmredi}
\begin{rawhtml}
<!-- CMIREDIR:gmredi: -->
\end{rawhtml}

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

The first GCM implementation of the Redi scheme was by \cite{Cox87} 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
stream-functions 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
\cite{gretal:98}: 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.

\subsubsection{Redi scheme: Isopycnal diffusion}

The Redi scheme diffuses tracers along isopycnals and introduces a
term in the tendency (rhs) of such a tracer (here $\tau$) of the form:
\begin{equation}
\bf{\nabla} \cdot \kappa_\rho \bf{K}_{Redi}  \bf{\nabla} \tau
\end{equation}
where $\kappa_\rho$ is the along isopycnal diffusivity and
$\bf{K}_{Redi}$ is a rank 2 tensor that projects the gradient of
$\tau$ onto the isopycnal surface. The unapproximated projection tensor is:
\begin{equation}
\bf{K}_{Redi} = \frac{1}{1 + |S|^2} \left(
\begin{array}{ccc}
1 + S_y^2& -S_x S_y & S_x \\
-S_x S_y  & 1 + S_x^2 & S_y \\
S_x & S_y & |S|^2 \\
\end{array}
\right)
\end{equation}
Here, $S_x = -\partial_x \sigma / \partial_z \sigma$ and $S_y =
-\partial_y \sigma / \partial_z \sigma$ 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:
\begin{equation}
\bf{K}_{Redi} = \left(
\begin{array}{ccc}
1 & 0 & S_x \\
0 & 1 & S_y \\
S_x & S_y & |S|^2 \\
\end{array}
\right)
\end{equation}


\subsubsection{GM parameterization}

The GM parameterization aims to represent the ``advective'' or
``transport'' effect of geostrophic eddies by means of a ``bolus''
velocity, $\bf{u}^\star$. The divergence of this advective flux is added
to the tracer tendency equation (on the rhs):
\begin{equation}
- \bf{\nabla} \cdot \tau \bf{u}^\star
\end{equation}

The bolus velocity $\bf{u}^\star$ is defined as the rotational of a
streamfunction $\bf{F}^\star$=$(F_x^\star,F_y^\star,0)$:
\begin{equation}
\bf{u}^\star = \nabla \times \bf{F}^\star =
\left( \begin{array}{c}
- \partial_z  F_y^\star \\
+ \partial_z  F_x^\star \\
\partial_x F_y^\star - \partial_y F_x^\star
\end{array} \right),
\end{equation}
and thus is automatically non-divergent. In the GM parameterization, the streamfunction is
specified in terms of the isoneutral slopes $S_x$ and $S_y$:
\begin{eqnarray}
F_x^\star & = & -\kappa_{GM} S_y \\
F_y^\star & = &  \kappa_{GM} S_x
\end{eqnarray}
with boundary conditions $F_x^\star=F_y^\star=0$ on upper and lower boundaries.
In the end, the bolus transport in the GM parameterization is given by:
\begin{equation}
\bf{u}^\star = \left(
\begin{array}{c}
u^\star \\
v^\star \\
w^\star
\end{array}
\right) = \left(
\begin{array}{c}
- \partial_z (\kappa_{GM} S_x) \\
- \partial_z (\kappa_{GM} S_y) \\
\partial_x  (\kappa_{GM} S_x) + \partial_y (\kappa_{GM} S_y)
\end{array}
\right)
\end{equation}

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

Note that in the MITgcm, the variables containing the GM bolus streamfunction are:
\begin{equation}
\left(
\begin{array}{c}
GM\_PsiX \\
GM\_PsiY
\end{array}
\right) = \left(
\begin{array}{c}
\kappa_{GM} S_x \\
\kappa_{GM} S_y
\end{array}
\right)= \left(
\begin{array}{c}
F_y^\star \\
-F_x^\star
\end{array}
\right).
\end{equation}
  
\subsubsection{Griffies Skew Flux}

\cite{gr:98} 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:
\begin{eqnarray*}
\bf{u}^\star \tau
& = &
\left( \begin{array}{c}
- \partial_z ( \kappa_{GM} S_x ) \tau \\
- \partial_z ( \kappa_{GM} S_y ) \tau \\
(\partial_x \kappa_{GM} S_x + \partial_y \kappa_{GM} S_y)\tau
\end{array} \right)
\\
& = &
\left( \begin{array}{c}
- \partial_z ( \kappa_{GM} S_x \tau) \\
- \partial_z ( \kappa_{GM} S_y \tau) \\
\partial_x ( \kappa_{GM} S_x \tau) + \partial_y ( \kappa_{GM} S_y \tau)
\end{array} \right)
+ \left( \begin{array}{c}
 \kappa_{GM} S_x \partial_z \tau \\
 \kappa_{GM} S_y \partial_z \tau \\
- \kappa_{GM} S_x \partial_x \tau - \kappa_{GM} S_y \partial_y \tau
\end{array} \right)
\end{eqnarray*}
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:
\begin{equation}
\bf{u}^\star \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau
\end{equation}
where
\begin{equation}
\bf{K}_{GM} =
\left(
\begin{array}{ccc}
0 & 0 & -S_x \\
0 & 0 & -S_y \\
S_x & S_y & 0
\end{array}
\right)
\end{equation}
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:
\begin{equation}
\kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
- u^\star \tau = 
( \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} ) \bf{\nabla} \tau
\end{equation}
In the instance that $\kappa_{GM} = \kappa_{\rho}$ then
\begin{equation}
\kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} =
\kappa_\rho
\left( \begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
2 S_x & 2 S_y & |S|^2 
\end{array}
\right)
\end{equation}
which differs from the variable Laplacian diffusion tensor by only
two non-zero elements in the $z$-row.

\fbox{ \begin{minipage}{4.75in}
{\em S/R GMREDI\_CALC\_TENSOR} ({\em pkg/gmredi/gmredi\_calc\_tensor.F})

$\sigma_x$: {\bf SlopeX} (argument on entry)

$\sigma_y$: {\bf SlopeY} (argument on entry)

$\sigma_z$: {\bf SlopeY} (argument)

$S_x$: {\bf SlopeX} (argument on exit)

$S_y$: {\bf SlopeY} (argument on exit)

\end{minipage} }


\subsubsection{Variable $\kappa_{GM}$}

\cite{visbeck:97} suggest making the eddy coefficient,
$\kappa_{GM}$, a function of the Eady growth rate,
$|f|/\sqrt{Ri}$. The formula involves a non-dimensional constant,
$\alpha$, and a length-scale $L$:
\begin{displaymath}
\kappa_{GM} = \alpha L^2 \overline{ \frac{|f|}{\sqrt{Ri}} }^z
\end{displaymath}
where the Eady growth rate has been depth averaged (indicated by the
over-line). A local Richardson number is defined $Ri = N^2 / (\partial
u/\partial z)^2$ which, when combined with thermal wind gives:
\begin{displaymath}
\frac{1}{Ri} = \frac{(\frac{\partial u}{\partial z})^2}{N^2} =
\frac{ ( \frac{g}{f \rho_o} | {\bf \nabla} \sigma | )^2 }{N^2} =
\frac{ M^4 }{ |f|^2 N^2 }
\end{displaymath}
where $M^2$ is defined $M^2 = \frac{g}{\rho_o} |{\bf \nabla} \sigma|$.
Substituting into the formula for $\kappa_{GM}$ gives:
\begin{displaymath}
\kappa_{GM} = \alpha L^2 \overline{ \frac{M^2}{N} }^z =
\alpha L^2 \overline{ \frac{M^2}{N^2} N }^z =
\alpha L^2 \overline{ |S| N }^z
\end{displaymath}


\subsubsection{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 \citep{lar-eta:94} but in fact, as subsequent experience
with the MIT model has found, is necessary for any convective
parameterization.

\fbox{ \begin{minipage}{4.75in}
{\em S/R GMREDI\_SLOPE\_LIMIT} ({\em
pkg/gmredi/gmredi\_slope\_limit.F})

$\sigma_x, s_x$: {\bf SlopeX} (argument)

$\sigma_y, s_y$: {\bf SlopeY} (argument)

$\sigma_z$: {\bf dSigmadRReal} (argument)

$z_\sigma^{*}$: {\bf dRdSigmaLtd} (argument)

\end{minipage} }

\begin{figure}
\begin{center}
\resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/tapers.eps}}
\end{center}
\caption{Taper functions used in GKW91 and DM95.}
\label{fig:tapers}
\end{figure}

\begin{figure}
\begin{center}
\resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/effective_slopes.eps}}
\end{center}
\caption{Effective slope as a function of ``true'' slope using Cox
slope clipping, GKW91 limiting and DM95 limiting.}
\label{fig:effective_slopes}
\end{figure}


\subsubsection*{Slope clipping}

Deep convection sites and the mixed layer are indicated by
homogenized, unstable or nearly unstable stratification. The slopes in
such regions can be either 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 recognized by
\cite{Cox87} who implemented ``slope clipping'' in the isopycnal mixing
tensor. Here, the slope magnitude is simply restricted by an upper
limit:
\begin{eqnarray}
|\nabla \sigma| & = & \sqrt{ \sigma_x^2 + \sigma_y^2 } \\
S_{lim} & = & - \frac{|\nabla \sigma|}{ S_{max} }
\;\;\;\;\;\;\;\; \mbox{where $S_{max}$ is a parameter} \\
\sigma_z^\star & = & \min( \sigma_z , S_{lim} ) \\
{[s_x,s_y]} & = & - \frac{ [\sigma_x,\sigma_y] }{\sigma_z^\star}
\end{eqnarray}
Notice that this algorithm assumes stable stratification through the
``min'' function. In the case where the fluid is well stratified ($\sigma_z < S_{lim}$) then the slopes evaluate to:
\begin{equation}
{[s_x,s_y]} = - \frac{ [\sigma_x,\sigma_y] }{\sigma_z}
\end{equation}
while in the limited regions ($\sigma_z > S_{lim}$) the slopes become:
\begin{equation}
{[s_x,s_y]} = \frac{ [\sigma_x,\sigma_y] }{|\nabla \sigma|/S_{max}}
\end{equation}
so that the slope magnitude is limited $\sqrt{s_x^2 + s_y^2} =
S_{max}$.

The slope clipping scheme is activated in the model by setting {\bf
GM\_tap\-er\_scheme = 'clipping'} in {\em data.gmredi}.

Even using slope clipping, it is normally the case that the vertical
diffusion term (with coefficient $\kappa_\rho{\bf K}_{33} =
\kappa_\rho S_{max}^2$) is large and must be time-stepped using an
implicit procedure (see section on discretisation and code later).
Fig. \ref{fig-mixedlayer} 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 reduced mixed layers is
evident. Indeed, the deep convection sites to just one or two points
each and are much shallower than we might prefer. This, it turns out,
is due to the over zealous re-stratification due to the bolus transport
parameterization. Limiting the slopes also breaks the adiabatic nature
of the GM/Redi parameterization, re-introducing diabatic fluxes in
regions where the limiting is in effect.

\subsubsection*{Tapering: Gerdes, Koberle and Willebrand, Clim. Dyn. 1991}

The tapering scheme used in \cite{gkw:91}
addressed two issues with the clipping method: the introduction of
large vertical fluxes in addition to convective adjustment fluxes is
avoided by tapering the GM/Redi slopes back to zero in
low-stratification regions; the adjustment of slopes is replaced by a
tapering of the entire GM/Redi tensor. This means the direction of
fluxes is unaffected as the amplitude is scaled.

The scheme inserts a tapering function, $f_1(S)$, in front of the
GM/Redi tensor:
\begin{equation}
f_1(S) = \min \left[ 1, \left( \frac{S_{max}}{|S|}\right)^2 \right]
\end{equation}
where $S_{max}$ is the maximum slope you want allowed. Where the
slopes, $|S|<S_{max}$ then $f_1(S) = 1$ and the tensor is un-tapered
but where $|S| \ge S_{max}$ then $f_1(S)$ scales down the tensor so
that the effective vertical diffusivity term $\kappa f_1(S) |S|^2 =
\kappa S_{max}^2$.

The GKW91 tapering scheme is activated in the model by setting {\bf
GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.

\subsubsection*{Tapering: Danabasoglu and McWilliams, J. Clim. 1995}

The tapering scheme used by \cite{dm:95} followed a similar procedure but used a different
tapering function, $f_1(S)$:
\begin{equation}
f_1(S) = \frac{1}{2} \left( 1+\tanh \left[ \frac{S_c - |S|}{S_d} \right] \right)
\end{equation}
where $S_c = 0.004$ is a cut-off slope and $S_d=0.001$ is a scale over
which the slopes are smoothly tapered. Functionally, the operates in
the same way as the GKW91 scheme but has a substantially lower
cut-off, turning off the GM/Redi SGS parameterization for weaker
slopes.

The DM95 tapering scheme is activated in the model by setting {\bf
GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.

\subsubsection*{Tapering: Large, Danabasoglu and Doney, JPO 1997}

The tapering used in \cite{lar-eta:97} is based on the
DM95 tapering scheme, but also tapers the scheme with an additional
function of height, $f_2(z)$, so that the GM/Redi SGS fluxes are
reduced near the surface:
\begin{equation}
f_2(z) = \frac{1}{2} \left( 1 + \sin(\pi \frac{z}{D} - \frac{\pi}{2})\right)
\end{equation}
where $D = L_\rho |S|$ is a depth-scale and $L_\rho=c/f$ with
$c=2$~m~s$^{-1}$.  This tapering with height was introduced to fix
some spurious interaction with the mixed-layer KPP parameterization.

The LDD97 tapering scheme is activated in the model by setting {\bf
GM\_tap\-er\_scheme = 'ldd97'} in {\em data.gmredi}.


\begin{figure}
\begin{center}
%\includegraphics{mixedlayer-cox.eps}
%\includegraphics{mixedlayer-diff.eps}
Figure missing.
\end{center}
\caption{Mixed layer depth using GM parameterization with a) Cox slope
clipping and for comparison b) using horizontal constant diffusion.}
\label{fig-mixedlayer}
\end{figure}

\subsubsection{Package Reference}
\label{sec:pkg:gmredi:diagnostics}

{\footnotesize
\begin{verbatim}
------------------------------------------------------------------------
<-Name->|Levs|<-parsing code->|<--  Units   -->|<- Tile (max=80c) 
------------------------------------------------------------------------
GM_VisbK|  1 |SM P    M1      |m^2/s           |Mixing coefficient from Visbeck etal parameterization
GM_Kux  | 15 |UU P 177MR      |m^2/s           |K_11 element (U.point, X.dir) of GM-Redi tensor
GM_Kvy  | 15 |VV P 176MR      |m^2/s           |K_22 element (V.point, Y.dir) of GM-Redi tensor
GM_Kuz  | 15 |UU   179MR      |m^2/s           |K_13 element (U.point, Z.dir) of GM-Redi tensor
GM_Kvz  | 15 |VV   178MR      |m^2/s           |K_23 element (V.point, Z.dir) of GM-Redi tensor
GM_Kwx  | 15 |UM   181LR      |m^2/s           |K_31 element (W.point, X.dir) of GM-Redi tensor
GM_Kwy  | 15 |VM   180LR      |m^2/s           |K_32 element (W.point, Y.dir) of GM-Redi tensor
GM_Kwz  | 15 |WM P    LR      |m^2/s           |K_33 element (W.point, Z.dir) of GM-Redi tensor
GM_PsiX | 15 |UU   184LR      |m^2/s           |GM Bolus transport stream-function : X component
GM_PsiY | 15 |VV   183LR      |m^2/s           |GM Bolus transport stream-function : Y component
GM_KuzTz| 15 |UU   186MR      |degC.m^3/s      |Redi Off-diagonal Tempetature flux: X component
GM_KvzTz| 15 |VV   185MR      |degC.m^3/s      |Redi Off-diagonal Tempetature flux: Y component
\end{verbatim}
}

\subsubsection{Experiments and tutorials that use gmredi}
\label{sec:pkg:gmredi:experiments}

\begin{itemize}
\item{Global Ocean tutorial, in tutorial\_global\_oce\_latlon verification directory,
described in section \ref{sec:eg-global} }
\item{ Front Relax experiment, in front\_relax verification directory.}
\item{ Ideal 2D Ocean experiment, in ideal\_2D\_oce verification directory.}
\end{itemize}

% DO NOT EDIT HERE


1	dfer	1.16	\subsection{GMREDI: Gent-McWilliams/Redi SGS Eddy Parameterization}
2	edhill	1.8	\label{sec:pkg:gmredi}
3	edhill	1.7	\begin{rawhtml}
4			<!-- CMIREDIR:gmredi: -->
5			\end{rawhtml}
6	adcroft	1.1
7			There are two parts to the Redi/GM parameterization of geostrophic
8	jmc	1.17	eddies. The first, the Redi scheme \citep{re82}, aims to mix tracer properties
9			along isentropes (neutral surfaces) by means of a diffusion operator oriented
10			along the local isentropic surface.
11			The second part, GM \citep{gen-mcw:90,gen-eta:95}, adiabatically
12	adcroft	1.1	re-arranges tracers through an advective flux where the advecting flow
13	dfer	1.16	is a function of slope of the isentropic surfaces.
14	adcroft	1.1
15	dfer	1.16	The first GCM implementation of the Redi scheme was by \cite{Cox87} in the
16	adcroft	1.1	GFDL ocean circulation model. The original approach failed to
17			distinguish between isopycnals and surfaces of locally referenced
18			potential density (now called neutral surfaces) which are proper
19			isentropes for the ocean. As will be discussed later, it also appears
20			that the Cox implementation is susceptible to a computational mode.
21			Due to this mode, the Cox scheme requires a background lateral
22			diffusion to be present to conserve the integrity of the model fields.
23
24			The GM parameterization was then added to the GFDL code in the form of
25			a non-divergent bolus velocity. The method defines two
26			stream-functions expressed in terms of the isoneutral slopes subject
27			to the boundary condition of zero value on upper and lower
28			boundaries. The horizontal bolus velocities are then the vertical
29			derivative of these functions. Here in lies a problem highlighted by
30	dfer	1.16	\cite{gretal:98}: the bolus velocities involve multiple
31	adcroft	1.1	derivatives on the potential density field, which can consequently
32			give rise to noise. Griffies et al. point out that the GM bolus fluxes
33			can be identically written as a skew flux which involves fewer
34			differential operators. Further, combining the skew flux formulation
35			and Redi scheme, substantial cancellations take place to the point
36			that the horizontal fluxes are unmodified from the lateral diffusion
37			parameterization.
38
39	molod	1.10	\subsubsection{Redi scheme: Isopycnal diffusion}
40	adcroft	1.1
41			The Redi scheme diffuses tracers along isopycnals and introduces a
42			term in the tendency (rhs) of such a tracer (here $\tau$) of the form:
43			\begin{equation}
44			\bf{\nabla} \cdot \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
45			\end{equation}
46			where $\kappa_\rho$ is the along isopycnal diffusivity and
47			$\bf{K}_{Redi}$ is a rank 2 tensor that projects the gradient of
48			$\tau$ onto the isopycnal surface. The unapproximated projection tensor is:
49			\begin{equation}
50	jmc	1.17	\bf{K}_{Redi} = \frac{1}{1 + \|S\|^2} \left(
51	adcroft	1.1	\begin{array}{ccc}
52	jmc	1.17	1 + S_y^2& -S_x S_y & S_x \\
53			-S_x S_y & 1 + S_x^2 & S_y \\
54	adcroft	1.1	S_x & S_y & \|S\|^2 \\
55			\end{array}
56			\right)
57			\end{equation}
58			Here, $S_x = -\partial_x \sigma / \partial_z \sigma$ and $S_y =
59			-\partial_y \sigma / \partial_z \sigma$ are the components of the
60			isoneutral slope.
61
62			The first point to note is that a typical slope in the ocean interior
63			is small, say of the order $10^{-4}$. A maximum slope might be of
64			order $10^{-2}$ and only exceeds such in unstratified regions where
65			the slope is ill defined. It is therefore justifiable, and
66			customary, to make the small slope approximation, $\|S\| << 1$. The Redi
67			projection tensor then becomes:
68			\begin{equation}
69			\bf{K}_{Redi} = \left(
70			\begin{array}{ccc}
71			1 & 0 & S_x \\
72			0 & 1 & S_y \\
73			S_x & S_y & \|S\|^2 \\
74			\end{array}
75			\right)
76			\end{equation}
77
78
79	molod	1.10	\subsubsection{GM parameterization}
80	adcroft	1.1
81	dfer	1.16	The GM parameterization aims to represent the ``advective'' or
82	adcroft	1.1	``transport'' effect of geostrophic eddies by means of a ``bolus''
83	dfer	1.16	velocity, $\bf{u}^\star$. The divergence of this advective flux is added
84	adcroft	1.1	to the tracer tendency equation (on the rhs):
85			\begin{equation}
86	dfer	1.16	- \bf{\nabla} \cdot \tau \bf{u}^\star
87	adcroft	1.1	\end{equation}
88
89	dfer	1.16	The bolus velocity $\bf{u}^\star$ is defined as the rotational of a
90			streamfunction $\bf{F}^\star$=$(F_x^\star,F_y^\star,0)$:
91			\begin{equation}
92			\bf{u}^\star = \nabla \times \bf{F}^\star =
93			\left( \begin{array}{c}
94			- \partial_z F_y^\star \\
95			+ \partial_z F_x^\star \\
96			\partial_x F_y^\star - \partial_y F_x^\star
97			\end{array} \right),
98			\end{equation}
99			and thus is automatically non-divergent. In the GM parameterization, the streamfunction is
100			specified in terms of the isoneutral slopes $S_x$ and $S_y$:
101	adcroft	1.1	\begin{eqnarray}
102	dfer	1.16	F_x^\star & = & -\kappa_{GM} S_y \\
103			F_y^\star & = & \kappa_{GM} S_x
104	adcroft	1.1	\end{eqnarray}
105	dfer	1.16	with boundary conditions $F_x^\star=F_y^\star=0$ on upper and lower boundaries.
106			In the end, the bolus transport in the GM parameterization is given by:
107			\begin{equation}
108			\bf{u}^\star = \left(
109			\begin{array}{c}
110			u^\star \\
111			v^\star \\
112			w^\star
113			\end{array}
114			\right) = \left(
115			\begin{array}{c}
116			- \partial_z (\kappa_{GM} S_x) \\
117			- \partial_z (\kappa_{GM} S_y) \\
118			\partial_x (\kappa_{GM} S_x) + \partial_y (\kappa_{GM} S_y)
119			\end{array}
120			\right)
121			\end{equation}
122
123	adcroft	1.1	This is the form of the GM parameterization as applied by Donabasaglu,
124			1997, in MOM versions 1 and 2.
125
126	dfer	1.16	Note that in the MITgcm, the variables containing the GM bolus streamfunction are:
127			\begin{equation}
128			\left(
129			\begin{array}{c}
130			GM\_PsiX \\
131			GM\_PsiY
132			\end{array}
133			\right) = \left(
134			\begin{array}{c}
135			\kappa_{GM} S_x \\
136			\kappa_{GM} S_y
137			\end{array}
138			\right)= \left(
139			\begin{array}{c}
140			F_y^\star \\
141			-F_x^\star
142			\end{array}
143			\right).
144			\end{equation}
145
146	molod	1.10	\subsubsection{Griffies Skew Flux}
147	adcroft	1.1
148	dfer	1.16	\cite{gr:98} notes that the discretisation of bolus velocities involves
149	adcroft	1.1	multiple layers of differencing and interpolation that potentially
150			lead to noisy fields and computational modes. He pointed out that the
151			bolus flux can be re-written in terms of a non-divergent flux and a
152			skew-flux:
153			\begin{eqnarray*}
154	dfer	1.16	\bf{u}^\star \tau
155	adcroft	1.1	& = &
156			\left( \begin{array}{c}
157			- \partial_z ( \kappa_{GM} S_x ) \tau \\
158			- \partial_z ( \kappa_{GM} S_y ) \tau \\
159			(\partial_x \kappa_{GM} S_x + \partial_y \kappa_{GM} S_y)\tau
160			\end{array} \right)
161			\\
162			& = &
163			\left( \begin{array}{c}
164			- \partial_z ( \kappa_{GM} S_x \tau) \\
165			- \partial_z ( \kappa_{GM} S_y \tau) \\
166	dfer	1.16	\partial_x ( \kappa_{GM} S_x \tau) + \partial_y ( \kappa_{GM} S_y \tau)
167	adcroft	1.1	\end{array} \right)
168			+ \left( \begin{array}{c}
169			\kappa_{GM} S_x \partial_z \tau \\
170			\kappa_{GM} S_y \partial_z \tau \\
171	dfer	1.16	- \kappa_{GM} S_x \partial_x \tau - \kappa_{GM} S_y \partial_y \tau
172	adcroft	1.1	\end{array} \right)
173			\end{eqnarray*}
174			The first vector is non-divergent and thus has no effect on the tracer
175			field and can be dropped. The remaining flux can be written:
176			\begin{equation}
177	dfer	1.16	\bf{u}^\star \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau
178	adcroft	1.1	\end{equation}
179			where
180			\begin{equation}
181			\bf{K}_{GM} =
182			\left(
183			\begin{array}{ccc}
184			0 & 0 & -S_x \\
185			0 & 0 & -S_y \\
186			S_x & S_y & 0
187			\end{array}
188			\right)
189			\end{equation}
190			is an anti-symmetric tensor.
191
192			This formulation of the GM parameterization involves fewer derivatives
193			than the original and also involves only terms that already appear in
194			the Redi mixing scheme. Indeed, a somewhat fortunate cancellation
195			becomes apparent when we use the GM parameterization in conjunction
196			with the Redi isoneutral mixing scheme:
197			\begin{equation}
198			\kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
199	dfer	1.16	- u^\star \tau =
200	adcroft	1.1	( \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} ) \bf{\nabla} \tau
201			\end{equation}
202			In the instance that $\kappa_{GM} = \kappa_{\rho}$ then
203			\begin{equation}
204			\kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} =
205			\kappa_\rho
206			\left( \begin{array}{ccc}
207			1 & 0 & 0 \\
208			0 & 1 & 0 \\
209			2 S_x & 2 S_y & \|S\|^2
210			\end{array}
211			\right)
212			\end{equation}
213	cnh	1.3	which differs from the variable Laplacian diffusion tensor by only
214	adcroft	1.1	two non-zero elements in the $z$-row.
215
216	adcroft	1.2	\fbox{ \begin{minipage}{4.75in}
217			{\em S/R GMREDI\_CALC\_TENSOR} ({\em pkg/gmredi/gmredi\_calc\_tensor.F})
218
219			$\sigma_x$: {\bf SlopeX} (argument on entry)
220
221			$\sigma_y$: {\bf SlopeY} (argument on entry)
222
223			$\sigma_z$: {\bf SlopeY} (argument)
224
225			$S_x$: {\bf SlopeX} (argument on exit)
226
227			$S_y$: {\bf SlopeY} (argument on exit)
228
229			\end{minipage} }
230
231
232
233	molod	1.10	\subsubsection{Variable $\kappa_{GM}$}
234	adcroft	1.1
235	dfer	1.16	\cite{visbeck:97} suggest making the eddy coefficient,
236	adcroft	1.1	$\kappa_{GM}$, a function of the Eady growth rate,
237			$\|f\|/\sqrt{Ri}$. The formula involves a non-dimensional constant,
238			$\alpha$, and a length-scale $L$:
239			\begin{displaymath}
240			\kappa_{GM} = \alpha L^2 \overline{ \frac{\|f\|}{\sqrt{Ri}} }^z
241			\end{displaymath}
242			where the Eady growth rate has been depth averaged (indicated by the
243			over-line). A local Richardson number is defined $Ri = N^2 / (\partial
244			u/\partial z)^2$ which, when combined with thermal wind gives:
245			\begin{displaymath}
246			\frac{1}{Ri} = \frac{(\frac{\partial u}{\partial z})^2}{N^2} =
247			\frac{ ( \frac{g}{f \rho_o} \| {\bf \nabla} \sigma \| )^2 }{N^2} =
248			\frac{ M^4 }{ \|f\|^2 N^2 }
249			\end{displaymath}
250			where $M^2$ is defined $M^2 = \frac{g}{\rho_o} \|{\bf \nabla} \sigma\|$.
251			Substituting into the formula for $\kappa_{GM}$ gives:
252			\begin{displaymath}
253			\kappa_{GM} = \alpha L^2 \overline{ \frac{M^2}{N} }^z =
254			\alpha L^2 \overline{ \frac{M^2}{N^2} N }^z =
255			\alpha L^2 \overline{ \|S\| N }^z
256			\end{displaymath}
257
258
259	molod	1.10	\subsubsection{Tapering and stability}
260	adcroft	1.1
261			Experience with the GFDL model showed that the GM scheme has to be
262			matched to the convective parameterization. This was originally
263			expressed in connection with the introduction of the KPP boundary
264	dfer	1.16	layer scheme \citep{lar-eta:94} but in fact, as subsequent experience
265	adcroft	1.1	with the MIT model has found, is necessary for any convective
266			parameterization.
267
268			\fbox{ \begin{minipage}{4.75in}
269			{\em S/R GMREDI\_SLOPE\_LIMIT} ({\em
270			pkg/gmredi/gmredi\_slope\_limit.F})
271
272			$\sigma_x, s_x$: {\bf SlopeX} (argument)
273
274			$\sigma_y, s_y$: {\bf SlopeY} (argument)
275
276			$\sigma_z$: {\bf dSigmadRReal} (argument)
277
278			$z_\sigma^{*}$: {\bf dRdSigmaLtd} (argument)
279
280			\end{minipage} }
281
282	adcroft	1.2	\begin{figure}
283			\begin{center}
284	jmc	1.14	\resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/tapers.eps}}
285	adcroft	1.2	\end{center}
286	dfer	1.16	\caption{Taper functions used in GKW91 and DM95.}
287	adcroft	1.2	\label{fig:tapers}
288			\end{figure}
289
290			\begin{figure}
291			\begin{center}
292	jmc	1.14	\resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/effective_slopes.eps}}
293	adcroft	1.2	\end{center}
294			\caption{Effective slope as a function of ``true'' slope using Cox
295			slope clipping, GKW91 limiting and DM95 limiting.}
296			\label{fig:effective_slopes}
297			\end{figure}
298
299	adcroft	1.1
300	dfer	1.16	\subsubsection*{Slope clipping}
301	adcroft	1.1
302			Deep convection sites and the mixed layer are indicated by
303			homogenized, unstable or nearly unstable stratification. The slopes in
304			such regions can be either infinite, very large with a sign reversal
305			or simply very large. From a numerical point of view, large slopes
306			lead to large variations in the tensor elements (implying large bolus
307	cnh	1.3	flow) and can be numerically unstable. This was first recognized by
308	dfer	1.16	\cite{Cox87} who implemented ``slope clipping'' in the isopycnal mixing
309	adcroft	1.1	tensor. Here, the slope magnitude is simply restricted by an upper
310			limit:
311			\begin{eqnarray}
312			\|\nabla \sigma\| & = & \sqrt{ \sigma_x^2 + \sigma_y^2 } \\
313			S_{lim} & = & - \frac{\|\nabla \sigma\|}{ S_{max} }
314			\;\;\;\;\;\;\;\; \mbox{where $S_{max}$ is a parameter} \\
315			\sigma_z^\star & = & \min( \sigma_z , S_{lim} ) \\
316			{[s_x,s_y]} & = & - \frac{ [\sigma_x,\sigma_y] }{\sigma_z^\star}
317			\end{eqnarray}
318			Notice that this algorithm assumes stable stratification through the
319			``min'' function. In the case where the fluid is well stratified ($\sigma_z < S_{lim}$) then the slopes evaluate to:
320			\begin{equation}
321			{[s_x,s_y]} = - \frac{ [\sigma_x,\sigma_y] }{\sigma_z}
322			\end{equation}
323			while in the limited regions ($\sigma_z > S_{lim}$) the slopes become:
324			\begin{equation}
325			{[s_x,s_y]} = \frac{ [\sigma_x,\sigma_y] }{\|\nabla \sigma\|/S_{max}}
326			\end{equation}
327			so that the slope magnitude is limited $\sqrt{s_x^2 + s_y^2} =
328			S_{max}$.
329
330			The slope clipping scheme is activated in the model by setting {\bf
331			GM\_tap\-er\_scheme = 'clipping'} in {\em data.gmredi}.
332
333			Even using slope clipping, it is normally the case that the vertical
334			diffusion term (with coefficient $\kappa_\rho{\bf K}_{33} =
335			\kappa_\rho S_{max}^2$) is large and must be time-stepped using an
336			implicit procedure (see section on discretisation and code later).
337			Fig. \ref{fig-mixedlayer} shows the mixed layer depth resulting from
338			a) using the GM scheme with clipping and b) no GM scheme (horizontal
339			diffusion). The classic result of dramatically reduced mixed layers is
340			evident. Indeed, the deep convection sites to just one or two points
341			each and are much shallower than we might prefer. This, it turns out,
342	cnh	1.3	is due to the over zealous re-stratification due to the bolus transport
343	adcroft	1.1	parameterization. Limiting the slopes also breaks the adiabatic nature
344			of the GM/Redi parameterization, re-introducing diabatic fluxes in
345			regions where the limiting is in effect.
346
347	dfer	1.16	\subsubsection*{Tapering: Gerdes, Koberle and Willebrand, Clim. Dyn. 1991}
348	adcroft	1.1
349	dfer	1.16	The tapering scheme used in \cite{gkw:91}
350	adcroft	1.1	addressed two issues with the clipping method: the introduction of
351			large vertical fluxes in addition to convective adjustment fluxes is
352			avoided by tapering the GM/Redi slopes back to zero in
353			low-stratification regions; the adjustment of slopes is replaced by a
354			tapering of the entire GM/Redi tensor. This means the direction of
355			fluxes is unaffected as the amplitude is scaled.
356
357			The scheme inserts a tapering function, $f_1(S)$, in front of the
358			GM/Redi tensor:
359			\begin{equation}
360			f_1(S) = \min \left[ 1, \left( \frac{S_{max}}{\|S\|}\right)^2 \right]
361			\end{equation}
362			where $S_{max}$ is the maximum slope you want allowed. Where the
363			slopes, $\|S\|<S_{max}$ then $f_1(S) = 1$ and the tensor is un-tapered
364			but where $\|S\| \ge S_{max}$ then $f_1(S)$ scales down the tensor so
365			that the effective vertical diffusivity term $\kappa f_1(S) \|S\|^2 =
366			\kappa S_{max}^2$.
367
368	dfer	1.16	The GKW91 tapering scheme is activated in the model by setting {\bf
369	adcroft	1.1	GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.
370
371	dfer	1.16	\subsubsection*{Tapering: Danabasoglu and McWilliams, J. Clim. 1995}
372	adcroft	1.1
373	dfer	1.16	The tapering scheme used by \cite{dm:95} followed a similar procedure but used a different
374	adcroft	1.1	tapering function, $f_1(S)$:
375			\begin{equation}
376			f_1(S) = \frac{1}{2} \left( 1+\tanh \left[ \frac{S_c - \|S\|}{S_d} \right] \right)
377			\end{equation}
378			where $S_c = 0.004$ is a cut-off slope and $S_d=0.001$ is a scale over
379			which the slopes are smoothly tapered. Functionally, the operates in
380			the same way as the GKW91 scheme but has a substantially lower
381			cut-off, turning off the GM/Redi SGS parameterization for weaker
382			slopes.
383
384	dfer	1.16	The DM95 tapering scheme is activated in the model by setting {\bf
385	adcroft	1.1	GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.
386
387	dfer	1.16	\subsubsection*{Tapering: Large, Danabasoglu and Doney, JPO 1997}
388	adcroft	1.1
389	jmc	1.18	The tapering used in \cite{lar-eta:97} is based on the
390	adcroft	1.1	DM95 tapering scheme, but also tapers the scheme with an additional
391			function of height, $f_2(z)$, so that the GM/Redi SGS fluxes are
392			reduced near the surface:
393			\begin{equation}
394	dfer	1.16	f_2(z) = \frac{1}{2} \left( 1 + \sin(\pi \frac{z}{D} - \frac{\pi}{2})\right)
395	adcroft	1.1	\end{equation}
396			where $D = L_\rho \|S\|$ is a depth-scale and $L_\rho=c/f$ with
397			$c=2$~m~s$^{-1}$. This tapering with height was introduced to fix
398			some spurious interaction with the mixed-layer KPP parameterization.
399
400	dfer	1.16	The LDD97 tapering scheme is activated in the model by setting {\bf
401	adcroft	1.1	GM\_tap\-er\_scheme = 'ldd97'} in {\em data.gmredi}.
402
403
404
405	adcroft	1.2
406	adcroft	1.1	\begin{figure}
407	adcroft	1.4	\begin{center}
408	adcroft	1.1	%\includegraphics{mixedlayer-cox.eps}
409			%\includegraphics{mixedlayer-diff.eps}
410	adcroft	1.4	Figure missing.
411			\end{center}
412	adcroft	1.1	\caption{Mixed layer depth using GM parameterization with a) Cox slope
413			clipping and for comparison b) using horizontal constant diffusion.}
414	adcroft	1.4	\label{fig-mixedlayer}
415	adcroft	1.1	\end{figure}
416
417	molod	1.10	\subsubsection{Package Reference}
418	molod	1.11	\label{sec:pkg:gmredi:diagnostics}
419
420	edhill	1.12	{\footnotesize
421	molod	1.11	\begin{verbatim}
422			------------------------------------------------------------------------
423			<-Name->\|Levs\|<-parsing code->\|<-- Units -->\|<- Tile (max=80c)
424			------------------------------------------------------------------------
425			GM_VisbK\| 1 \|SM P M1 \|m^2/s \|Mixing coefficient from Visbeck etal parameterization
426			GM_Kux \| 15 \|UU P 177MR \|m^2/s \|K_11 element (U.point, X.dir) of GM-Redi tensor
427			GM_Kvy \| 15 \|VV P 176MR \|m^2/s \|K_22 element (V.point, Y.dir) of GM-Redi tensor
428			GM_Kuz \| 15 \|UU 179MR \|m^2/s \|K_13 element (U.point, Z.dir) of GM-Redi tensor
429			GM_Kvz \| 15 \|VV 178MR \|m^2/s \|K_23 element (V.point, Z.dir) of GM-Redi tensor
430			GM_Kwx \| 15 \|UM 181LR \|m^2/s \|K_31 element (W.point, X.dir) of GM-Redi tensor
431			GM_Kwy \| 15 \|VM 180LR \|m^2/s \|K_32 element (W.point, Y.dir) of GM-Redi tensor
432			GM_Kwz \| 15 \|WM P LR \|m^2/s \|K_33 element (W.point, Z.dir) of GM-Redi tensor
433			GM_PsiX \| 15 \|UU 184LR \|m^2/s \|GM Bolus transport stream-function : X component
434			GM_PsiY \| 15 \|VV 183LR \|m^2/s \|GM Bolus transport stream-function : Y component
435			GM_KuzTz\| 15 \|UU 186MR \|degC.m^3/s \|Redi Off-diagonal Tempetature flux: X component
436			GM_KvzTz\| 15 \|VV 185MR \|degC.m^3/s \|Redi Off-diagonal Tempetature flux: Y component
437			\end{verbatim}
438	edhill	1.12	}
439	molod	1.11
440	molod	1.13	\subsubsection{Experiments and tutorials that use gmredi}
441			\label{sec:pkg:gmredi:experiments}
442
443			\begin{itemize}
444			\item{Global Ocean tutorial, in tutorial\_global\_oce\_latlon verification directory,
445	jmc	1.15	described in section \ref{sec:eg-global} }
446	molod	1.13	\item{ Front Relax experiment, in front\_relax verification directory.}
447			\item{ Ideal 2D Ocean experiment, in ideal\_2D\_oce verification directory.}
448			\end{itemize}
449
450	cnh	1.6	% DO NOT EDIT HERE
451	adcroft	1.1
452
453