/[MITgcm]/manual/s_phys_pkgs/text/gmredi.tex
ViewVC logotype

Diff of /manual/s_phys_pkgs/text/gmredi.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by adcroft, Thu Sep 27 19:43:36 2001 UTC revision 1.18 by jmc, Tue Mar 26 15:10:46 2013 UTC
# Line 1  Line 1 
1  \section{Gent/McWiliams/Redi SGS Eddy parameterization}  \subsection{GMREDI: Gent-McWilliams/Redi SGS Eddy Parameterization}
2    \label{sec:pkg:gmredi}
3    \begin{rawhtml}
4    <!-- CMIREDIR:gmredi: -->
5    \end{rawhtml}
6    
7  There are two parts to the Redi/GM parameterization of geostrophic  There are two parts to the Redi/GM parameterization of geostrophic
8  eddies. The first aims to mix tracer properties along isentropes  eddies. The first, the Redi scheme \citep{re82}, aims to mix tracer properties
9  (neutral surfaces) by means of a diffusion operator oriented along the  along isentropes (neutral surfaces) by means of a diffusion operator oriented
10  local isentropic surface (Redi). The second part, adiabatically  along the local isentropic surface.
11    The second part, GM \citep{gen-mcw:90,gen-eta:95}, adiabatically
12  re-arranges tracers through an advective flux where the advecting flow  re-arranges tracers through an advective flux where the advecting flow
13  is a function of slope of the isentropic surfaces (GM).  is a function of slope of the isentropic surfaces.
14    
15  The first GCM implementation of the Redi scheme was by Cox 1987 in the  The first GCM implementation of the Redi scheme was by \cite{Cox87} in the
16  GFDL ocean circulation model. The original approach failed to  GFDL ocean circulation model. The original approach failed to
17  distinguish between isopycnals and surfaces of locally referenced  distinguish between isopycnals and surfaces of locally referenced
18  potential density (now called neutral surfaces) which are proper  potential density (now called neutral surfaces) which are proper
# Line 22  stream-functions expressed in terms of t Line 27  stream-functions expressed in terms of t
27  to the boundary condition of zero value on upper and lower  to the boundary condition of zero value on upper and lower
28  boundaries. The horizontal bolus velocities are then the vertical  boundaries. The horizontal bolus velocities are then the vertical
29  derivative of these functions. Here in lies a problem highlighted by  derivative of these functions. Here in lies a problem highlighted by
30  Griffies et al., 1997: the bolus velocities involve multiple  \cite{gretal:98}: the bolus velocities involve multiple
31  derivatives on the potential density field, which can consequently  derivatives on the potential density field, which can consequently
32  give rise to noise. Griffies et al. point out that the GM bolus fluxes  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  can be identically written as a skew flux which involves fewer
# Line 31  and Redi scheme, substantial cancellatio Line 36  and Redi scheme, substantial cancellatio
36  that the horizontal fluxes are unmodified from the lateral diffusion  that the horizontal fluxes are unmodified from the lateral diffusion
37  parameterization.  parameterization.
38    
39  \subsection{Redi scheme: Isopycnal diffusion}  \subsubsection{Redi scheme: Isopycnal diffusion}
40    
41  The Redi scheme diffuses tracers along isopycnals and introduces a  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:  term in the tendency (rhs) of such a tracer (here $\tau$) of the form:
# Line 42  where $\kappa_\rho$ is the along isopycn Line 47  where $\kappa_\rho$ is the along isopycn
47  $\bf{K}_{Redi}$ is a rank 2 tensor that projects the gradient of  $\bf{K}_{Redi}$ is a rank 2 tensor that projects the gradient of
48  $\tau$ onto the isopycnal surface. The unapproximated projection tensor is:  $\tau$ onto the isopycnal surface. The unapproximated projection tensor is:
49  \begin{equation}  \begin{equation}
50  \bf{K}_{Redi} = \left(  \bf{K}_{Redi} = \frac{1}{1 + |S|^2} \left(
51  \begin{array}{ccc}  \begin{array}{ccc}
52  1 + S_x& S_x S_y & S_x \\  1 + S_y^2& -S_x S_y & S_x \\
53  S_x S_y  & 1 + S_y & S_y \\  -S_x S_y  & 1 + S_x^2 & S_y \\
54  S_x & S_y & |S|^2 \\  S_x & S_y & |S|^2 \\
55  \end{array}  \end{array}
56  \right)  \right)
# Line 71  S_x & S_y & |S|^2 \\ Line 76  S_x & S_y & |S|^2 \\
76  \end{equation}  \end{equation}
77    
78    
79  \subsection{GM parameterization}  \subsubsection{GM parameterization}
80    
81  The GM parameterization aims to parameterise the ``advective'' or  The GM parameterization aims to represent the ``advective'' or
82  ``transport'' effect of geostrophic eddies by means of a ``bolus''  ``transport'' effect of geostrophic eddies by means of a ``bolus''
83  velocity, $\bf{u}^*$. The divergence of this advective flux is added  velocity, $\bf{u}^\star$. The divergence of this advective flux is added
84  to the tracer tendency equation (on the rhs):  to the tracer tendency equation (on the rhs):
85  \begin{equation}  \begin{equation}
86  - \bf{\nabla} \cdot \tau \bf{u}^*  - \bf{\nabla} \cdot \tau \bf{u}^\star
87  \end{equation}  \end{equation}
88    
89  The bolus velocity is defined as:  The bolus velocity $\bf{u}^\star$ is defined as the rotational of a
90  \begin{eqnarray}  streamfunction $\bf{F}^\star$=$(F_x^\star,F_y^\star,0)$:
91  u^* & = & - \partial_z F_x \\  \begin{equation}
92  v^* & = & - \partial_z F_y \\  \bf{u}^\star = \nabla \times \bf{F}^\star =
93  w^* & = & \partial_x F_x + \partial_y F_y  \left( \begin{array}{c}
94  \end{eqnarray}  - \partial_z  F_y^\star \\
95  where $F_x$ and $F_y$ are stream-functions with boundary conditions  + \partial_z  F_x^\star \\
96  $F_x=F_y=0$ on upper and lower boundaries. The virtue of casting the  \partial_x F_y^\star - \partial_y F_x^\star
97  bolus velocity in terms of these stream-functions is that they are  \end{array} \right),
98  automatically non-divergent ($\partial_x u^* + \partial_y v^* = -  \end{equation}
99  \partial_{xz} F_x - \partial_{yz} F_y = - \partial_z w^*$). $F_x$ and  and thus is automatically non-divergent. In the GM parameterization, the streamfunction is
100  $F_y$ are specified in terms of the isoneutral slopes $S_x$ and $S_y$:  specified in terms of the isoneutral slopes $S_x$ and $S_y$:
101  \begin{eqnarray}  \begin{eqnarray}
102  F_x & = & \kappa_{GM} S_x \\  F_x^\star & = & -\kappa_{GM} S_y \\
103  F_y & = & \kappa_{GM} S_y  F_y^\star & = &  \kappa_{GM} S_x
104  \end{eqnarray}  \end{eqnarray}
105    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  This is the form of the GM parameterization as applied by Donabasaglu,  This is the form of the GM parameterization as applied by Donabasaglu,
124  1997, in MOM versions 1 and 2.  1997, in MOM versions 1 and 2.
125    
126  \subsection{Griffies Skew Flux}  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    \subsubsection{Griffies Skew Flux}
147    
148  Griffies notes that the discretisation of bolus velocities involves  \cite{gr:98} notes that the discretisation of bolus velocities involves
149  multiple layers of differencing and interpolation that potentially  multiple layers of differencing and interpolation that potentially
150  lead to noisy fields and computational modes. He pointed out that the  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  bolus flux can be re-written in terms of a non-divergent flux and a
152  skew-flux:  skew-flux:
153  \begin{eqnarray*}  \begin{eqnarray*}
154  \bf{u}^* \tau  \bf{u}^\star \tau
155  & = &  & = &
156  \left( \begin{array}{c}  \left( \begin{array}{c}
157  - \partial_z ( \kappa_{GM} S_x ) \tau \\  - \partial_z ( \kappa_{GM} S_x ) \tau \\
# Line 120  skew-flux: Line 163  skew-flux:
163  \left( \begin{array}{c}  \left( \begin{array}{c}
164  - \partial_z ( \kappa_{GM} S_x \tau) \\  - \partial_z ( \kappa_{GM} S_x \tau) \\
165  - \partial_z ( \kappa_{GM} S_y \tau) \\  - \partial_z ( \kappa_{GM} S_y \tau) \\
166  \partial_x ( \kappa_{GM} S_x \tau) + \partial_y ( \kappa_{GM} S_y) \tau)  \partial_x ( \kappa_{GM} S_x \tau) + \partial_y ( \kappa_{GM} S_y \tau)
167  \end{array} \right)  \end{array} \right)
168  + \left( \begin{array}{c}  + \left( \begin{array}{c}
169   \kappa_{GM} S_x \partial_z \tau \\   \kappa_{GM} S_x \partial_z \tau \\
170   \kappa_{GM} S_y \partial_z \tau \\   \kappa_{GM} S_y \partial_z \tau \\
171  - \kappa_{GM} S_x \partial_x \tau - \kappa_{GM} S_y) \partial_y \tau  - \kappa_{GM} S_x \partial_x \tau - \kappa_{GM} S_y \partial_y \tau
172  \end{array} \right)  \end{array} \right)
173  \end{eqnarray*}  \end{eqnarray*}
174  The first vector is non-divergent and thus has no effect on the tracer  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:  field and can be dropped. The remaining flux can be written:
176  \begin{equation}  \begin{equation}
177  \bf{u}^* \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau  \bf{u}^\star \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau
178  \end{equation}  \end{equation}
179  where  where
180  \begin{equation}  \begin{equation}
# Line 153  becomes apparent when we use the GM para Line 196  becomes apparent when we use the GM para
196  with the Redi isoneutral mixing scheme:  with the Redi isoneutral mixing scheme:
197  \begin{equation}  \begin{equation}
198  \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau  \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
199  - u^* \tau =  - u^\star \tau =
200  ( \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} ) \bf{\nabla} \tau  ( \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} ) \bf{\nabla} \tau
201  \end{equation}  \end{equation}
202  In the instance that $\kappa_{GM} = \kappa_{\rho}$ then  In the instance that $\kappa_{GM} = \kappa_{\rho}$ then
# Line 167  In the instance that $\kappa_{GM} = \kap Line 210  In the instance that $\kappa_{GM} = \kap
210  \end{array}  \end{array}
211  \right)  \right)
212  \end{equation}  \end{equation}
213  which differs from the variable laplacian diffusion tensor by only  which differs from the variable Laplacian diffusion tensor by only
214  two non-zero elements in the $z$-row.  two non-zero elements in the $z$-row.
215    
216  \subsection{Variable $\kappa_{GM}$}  \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  Visbeck et al., 1996, suggest making the eddy coefficient,  
233    \subsubsection{Variable $\kappa_{GM}$}
234    
235    \cite{visbeck:97} suggest making the eddy coefficient,
236  $\kappa_{GM}$, a function of the Eady growth rate,  $\kappa_{GM}$, a function of the Eady growth rate,
237  $|f|/\sqrt{Ri}$. The formula involves a non-dimensional constant,  $|f|/\sqrt{Ri}$. The formula involves a non-dimensional constant,
238  $\alpha$, and a length-scale $L$:  $\alpha$, and a length-scale $L$:
# Line 196  Substituting into the formula for $\kapp Line 256  Substituting into the formula for $\kapp
256  \end{displaymath}  \end{displaymath}
257    
258    
259  \subsection{Tapering and stability}  \subsubsection{Tapering and stability}
260    
261  Experience with the GFDL model showed that the GM scheme has to be  Experience with the GFDL model showed that the GM scheme has to be
262  matched to the convective parameterization. This was originally  matched to the convective parameterization. This was originally
263  expressed in connection with the introduction of the KPP boundary  expressed in connection with the introduction of the KPP boundary
264  layer scheme (Large et al., 97) but infact, as subsequent experience  layer scheme \citep{lar-eta:94} but in fact, as subsequent experience
265  with the MIT model has found, is necessary for any convective  with the MIT model has found, is necessary for any convective
266  parameterization.  parameterization.
267    
# Line 219  $z_\sigma^{*}$: {\bf dRdSigmaLtd} (argum Line 279  $z_\sigma^{*}$: {\bf dRdSigmaLtd} (argum
279    
280  \end{minipage} }  \end{minipage} }
281    
282    \begin{figure}
283    \begin{center}
284    \resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/tapers.eps}}
285    \end{center}
286    \caption{Taper functions used in GKW91 and DM95.}
287    \label{fig:tapers}
288    \end{figure}
289    
290    \begin{figure}
291    \begin{center}
292    \resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/effective_slopes.eps}}
293    \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  \subsubsection{Slope clipping}  
300    \subsubsection*{Slope clipping}
301    
302  Deep convection sites and the mixed layer are indicated by  Deep convection sites and the mixed layer are indicated by
303  homogenized, unstable or nearly unstable stratification. The slopes in  homogenized, unstable or nearly unstable stratification. The slopes in
304  such regions can be either infinite, very large with a sign reversal  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  or simply very large. From a numerical point of view, large slopes
306  lead to large variations in the tensor elements (implying large bolus  lead to large variations in the tensor elements (implying large bolus
307  flow) and can be numerically unstable. This was first reognized by  flow) and can be numerically unstable. This was first recognized by
308  Cox, 1987, who implemented ``slope clipping'' in the isopycnal mixing  \cite{Cox87} who implemented ``slope clipping'' in the isopycnal mixing
309  tensor. Here, the slope magnitude is simply restricted by an upper  tensor. Here, the slope magnitude is simply restricted by an upper
310  limit:  limit:
311  \begin{eqnarray}  \begin{eqnarray}
# Line 262  a) using the GM scheme with clipping and Line 339  a) using the GM scheme with clipping and
339  diffusion). The classic result of dramatically reduced mixed layers is  diffusion). The classic result of dramatically reduced mixed layers is
340  evident. Indeed, the deep convection sites to just one or two points  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,  each and are much shallower than we might prefer. This, it turns out,
342  is due to the over zealous restratification due to the bolus transport  is due to the over zealous re-stratification due to the bolus transport
343  parameterization. Limiting the slopes also breaks the adiabatic nature  parameterization. Limiting the slopes also breaks the adiabatic nature
344  of the GM/Redi parameterization, re-introducing diabatic fluxes in  of the GM/Redi parameterization, re-introducing diabatic fluxes in
345  regions where the limiting is in effect.  regions where the limiting is in effect.
346    
347  \subsubsection{Tapering: Gerdes, Koberle and Willebrand, Clim. Dyn. 1991}  \subsubsection*{Tapering: Gerdes, Koberle and Willebrand, Clim. Dyn. 1991}
348    
349  The tapering scheme used in Gerdes et al., 1991, (\cite{gkw91})  The tapering scheme used in \cite{gkw:91}
350  addressed two issues with the clipping method: the introduction of  addressed two issues with the clipping method: the introduction of
351  large vertical fluxes in addition to convective adjustment fluxes is  large vertical fluxes in addition to convective adjustment fluxes is
352  avoided by tapering the GM/Redi slopes back to zero in  avoided by tapering the GM/Redi slopes back to zero in
# Line 288  but where $|S| \ge S_{max}$ then $f_1(S) Line 365  but where $|S| \ge S_{max}$ then $f_1(S)
365  that the effective vertical diffusivity term $\kappa f_1(S) |S|^2 =  that the effective vertical diffusivity term $\kappa f_1(S) |S|^2 =
366  \kappa S_{max}^2$.  \kappa S_{max}^2$.
367    
368  The GKW tapering scheme is activated in the model by setting {\bf  The GKW91 tapering scheme is activated in the model by setting {\bf
369  GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.  GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.
370    
371  \subsection{Tapering: Danabasoglu and McWilliams, J. Clim. 1995}  \subsubsection*{Tapering: Danabasoglu and McWilliams, J. Clim. 1995}
372    
373  The tapering scheme used by Danabasoglu and McWilliams, 1995,  The tapering scheme used by \cite{dm:95} followed a similar procedure but used a different
 \cite{DM95}, followed a similar procedure but used a different  
374  tapering function, $f_1(S)$:  tapering function, $f_1(S)$:
375  \begin{equation}  \begin{equation}
376  f_1(S) = \frac{1}{2} \left( 1+\tanh \left[ \frac{S_c - |S|}{S_d} \right] \right)  f_1(S) = \frac{1}{2} \left( 1+\tanh \left[ \frac{S_c - |S|}{S_d} \right] \right)
# Line 305  the same way as the GKW91 scheme but has Line 381  the same way as the GKW91 scheme but has
381  cut-off, turning off the GM/Redi SGS parameterization for weaker  cut-off, turning off the GM/Redi SGS parameterization for weaker
382  slopes.  slopes.
383    
384  The DM tapering scheme is activated in the model by setting {\bf  The DM95 tapering scheme is activated in the model by setting {\bf
385  GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.  GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.
386    
387  \subsection{Tapering: Large, Danabasoglu and Doney, JPO 1997}  \subsubsection*{Tapering: Large, Danabasoglu and Doney, JPO 1997}
388    
389  The tapering used in Large et al., 1997, \cite{ldd97}, is based on the  The tapering used in \cite{lar-eta:97} is based on the
390  DM95 tapering scheme, but also tapers the scheme with an additional  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  function of height, $f_2(z)$, so that the GM/Redi SGS fluxes are
392  reduced near the surface:  reduced near the surface:
393  \begin{equation}  \begin{equation}
394  f_2(S) = \frac{1}{2} \left( 1 + \sin(\pi \frac{z}{D} - \pi/2)\right)  f_2(z) = \frac{1}{2} \left( 1 + \sin(\pi \frac{z}{D} - \frac{\pi}{2})\right)
395  \end{equation}  \end{equation}
396  where $D = L_\rho |S|$ is a depth-scale and $L_\rho=c/f$ with  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  $c=2$~m~s$^{-1}$.  This tapering with height was introduced to fix
398  some spurious interaction with the mixed-layer KPP parameterization.  some spurious interaction with the mixed-layer KPP parameterization.
399    
400  The LDD tapering scheme is activated in the model by setting {\bf  The LDD97 tapering scheme is activated in the model by setting {\bf
401  GM\_tap\-er\_scheme = 'ldd97'} in {\em data.gmredi}.  GM\_tap\-er\_scheme = 'ldd97'} in {\em data.gmredi}.
402    
403    
404    
405    
406  \begin{figure}  \begin{figure}
407    \begin{center}
408  %\includegraphics{mixedlayer-cox.eps}  %\includegraphics{mixedlayer-cox.eps}
409  %\includegraphics{mixedlayer-diff.eps}  %\includegraphics{mixedlayer-diff.eps}
410    Figure missing.
411    \end{center}
412  \caption{Mixed layer depth using GM parameterization with a) Cox slope  \caption{Mixed layer depth using GM parameterization with a) Cox slope
413  clipping and for comparison b) using horizontal constant diffusion.}  clipping and for comparison b) using horizontal constant diffusion.}
414  \ref{fig-mixedlayer}  \label{fig-mixedlayer}
415  \end{figure}  \end{figure}
416    
417  \begin{figure}  \subsubsection{Package Reference}
418  %\includegraphics{slopelimits.eps}  \label{sec:pkg:gmredi:diagnostics}
 \caption{Effective slope as a function of ``true'' slope using a) Cox  
 slope clipping, b) GKW91 limiting, c) DM95 limiting and d) LDD97  
 limiting.}  
 \end{figure}  
   
   
 %\begin{figure}  
 %\includegraphics{coxslope.eps}  
 %\includegraphics{gkw91slope.eps}  
 %\includegraphics{dm95slope.eps}  
 %\includegraphics{ldd97slope.eps}  
 %\caption{Effective slope magnitude at 100~m depth evaluated using a)  
 %Cox slope clipping, b) GKW91 limiting, c) DM95 limiting and d) LDD97  
 %limiting.}  
 %\end{figure}  
   
 \section{Discretisation and code}  
   
 This is the old documentation.....has to be brought upto date with MITgcm.  
   
   
 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 ({\bf gmredi\_calc\_tensor})  
 and 2) evaluation of the lateral and vertical fluxes due to gradients  
 along isopycnals or bolus transport ({\bf gmredi\_xtransport}, {\bf  
 gmredi\_ytransport} and {\bf 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, each row corresponding to a 'u', 'v' or 'w' point  
 on the model grid.  
   
 The code that implements the Redi/GM/Griffies schemes involves an  
 original core routine {\bf inc\_tracer()} that is used to calculate  
 the tendency in the tracers (namely, salt and potential temperature)  
 and a new routine {\bf RediTensor()} that calculates the tensor  
 components and $\kappa_{GM}$.  
   
 \subsection{subroutine RediTensor()}  
419    
420  {\small  {\footnotesize
421  \begin{verbatim}  \begin{verbatim}
422  subroutine RediTensor(Temp,Salt,Kredigm,K31,K32,K33, nIter,DumpFlag)  ------------------------------------------------------------------------
423                        |---in--| |-------out-------|  <-Name->|Levs|<-parsing code->|<--  Units   -->|<- Tile (max=80c)
424  ! Input  ------------------------------------------------------------------------
425  real Temp(Nx,Ny,Nz)     ! Potential temperature  GM_VisbK|  1 |SM P    M1      |m^2/s           |Mixing coefficient from Visbeck etal parameterization
426  real Salt(Nx,Ny,Nz)     ! Salinity  GM_Kux  | 15 |UU P 177MR      |m^2/s           |K_11 element (U.point, X.dir) of GM-Redi tensor
427  ! Output  GM_Kvy  | 15 |VV P 176MR      |m^2/s           |K_22 element (V.point, Y.dir) of GM-Redi tensor
428  real Kredigm(Nx,Ny,Nz)  ! Redi/GM eddy coefficient  GM_Kuz  | 15 |UU   179MR      |m^2/s           |K_13 element (U.point, Z.dir) of GM-Redi tensor
429  real K31(Nx,Ny,Nz)      ! Redi/GM (3,1) tensor component  GM_Kvz  | 15 |VV   178MR      |m^2/s           |K_23 element (V.point, Z.dir) of GM-Redi tensor
430  real K32(Nx,Ny,Nz)      ! Redi/GM (3,2) tensor component  GM_Kwx  | 15 |UM   181LR      |m^2/s           |K_31 element (W.point, X.dir) of GM-Redi tensor
431  real K33(Nx,Ny,Nz)      ! Redi/GM (3,3) tensor component  GM_Kwy  | 15 |VM   180LR      |m^2/s           |K_32 element (W.point, Y.dir) of GM-Redi tensor
432  ! Auxiliary input  GM_Kwz  | 15 |WM P    LR      |m^2/s           |K_33 element (W.point, Z.dir) of GM-Redi tensor
433  integer nIter           ! interation/time-step number  GM_PsiX | 15 |UU   184LR      |m^2/s           |GM Bolus transport stream-function : X component
434  logical DumpFlag        ! flag to indicate routine should ``dump''  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}  \end{verbatim}
438  }  }
439    
440  The subroutine {\bf RediTensor()} is called from {\bf model()} with  \subsubsection{Experiments and tutorials that use gmredi}
441  input arguments $T$ and $S$. It returns the 3D-arrays {\tt Kredigm},  \label{sec:pkg:gmredi:experiments}
 {\t K31}, {\tt K32} and {\tt K33} which represent $\kappa_{GM}$ (at  
 $T/S$ points) and the three components of the bottom row in the  
 Redi/GM tensor; $2 S_x$, $2 S_y$ and $|S|^2$ respectively, all at $W$  
 points.  
   
 The discretisations and algorithm within {\bf RediTensor()} are as  
 follows.  The routine first calculates the locally reference potential  
 density $\sigma_\theta$ from $T$ and $S$ and calculates the potential  
 density gradients in subroutine {\bf gradSigma()}:  
   
 \centerline{\begin{tabular}{ccl}  
 & & \\  
 Array & Grid-point & Definition \\  
 {\tt SigX} & U &  
 $\sigma_x = \frac{1}{\Delta x} \delta_x \sigma|_{z(k)}$  
 \\  
 {\tt SigY} & V &  
 $\sigma_y = \frac{1}{\Delta y} \delta_y \sigma|_{z(k)}$  
 \\  
 {\tt SigZ} & W &  
 $\sigma_z = \frac{1}{\Delta z}  
 [ \sigma|_{z(k)}(k-1/2) - \sigma|_{z(k)}(k+1/2) ]$  
 \\  
 \\  
 \end{tabular}}  
   
 Note that $\sigma_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 {\tt K13}, {\tt  
 K23} and {\tt K33} in subroutine {\bf KtensorWface()}. First, the  
 lateral gradients $\sigma_x$ and $\sigma_y$ are interpolated to the  
 $W$ points and stored in intermediate variables:  
 \begin{eqnarray*}  
 \mbox{\tt Sx} & = & \overline{ \overline{ \sigma_x }^x }^z \\  
 \mbox{\tt Sy} & = & \overline{ \overline{ \sigma_y }^y }^z  
 \end{eqnarray*}  
 Next, the magnitude of ${\bf \nabla}_z \sigma$ is stored in an intermediate  
 variable:  
 \begin{displaymath}  
 \mbox{\tt Sxy2} = \sqrt{ {\tt Sx}^2 + {\tt Sy}^2 }  
 \end{displaymath}  
 The stratification ($\sigma_z$) is ``checked'' such that the slope  
 vector has magnitude less than or equal to {\tt Smax} and stored in  
 an intermediate variable:  
 \begin{displaymath}  
 \mbox{\tt Sz} = \max ( \sigma_z , - \mbox{\tt Sxy2/Smax} )  
 \end{displaymath}  
 This guarantees stability and at the same time retains the lateral  
 orientation of the slope vector. The tensor components are then calculated:  
 \begin{eqnarray*}  
 \mbox{\tt K13} & = & -2 {\tt Sx/Sz} \\  
 \mbox{\tt K23} & = & -2 {\tt Sx/Sz} \\  
 \mbox{\tt K33} & = & ({\tt Sx/Sz})^2 + ({\tt Sy/Sz})^2  
 \end{eqnarray*}  
   
 Finally, {\tt Kredigm} ($\kappa_{GM}$) is calculated in subroutine  
 {\bf GMRediCoefficient()}. First, all the gradients are interpolated  
 to the $T/S$ points and stored in intermediate variables:  
 \begin{eqnarray*}  
 \mbox{\tt Sx} & = & \overline{ \sigma_x }^x \\  
 \mbox{\tt Sy} & = & \overline{ \sigma_y }^y \\  
 \mbox{\tt Sz} & = & \overline{ \sigma_z }^z  
 \end{eqnarray*}  
 Again, a nominal stratification is found by ``check'' the magnitude of  
 the slope vector but here is converted to a Brunt-Vasala frequency:  
 \begin{eqnarray*}  
 {\tt M2} & = & \sqrt{ {\tt Sx}^2 + {\tt Sy}^2} \\  
 {\tt N2} & = & - \frac{g}{\rho_o} \max ( {\tt Sz} , -{\tt M2 / Smax}  
 \end{eqnarray*}  
 The magnitude of the slope is then $|S| = {\tt M2}/{\tt N2}$. The Eady  
 growth rate is defined as $|f|/\sqrt(Ri) = |S| N$ and is calculated  
 as:  
 \begin{displaymath}  
 {\tt FrRi} = \frac{\tt M2}{\tt N2} ( - \frac{g}{\rho} {\tt Sz} )  
 \end{displaymath}  
 The Eady growth rate is then averaged over the upper layers (about  
 1100m) and $\kappa_{GM}$ specified from this 2D-variable:  
 \begin{displaymath}  
 {\tt Kredigm} = 0.02 * (200d3 **2) * {\tt FrRi}  
 \end{displaymath}  
   
 \subsection{subroutine inc\_tracer()}  
442    
443  {\bf inc\-tracer()} is called from {\bf model()} and has {\em four  \begin{itemize}
444  new} arguments:  \item{Global Ocean tutorial, in tutorial\_global\_oce\_latlon verification directory,
445  \begin{verbatim}  described in section \ref{sec:eg-global} }
446  subroutine inc_tracer( ...,Kredigm,K31,K32,K33, ... )  \item{ Front Relax experiment, in front\_relax verification directory.}
447  real Kredigm(Nx,Ny,Nz)  ! Eddy coefficient  \item{ Ideal 2D Ocean experiment, in ideal\_2D\_oce verification directory.}
448  real K31(Nx,Ny,Nz)      ! (3,1) tensor coefficient  \end{itemize}
 real K32(Nx,Ny,Nz)      ! (3,2) tensor coefficient  
 real K33(Nx,Ny,Nz)      ! (3,3) tensor coefficient  
 \end{verbatim}  
449    
450  Within the routine, the lateral fluxes, {\tt fluxWest} and {\tt  % DO NOT EDIT HERE
 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:  
 \begin{eqnarray*}  
 {\tt fluxWest}(\tau) & = & \ldots +  
 \overline{\tt Kredigm}^x \partial_x \tau \\  
 {\tt fluxSouth}(\tau) & = & \ldots +  
 \overline{\tt Kredigm}^y \partial_y \tau  
 \end{eqnarray*}  
451    
 The Redi/GM/Griffies scheme adds three terms to the vertical flux  
 ({\tt fluxUpper}) in the tracer equation. It is discretise simply:  
 \begin{displaymath}  
 {\tt fluxUpper}(\tau) = \ldots + \overline{\tt Kredigm}^z  
 \left(  
 {\tt K13} \overline{\partial_x \tau}^{xz} +  
 {\tt K23} \overline{\partial_y \tau}^{yz} +  
 {\tt K33} \partial_z \tau  
 \right)  
 \end{displaymath}  
 On boundaries, {\tt fluxUpper} is set to zero.  
452    
453    

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.18

  ViewVC Help
Powered by ViewVC 1.1.22