s_phys_pkgs/text/gmredi.tex

\section{Gent/McWiliams/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 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
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
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.

\subsection{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} = \left(
\begin{array}{ccc}
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{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}


\subsection{GM parameterization}

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

The bolus velocity is defined as:
\begin{eqnarray}
u^* & = & - \partial_z F_x \\
v^* & = & - \partial_z F_y \\
w^* & = & \partial_x F_x + \partial_y F_y
\end{eqnarray}
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 ($\partial_x u^* + \partial_y v^* = -
\partial_{xz} F_x - \partial_{yz} F_y = - \partial_z w^*$). $F_x$ and
$F_y$ are specified in terms of the isoneutral slopes $S_x$ and $S_y$:
\begin{eqnarray}
F_x & = & \kappa_{GM} S_x \\
F_y & = & \kappa_{GM} S_y
\end{eqnarray}
This is the form of the GM parameterization as applied by Donabasaglu,
1997, in MOM versions 1 and 2.

\subsection{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:
\begin{eqnarray*}
\bf{u}^* \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}^* \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^* \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} }


\subsection{Variable $\kappa_{GM}$}

Visbeck et al., 1996, 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}


\subsection{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 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{part6/tapers.eps}}
\end{center}
\caption{Taper functions used in GKW99 and DM95.}
\label{fig:tapers}
\end{figure}

\begin{figure}
\begin{center}
\resizebox{5.0in}{3.0in}{\includegraphics{part6/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
Cox, 1987, 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 Gerdes et al., 1999, (\cite{gkw:99})
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 GKW tapering scheme is activated in the model by setting {\bf
GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.

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

The tapering scheme used by Danabasoglu and McWilliams, 1995,
\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 DM tapering scheme is activated in the model by setting {\bf
GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.

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

The tapering used in Large et al., 1997, \cite{ldd: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(S) = \frac{1}{2} \left( 1 + \sin(\pi \frac{z}{D} - \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 LDD 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}

\subsection{Package Reference}
% DO NOT EDIT HERE


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