/[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.16 by dfer, Mon Jun 27 02:08:35 2011 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 along isentropes
9  (neutral surfaces) by means of a diffusion operator oriented along the  (neutral surfaces) by means of a diffusion operator oriented along the
10  local isentropic surface (Redi). The second part, adiabatically  local isentropic surface. The second part, GM \citep{gen-mcw:90,gen-eta:95}, adiabatically
11  re-arranges tracers through an advective flux where the advecting flow  re-arranges tracers through an advective flux where the advecting flow
12  is a function of slope of the isentropic surfaces (GM).  is a function of slope of the isentropic surfaces.
13    
14  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
15  GFDL ocean circulation model. The original approach failed to  GFDL ocean circulation model. The original approach failed to
16  distinguish between isopycnals and surfaces of locally referenced  distinguish between isopycnals and surfaces of locally referenced
17  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 26  stream-functions expressed in terms of t
26  to the boundary condition of zero value on upper and lower  to the boundary condition of zero value on upper and lower
27  boundaries. The horizontal bolus velocities are then the vertical  boundaries. The horizontal bolus velocities are then the vertical
28  derivative of these functions. Here in lies a problem highlighted by  derivative of these functions. Here in lies a problem highlighted by
29  Griffies et al., 1997: the bolus velocities involve multiple  \cite{gretal:98}: the bolus velocities involve multiple
30  derivatives on the potential density field, which can consequently  derivatives on the potential density field, which can consequently
31  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
32  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 35  and Redi scheme, substantial cancellatio
35  that the horizontal fluxes are unmodified from the lateral diffusion  that the horizontal fluxes are unmodified from the lateral diffusion
36  parameterization.  parameterization.
37    
38  \subsection{Redi scheme: Isopycnal diffusion}  \subsubsection{Redi scheme: Isopycnal diffusion}
39    
40  The Redi scheme diffuses tracers along isopycnals and introduces a  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:  term in the tendency (rhs) of such a tracer (here $\tau$) of the form:
# Line 71  S_x & S_y & |S|^2 \\ Line 75  S_x & S_y & |S|^2 \\
75  \end{equation}  \end{equation}
76    
77    
78  \subsection{GM parameterization}  \subsubsection{GM parameterization}
79    
80  The GM parameterization aims to parameterise the ``advective'' or  The GM parameterization aims to represent the ``advective'' or
81  ``transport'' effect of geostrophic eddies by means of a ``bolus''  ``transport'' effect of geostrophic eddies by means of a ``bolus''
82  velocity, $\bf{u}^*$. The divergence of this advective flux is added  velocity, $\bf{u}^\star$. The divergence of this advective flux is added
83  to the tracer tendency equation (on the rhs):  to the tracer tendency equation (on the rhs):
84  \begin{equation}  \begin{equation}
85  - \bf{\nabla} \cdot \tau \bf{u}^*  - \bf{\nabla} \cdot \tau \bf{u}^\star
86  \end{equation}  \end{equation}
87    
88  The bolus velocity is defined as:  The bolus velocity $\bf{u}^\star$ is defined as the rotational of a
89  \begin{eqnarray}  streamfunction $\bf{F}^\star$=$(F_x^\star,F_y^\star,0)$:
90  u^* & = & - \partial_z F_x \\  \begin{equation}
91  v^* & = & - \partial_z F_y \\  \bf{u}^\star = \nabla \times \bf{F}^\star =
92  w^* & = & \partial_x F_x + \partial_y F_y  \left( \begin{array}{c}
93  \end{eqnarray}  - \partial_z  F_y^\star \\
94  where $F_x$ and $F_y$ are stream-functions with boundary conditions  + \partial_z  F_x^\star \\
95  $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
96  bolus velocity in terms of these stream-functions is that they are  \end{array} \right),
97  automatically non-divergent ($\partial_x u^* + \partial_y v^* = -  \end{equation}
98  \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
99  $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$:
100  \begin{eqnarray}  \begin{eqnarray}
101  F_x & = & \kappa_{GM} S_x \\  F_x^\star & = & -\kappa_{GM} S_y \\
102  F_y & = & \kappa_{GM} S_y  F_y^\star & = &  \kappa_{GM} S_x
103  \end{eqnarray}  \end{eqnarray}
104    with boundary conditions $F_x^\star=F_y^\star=0$ on upper and lower boundaries.
105    In the end, the bolus transport in the GM parameterization is given by:
106    \begin{equation}
107    \bf{u}^\star = \left(
108    \begin{array}{c}
109    u^\star \\
110    v^\star \\
111    w^\star
112    \end{array}
113    \right) = \left(
114    \begin{array}{c}
115    - \partial_z (\kappa_{GM} S_x) \\
116    - \partial_z (\kappa_{GM} S_y) \\
117    \partial_x  (\kappa_{GM} S_x) + \partial_y (\kappa_{GM} S_y)
118    \end{array}
119    \right)
120    \end{equation}
121    
122  This is the form of the GM parameterization as applied by Donabasaglu,  This is the form of the GM parameterization as applied by Donabasaglu,
123  1997, in MOM versions 1 and 2.  1997, in MOM versions 1 and 2.
124    
125  \subsection{Griffies Skew Flux}  Note that in the MITgcm, the variables containing the GM bolus streamfunction are:
126    \begin{equation}
127    \left(
128    \begin{array}{c}
129    GM\_PsiX \\
130    GM\_PsiY
131    \end{array}
132    \right) = \left(
133    \begin{array}{c}
134    \kappa_{GM} S_x \\
135    \kappa_{GM} S_y
136    \end{array}
137    \right)= \left(
138    \begin{array}{c}
139    F_y^\star \\
140    -F_x^\star
141    \end{array}
142    \right).
143    \end{equation}
144      
145    \subsubsection{Griffies Skew Flux}
146    
147  Griffies notes that the discretisation of bolus velocities involves  \cite{gr:98} notes that the discretisation of bolus velocities involves
148  multiple layers of differencing and interpolation that potentially  multiple layers of differencing and interpolation that potentially
149  lead to noisy fields and computational modes. He pointed out that the  lead to noisy fields and computational modes. He pointed out that the
150  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
151  skew-flux:  skew-flux:
152  \begin{eqnarray*}  \begin{eqnarray*}
153  \bf{u}^* \tau  \bf{u}^\star \tau
154  & = &  & = &
155  \left( \begin{array}{c}  \left( \begin{array}{c}
156  - \partial_z ( \kappa_{GM} S_x ) \tau \\  - \partial_z ( \kappa_{GM} S_x ) \tau \\
# Line 120  skew-flux: Line 162  skew-flux:
162  \left( \begin{array}{c}  \left( \begin{array}{c}
163  - \partial_z ( \kappa_{GM} S_x \tau) \\  - \partial_z ( \kappa_{GM} S_x \tau) \\
164  - \partial_z ( \kappa_{GM} S_y \tau) \\  - \partial_z ( \kappa_{GM} S_y \tau) \\
165  \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)
166  \end{array} \right)  \end{array} \right)
167  + \left( \begin{array}{c}  + \left( \begin{array}{c}
168   \kappa_{GM} S_x \partial_z \tau \\   \kappa_{GM} S_x \partial_z \tau \\
169   \kappa_{GM} S_y \partial_z \tau \\   \kappa_{GM} S_y \partial_z \tau \\
170  - \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
171  \end{array} \right)  \end{array} \right)
172  \end{eqnarray*}  \end{eqnarray*}
173  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
174  field and can be dropped. The remaining flux can be written:  field and can be dropped. The remaining flux can be written:
175  \begin{equation}  \begin{equation}
176  \bf{u}^* \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau  \bf{u}^\star \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau
177  \end{equation}  \end{equation}
178  where  where
179  \begin{equation}  \begin{equation}
# Line 153  becomes apparent when we use the GM para Line 195  becomes apparent when we use the GM para
195  with the Redi isoneutral mixing scheme:  with the Redi isoneutral mixing scheme:
196  \begin{equation}  \begin{equation}
197  \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau  \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
198  - u^* \tau =  - u^\star \tau =
199  ( \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
200  \end{equation}  \end{equation}
201  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 209  In the instance that $\kappa_{GM} = \kap
209  \end{array}  \end{array}
210  \right)  \right)
211  \end{equation}  \end{equation}
212  which differs from the variable laplacian diffusion tensor by only  which differs from the variable Laplacian diffusion tensor by only
213  two non-zero elements in the $z$-row.  two non-zero elements in the $z$-row.
214    
215  \subsection{Variable $\kappa_{GM}$}  \fbox{ \begin{minipage}{4.75in}
216    {\em S/R GMREDI\_CALC\_TENSOR} ({\em pkg/gmredi/gmredi\_calc\_tensor.F})
217    
218    $\sigma_x$: {\bf SlopeX} (argument on entry)
219    
220    $\sigma_y$: {\bf SlopeY} (argument on entry)
221    
222    $\sigma_z$: {\bf SlopeY} (argument)
223    
224    $S_x$: {\bf SlopeX} (argument on exit)
225    
226    $S_y$: {\bf SlopeY} (argument on exit)
227    
228    \end{minipage} }
229    
230    
231    
232  Visbeck et al., 1996, suggest making the eddy coefficient,  \subsubsection{Variable $\kappa_{GM}$}
233    
234    \cite{visbeck:97} suggest making the eddy coefficient,
235  $\kappa_{GM}$, a function of the Eady growth rate,  $\kappa_{GM}$, a function of the Eady growth rate,
236  $|f|/\sqrt{Ri}$. The formula involves a non-dimensional constant,  $|f|/\sqrt{Ri}$. The formula involves a non-dimensional constant,
237  $\alpha$, and a length-scale $L$:  $\alpha$, and a length-scale $L$:
# Line 196  Substituting into the formula for $\kapp Line 255  Substituting into the formula for $\kapp
255  \end{displaymath}  \end{displaymath}
256    
257    
258  \subsection{Tapering and stability}  \subsubsection{Tapering and stability}
259    
260  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
261  matched to the convective parameterization. This was originally  matched to the convective parameterization. This was originally
262  expressed in connection with the introduction of the KPP boundary  expressed in connection with the introduction of the KPP boundary
263  layer scheme (Large et al., 97) but infact, as subsequent experience  layer scheme \citep{lar-eta:94} but in fact, as subsequent experience
264  with the MIT model has found, is necessary for any convective  with the MIT model has found, is necessary for any convective
265  parameterization.  parameterization.
266    
# Line 219  $z_\sigma^{*}$: {\bf dRdSigmaLtd} (argum Line 278  $z_\sigma^{*}$: {\bf dRdSigmaLtd} (argum
278    
279  \end{minipage} }  \end{minipage} }
280    
281    \begin{figure}
282    \begin{center}
283    \resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/tapers.eps}}
284    \end{center}
285    \caption{Taper functions used in GKW91 and DM95.}
286    \label{fig:tapers}
287    \end{figure}
288    
289  \subsubsection{Slope clipping}  \begin{figure}
290    \begin{center}
291    \resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/effective_slopes.eps}}
292    \end{center}
293    \caption{Effective slope as a function of ``true'' slope using Cox
294    slope clipping, GKW91 limiting and DM95 limiting.}
295    \label{fig:effective_slopes}
296    \end{figure}
297    
298    
299    \subsubsection*{Slope clipping}
300    
301  Deep convection sites and the mixed layer are indicated by  Deep convection sites and the mixed layer are indicated by
302  homogenized, unstable or nearly unstable stratification. The slopes in  homogenized, unstable or nearly unstable stratification. The slopes in
303  such regions can be either infinite, very large with a sign reversal  such regions can be either infinite, very large with a sign reversal
304  or simply very large. From a numerical point of view, large slopes  or simply very large. From a numerical point of view, large slopes
305  lead to large variations in the tensor elements (implying large bolus  lead to large variations in the tensor elements (implying large bolus
306  flow) and can be numerically unstable. This was first reognized by  flow) and can be numerically unstable. This was first recognized by
307  Cox, 1987, who implemented ``slope clipping'' in the isopycnal mixing  \cite{Cox87} who implemented ``slope clipping'' in the isopycnal mixing
308  tensor. Here, the slope magnitude is simply restricted by an upper  tensor. Here, the slope magnitude is simply restricted by an upper
309  limit:  limit:
310  \begin{eqnarray}  \begin{eqnarray}
# Line 262  a) using the GM scheme with clipping and Line 338  a) using the GM scheme with clipping and
338  diffusion). The classic result of dramatically reduced mixed layers is  diffusion). The classic result of dramatically reduced mixed layers is
339  evident. Indeed, the deep convection sites to just one or two points  evident. Indeed, the deep convection sites to just one or two points
340  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,
341  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
342  parameterization. Limiting the slopes also breaks the adiabatic nature  parameterization. Limiting the slopes also breaks the adiabatic nature
343  of the GM/Redi parameterization, re-introducing diabatic fluxes in  of the GM/Redi parameterization, re-introducing diabatic fluxes in
344  regions where the limiting is in effect.  regions where the limiting is in effect.
345    
346  \subsubsection{Tapering: Gerdes, Koberle and Willebrand, Clim. Dyn. 1991}  \subsubsection*{Tapering: Gerdes, Koberle and Willebrand, Clim. Dyn. 1991}
347    
348  The tapering scheme used in Gerdes et al., 1991, (\cite{gkw91})  The tapering scheme used in \cite{gkw:91}
349  addressed two issues with the clipping method: the introduction of  addressed two issues with the clipping method: the introduction of
350  large vertical fluxes in addition to convective adjustment fluxes is  large vertical fluxes in addition to convective adjustment fluxes is
351  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 364  but where $|S| \ge S_{max}$ then $f_1(S)
364  that the effective vertical diffusivity term $\kappa f_1(S) |S|^2 =  that the effective vertical diffusivity term $\kappa f_1(S) |S|^2 =
365  \kappa S_{max}^2$.  \kappa S_{max}^2$.
366    
367  The GKW tapering scheme is activated in the model by setting {\bf  The GKW91 tapering scheme is activated in the model by setting {\bf
368  GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.  GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.
369    
370  \subsection{Tapering: Danabasoglu and McWilliams, J. Clim. 1995}  \subsubsection*{Tapering: Danabasoglu and McWilliams, J. Clim. 1995}
371    
372  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  
373  tapering function, $f_1(S)$:  tapering function, $f_1(S)$:
374  \begin{equation}  \begin{equation}
375  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 380  the same way as the GKW91 scheme but has
380  cut-off, turning off the GM/Redi SGS parameterization for weaker  cut-off, turning off the GM/Redi SGS parameterization for weaker
381  slopes.  slopes.
382    
383  The DM tapering scheme is activated in the model by setting {\bf  The DM95 tapering scheme is activated in the model by setting {\bf
384  GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.  GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.
385    
386  \subsection{Tapering: Large, Danabasoglu and Doney, JPO 1997}  \subsubsection*{Tapering: Large, Danabasoglu and Doney, JPO 1997}
387    
388  The tapering used in Large et al., 1997, \cite{ldd97}, is based on the  The tapering used in \cite{ldd:97} is based on the
389  DM95 tapering scheme, but also tapers the scheme with an additional  DM95 tapering scheme, but also tapers the scheme with an additional
390  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
391  reduced near the surface:  reduced near the surface:
392  \begin{equation}  \begin{equation}
393  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)
394  \end{equation}  \end{equation}
395  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
396  $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
397  some spurious interaction with the mixed-layer KPP parameterization.  some spurious interaction with the mixed-layer KPP parameterization.
398    
399  The LDD tapering scheme is activated in the model by setting {\bf  The LDD97 tapering scheme is activated in the model by setting {\bf
400  GM\_tap\-er\_scheme = 'ldd97'} in {\em data.gmredi}.  GM\_tap\-er\_scheme = 'ldd97'} in {\em data.gmredi}.
401    
402    
403    
404    
405  \begin{figure}  \begin{figure}
406    \begin{center}
407  %\includegraphics{mixedlayer-cox.eps}  %\includegraphics{mixedlayer-cox.eps}
408  %\includegraphics{mixedlayer-diff.eps}  %\includegraphics{mixedlayer-diff.eps}
409    Figure missing.
410    \end{center}
411  \caption{Mixed layer depth using GM parameterization with a) Cox slope  \caption{Mixed layer depth using GM parameterization with a) Cox slope
412  clipping and for comparison b) using horizontal constant diffusion.}  clipping and for comparison b) using horizontal constant diffusion.}
413  \ref{fig-mixedlayer}  \label{fig-mixedlayer}
414  \end{figure}  \end{figure}
415    
416  \begin{figure}  \subsubsection{Package Reference}
417  %\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}$.  
418    
419  \subsection{subroutine RediTensor()}  {\footnotesize
   
 {\small  
420  \begin{verbatim}  \begin{verbatim}
421  subroutine RediTensor(Temp,Salt,Kredigm,K31,K32,K33, nIter,DumpFlag)  ------------------------------------------------------------------------
422                        |---in--| |-------out-------|  <-Name->|Levs|<-parsing code->|<--  Units   -->|<- Tile (max=80c)
423  ! Input  ------------------------------------------------------------------------
424  real Temp(Nx,Ny,Nz)     ! Potential temperature  GM_VisbK|  1 |SM P    M1      |m^2/s           |Mixing coefficient from Visbeck etal parameterization
425  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
426  ! Output  GM_Kvy  | 15 |VV P 176MR      |m^2/s           |K_22 element (V.point, Y.dir) of GM-Redi tensor
427  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
428  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
429  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
430  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
431  ! Auxiliary input  GM_Kwz  | 15 |WM P    LR      |m^2/s           |K_33 element (W.point, Z.dir) of GM-Redi tensor
432  integer nIter           ! interation/time-step number  GM_PsiX | 15 |UU   184LR      |m^2/s           |GM Bolus transport stream-function : X component
433  logical DumpFlag        ! flag to indicate routine should ``dump''  GM_PsiY | 15 |VV   183LR      |m^2/s           |GM Bolus transport stream-function : Y component
434    GM_KuzTz| 15 |UU   186MR      |degC.m^3/s      |Redi Off-diagonal Tempetature flux: X component
435    GM_KvzTz| 15 |VV   185MR      |degC.m^3/s      |Redi Off-diagonal Tempetature flux: Y component
436  \end{verbatim}  \end{verbatim}
437  }  }
438    
439  The subroutine {\bf RediTensor()} is called from {\bf model()} with  \subsubsection{Experiments and tutorials that use gmredi}
440  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}}  
441    
442  Note that $\sigma_z$ is the static stability because the potential  \begin{itemize}
443  densities are referenced to the same reference level ($W$-level).  \item{Global Ocean tutorial, in tutorial\_global\_oce\_latlon verification directory,
444    described in section \ref{sec:eg-global} }
445    \item{ Front Relax experiment, in front\_relax verification directory.}
446    \item{ Ideal 2D Ocean experiment, in ideal\_2D\_oce verification directory.}
447    \end{itemize}
448    
449  The next step calculates the three tensor components {\tt K13}, {\tt  % DO NOT EDIT HERE
 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*}  
450    
 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()}  
   
 {\bf inc\-tracer()} is called from {\bf model()} and has {\em four  
 new} arguments:  
 \begin{verbatim}  
 subroutine inc_tracer( ...,Kredigm,K31,K32,K33, ... )  
 real Kredigm(Nx,Ny,Nz)  ! Eddy coefficient  
 real K31(Nx,Ny,Nz)      ! (3,1) tensor coefficient  
 real K32(Nx,Ny,Nz)      ! (3,2) tensor coefficient  
 real K33(Nx,Ny,Nz)      ! (3,3) tensor coefficient  
 \end{verbatim}  
   
 Within the routine, the lateral fluxes, {\tt fluxWest} and {\tt  
 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*}  
   
 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.  
451    
452    

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

  ViewVC Help
Powered by ViewVC 1.1.22