New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Chap_LDF.tex in branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters – NEMO

source: branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_LDF.tex @ 6391

Last change on this file since 6391 was 6391, checked in by acc, 8 years ago

Branch 2016/dev_r6325_SIMPLIF_1. Documentation changes for the Smagorinsky reimplementation (#1689. 2016WP-SIMPLIF-5). Note this documentation has been purposely submitted to the 2016WP-SIMPLIF-1 development branch at gm's request to avoid complications or loss when Chap_LDF.tex is rewritten on that branch. This submission affects Chapters/Chap_LDF.tex and Namelist/namdyn_ldf only.

File size: 32.4 KB
Line 
1
2% ================================================================
3% Chapter ———  Lateral Ocean Physics (LDF)
4% ================================================================
5\chapter{Lateral Ocean Physics (LDF)}
6\label{LDF}
7\minitoc
8
9
10\newpage
11$\ $\newline    % force a new ligne
12
13
14The lateral physics terms in the momentum and tracer equations have been
15described in \S\ref{PE_zdf} and their discrete formulation in \S\ref{TRA_ldf} 
16and \S\ref{DYN_ldf}). In this section we further discuss each lateral physics option.
17Choosing one lateral physics scheme means for the user defining,
18(1) the type of operator used (laplacian or bilaplacian operators, or no lateral mixing term) ;
19(2) the direction along which the lateral diffusive fluxes are evaluated (model level, geopotential or isopycnal surfaces) ; and
20(3) the space and time variations of the eddy coefficients.
21These three aspects of the lateral diffusion are set through namelist parameters
22(see the \textit{\ngn{nam\_traldf}} and \textit{\ngn{nam\_dynldf}} below).
23Note that this chapter describes the standard implementation of iso-neutral
24tracer mixing, and Griffies's implementation, which is used if
25\np{traldf\_grif}=true, is described in Appdx\ref{sec:triad}
26
27%-----------------------------------nam_traldf - nam_dynldf--------------------------------------------
28\namdisplay{namtra_ldf} 
29\namdisplay{namdyn_ldf} 
30%--------------------------------------------------------------------------------------------------------------
31
32
33% ================================================================
34% Direction of lateral Mixing
35% ================================================================
36\section  [Direction of Lateral Mixing (\textit{ldfslp})]
37      {Direction of Lateral Mixing (\mdl{ldfslp})}
38\label{LDF_slp}
39
40%%%
41\gmcomment{  we should emphasize here that the implementation is a rather old one.
42Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme. }
43
44A direction for lateral mixing has to be defined when the desired operator does
45not act along the model levels. This occurs when $(a)$ horizontal mixing is
46required on tracer or momentum (\np{ln\_traldf\_hor} or \np{ln\_dynldf\_hor})
47in $s$- or mixed $s$-$z$- coordinates, and $(b)$ isoneutral mixing is required
48whatever the vertical coordinate is. This direction of mixing is defined by its
49slopes in the \textbf{i}- and \textbf{j}-directions at the face of the cell of the
50quantity to be diffused. For a tracer, this leads to the following four slopes :
51$r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \eqref{Eq_tra_ldf_iso}), while
52for momentum the slopes are  $r_{1t}$, $r_{1uw}$, $r_{2f}$, $r_{2uw}$ for
53$u$ and  $r_{1f}$, $r_{1vw}$, $r_{2t}$, $r_{2vw}$ for $v$.
54
55%gm% add here afigure of the slope in i-direction
56
57\subsection{slopes for tracer geopotential mixing in the $s$-coordinate}
58
59In $s$-coordinates, geopotential mixing ($i.e.$ horizontal mixing) $r_1$ and
60$r_2$ are the slopes between the geopotential and computational surfaces.
61Their discrete formulation is found by locally solving \eqref{Eq_tra_ldf_iso} 
62when the diffusive fluxes in the three directions are set to zero and $T$ is
63assumed to be horizontally uniform, $i.e.$ a linear function of $z_T$, the
64depth of a $T$-point.
65%gm { Steven : My version is obviously wrong since I'm left with an arbitrary constant which is the local vertical temperature gradient}
66
67\begin{equation} \label{Eq_ldfslp_geo}
68\begin{aligned}
69 r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)}
70           \;\delta_{i+1/2}[z_t]
71      &\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t] \ \ \
72\\
73 r_{2v} &= \frac{e_{3v}}{\left( e_{2v}\;\overline{\overline{e_{3w}}}^{\,j+1/2,\,k} \right)} 
74           \;\delta_{j+1/2} [z_t]
75      &\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t] \ \ \
76\\
77 r_{1w} &= \frac{1}{e_{1w}}\;\overline{\overline{\delta_{i+1/2}[z_t]}}^{\,i,\,k+1/2}
78      &\approx \frac{1}{e_{1w}}\; \delta_{i+1/2}[z_{uw}]
79 \\
80 r_{2w} &= \frac{1}{e_{2w}}\;\overline{\overline{\delta_{j+1/2}[z_t]}}^{\,j,\,k+1/2}
81      &\approx \frac{1}{e_{2w}}\; \delta_{j+1/2}[z_{vw}]
82 \\
83\end{aligned}
84\end{equation}
85
86%gm%  caution I'm not sure the simplification was a good idea!
87
88These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}=True,
89and either \np{ln\_traldf\_hor}=True or \np{ln\_dynldf\_hor}=True.
90
91\subsection{Slopes for tracer iso-neutral mixing}\label{LDF_slp_iso}
92In iso-neutral mixing  $r_1$ and $r_2$ are the slopes between the iso-neutral
93and computational surfaces. Their formulation does not depend on the vertical
94coordinate used. Their discrete formulation is found using the fact that the
95diffusive fluxes of locally referenced potential density ($i.e.$ $in situ$ density)
96vanish. So, substituting $T$ by $\rho$ in \eqref{Eq_tra_ldf_iso} and setting the
97diffusive fluxes in the three directions to zero leads to the following definition for
98the neutral slopes:
99
100\begin{equation} \label{Eq_ldfslp_iso}
101\begin{split}
102 r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac{\delta_{i+1/2}[\rho]}
103                        {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,i+1/2,\,k}}
104\\
105 r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac{\delta_{j+1/2}\left[\rho \right]}
106                        {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,j+1/2,\,k}}
107\\
108 r_{1w} &= \frac{e_{3w}}{e_{1w}}\; 
109         \frac{\overline{\overline{\delta_{i+1/2}[\rho]}}^{\,i,\,k+1/2}}
110             {\delta_{k+1/2}[\rho]}
111\\
112 r_{2w} &= \frac{e_{3w}}{e_{2w}}\; 
113         \frac{\overline{\overline{\delta_{j+1/2}[\rho]}}^{\,j,\,k+1/2}}
114             {\delta_{k+1/2}[\rho]}
115\\
116\end{split}
117\end{equation}
118
119%gm% rewrite this as the explanation is not very clear !!!
120%In practice, \eqref{Eq_ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \eqref{Eq_ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth.
121
122%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{McDougall1987}, therefore in \eqref{Eq_ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
123
124%In the $z$-coordinate, the derivative of the  \eqref{Eq_ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation.
125
126As the mixing is performed along neutral surfaces, the gradient of $\rho$ in
127\eqref{Eq_ldfslp_iso} has to be evaluated at the same local pressure (which,
128in decibars, is approximated by the depth in meters in the model). Therefore
129\eqref{Eq_ldfslp_iso} cannot be used as such, but further transformation is
130needed depending on the vertical coordinate used:
131
132\begin{description}
133
134\item[$z$-coordinate with full step : ] in \eqref{Eq_ldfslp_iso} the densities
135appearing in the $i$ and $j$ derivatives  are taken at the same depth, thus
136the $in situ$ density can be used. This is not the case for the vertical
137derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$, where $N^2$ 
138is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following
139\citet{McDougall1987} (see \S\ref{TRA_bn2}).
140
141\item[$z$-coordinate with partial step : ] this case is identical to the full step
142case except that at partial step level, the \emph{horizontal} density gradient
143is evaluated as described in \S\ref{TRA_zpshde}.
144
145\item[$s$- or hybrid $s$-$z$- coordinate : ] in the current release of \NEMO,
146iso-neutral mixing is only employed for $s$-coordinates if the
147Griffies scheme is used (\np{traldf\_grif}=true; see Appdx \ref{sec:triad}).
148In other words, iso-neutral mixing will only be accurately represented with a
149linear equation of state (\np{nn\_eos}=1 or 2). In the case of a "true" equation
150of state, the evaluation of $i$ and $j$ derivatives in \eqref{Eq_ldfslp_iso} 
151will include a pressure dependent part, leading to the wrong evaluation of
152the neutral slopes.
153
154%gm%
155Note: The solution for $s$-coordinate passes trough the use of different
156(and better) expression for the constraint on iso-neutral fluxes. Following
157\citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral
158diffusive flux of locally referenced potential density, we stay in the $T$-$S$ 
159plane and consider the balance between the neutral direction diffusive fluxes
160of potential temperature and salinity:
161\begin{equation}
162\alpha \ \textbf{F}(T) = \beta \ \textbf{F}(S)
163\end{equation}
164%gm{  where vector F is ....}
165
166This constraint leads to the following definition for the slopes:
167
168\begin{equation} \label{Eq_ldfslp_iso2}
169\begin{split}
170 r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac
171      {\alpha_u \;\delta_{i+1/2}[T] - \beta_u \;\delta_{i+1/2}[S]}
172      {\alpha_u \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,i+1/2,\,k}
173       -\beta_\;\overline{\overline{\delta_{k+1/2}[S]}}^{\,i+1/2,\,k} }
174\\
175 r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac
176      {\alpha_v \;\delta_{j+1/2}[T] - \beta_v \;\delta_{j+1/2}[S]}
177      {\alpha_v \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,j+1/2,\,k}
178       -\beta_\;\overline{\overline{\delta_{k+1/2}[S]}}^{\,j+1/2,\,k} }
179\\
180 r_{1w} &= \frac{e_{3w}}{e_{1w}}\; \frac
181      {\alpha_w \;\overline{\overline{\delta_{i+1/2}[T]}}^{\,i,\,k+1/2}
182       -\beta_\;\overline{\overline{\delta_{i+1/2}[S]}}^{\,i,\,k+1/2} }
183      {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]}
184\\
185 r_{2w} &= \frac{e_{3w}}{e_{2w}}\; \frac
186      {\alpha_w \;\overline{\overline{\delta_{j+1/2}[T]}}^{\,j,\,k+1/2}
187       -\beta_\;\overline{\overline{\delta_{j+1/2}[S]}}^{\,j,\,k+1/2} }
188      {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]}
189\\
190\end{split}
191\end{equation}
192where $\alpha$ and $\beta$, the thermal expansion and saline contraction
193coefficients introduced in \S\ref{TRA_bn2}, have to be evaluated at the three
194velocity points. In order to save computation time, they should be approximated
195by the mean of their values at $T$-points (for example in the case of $\alpha$
196$\alpha_u=\overline{\alpha_T}^{i+1/2}$$\alpha_v=\overline{\alpha_T}^{j+1/2}$ 
197and $\alpha_w=\overline{\alpha_T}^{k+1/2}$).
198
199Note that such a formulation could be also used in the $z$-coordinate and
200$z$-coordinate with partial steps cases.
201
202\end{description}
203
204This implementation is a rather old one. It is similar to the one
205proposed by Cox [1987], except for the background horizontal
206diffusion. Indeed, the Cox implementation of isopycnal diffusion in
207GFDL-type models requires a minimum background horizontal diffusion
208for numerical stability reasons.  To overcome this problem, several
209techniques have been proposed in which the numerical schemes of the
210ocean model are modified \citep{Weaver_Eby_JPO97,
211  Griffies_al_JPO98}. Griffies's scheme is now available in \NEMO if
212\np{traldf\_grif\_iso} is set true; see Appdx \ref{sec:triad}. Here,
213another strategy is presented \citep{Lazar_PhD97}: a local
214filtering of the iso-neutral slopes (made on 9 grid-points) prevents
215the development of grid point noise generated by the iso-neutral
216diffusion operator (Fig.~\ref{Fig_LDF_ZDF1}). This allows an
217iso-neutral diffusion scheme without additional background horizontal
218mixing. This technique can be viewed as a diffusion operator that acts
219along large-scale (2~$\Delta$x) \gmcomment{2deltax doesnt seem very
220  large scale} iso-neutral surfaces. The diapycnal diffusion required
221for numerical stability is thus minimized and its net effect on the
222flow is quite small when compared to the effect of an horizontal
223background mixing.
224
225Nevertheless, this iso-neutral operator does not ensure that variance cannot increase,
226contrary to the \citet{Griffies_al_JPO98} operator which has that property.
227
228%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
229\begin{figure}[!ht]      \begin{center}
230\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_LDF_ZDF1.pdf}
231\caption {    \label{Fig_LDF_ZDF1}
232averaging procedure for isopycnal slope computation.}
233\end{center}    \end{figure}
234%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
235
236%There are three additional questions about the slope calculation.
237%First the expression for the rotation tensor has been obtain assuming the "small slope" approximation, so a bound has to be imposed on slopes.
238%Second, numerical stability issues also require a bound on slopes.
239%Third, the question of boundary condition specified on slopes...
240
241%from griffies: chapter 13.1....
242
243
244
245% In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04},
246% the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly
247% to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the
248% surface motivates this flattening of isopycnals near the surface).
249
250For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also
251be bounded by $1/100$ everywhere. This constraint is applied in a piecewise linear
252fashion, increasing from zero at the surface to $1/100$ at $70$ metres and thereafter
253decreasing to zero at the bottom of the ocean. (the fact that the eddies "feel" the
254surface motivates this flattening of isopycnals near the surface).
255
256%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
257\begin{figure}[!ht]     \begin{center}
258\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_eiv_slp.pdf}
259\caption {     \label{Fig_eiv_slp}
260Vertical profile of the slope used for lateral mixing in the mixed layer :
261\textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior,
262which has to be adjusted at the surface boundary (i.e. it must tend to zero at the
263surface since there is no mixing across the air-sea interface: wall boundary
264condition). Nevertheless, the profile between the surface zero value and the interior
265iso-neutral one is unknown, and especially the value at the base of the mixed layer ;
266\textit{(b)} profile of slope using a linear tapering of the slope near the surface and
267imposing a maximum slope of 1/100 ; \textit{(c)} profile of slope actually used in
268\NEMO: a linear decrease of the slope from zero at the surface to its ocean interior
269value computed just below the mixed layer. Note the huge change in the slope at the
270base of the mixed layer between  \textit{(b)}  and \textit{(c)}.}
271\end{center}   \end{figure}
272%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
273
274\colorbox{yellow}{add here a discussion about the flattening of the slopes, vs  tapering the coefficient.}
275
276\subsection{slopes for momentum iso-neutral mixing}
277
278The iso-neutral diffusion operator on momentum is the same as the one used on
279tracers but applied to each component of the velocity separately (see
280\eqref{Eq_dyn_ldf_iso} in section~\ref{DYN_ldf_iso}). The slopes between the
281surface along which the diffusion operator acts and the surface of computation
282($z$- or $s$-surfaces) are defined at $T$-, $f$-, and \textit{uw}- points for the
283$u$-component, and $T$-, $f$- and \textit{vw}- points for the $v$-component.
284They are computed from the slopes used for tracer diffusion, $i.e.$ 
285\eqref{Eq_ldfslp_geo} and \eqref{Eq_ldfslp_iso} :
286
287\begin{equation} \label{Eq_ldfslp_dyn}
288\begin{aligned}
289&r_{1t}\ \ = \overline{r_{1u}}^{\,i}       &&&    r_{1f}\ \ &= \overline{r_{1u}}^{\,i+1/2} \\
290&r_{2f} \ \ = \overline{r_{2v}}^{\,j+1/2} &&&   r_{2t}\ &= \overline{r_{2v}}^{\,j} \\
291&r_{1uw}  = \overline{r_{1w}}^{\,i+1/2} &&\ \ \text{and} \ \ &   r_{1vw}&= \overline{r_{1w}}^{\,j+1/2} \\
292&r_{2uw}= \overline{r_{2w}}^{\,j+1/2} &&&         r_{2vw}&= \overline{r_{2w}}^{\,j+1/2}\\
293\end{aligned}
294\end{equation}
295
296The major issue remaining is in the specification of the boundary conditions.
297The same boundary conditions are chosen as those used for lateral
298diffusion along model level surfaces, i.e. using the shear computed along
299the model levels and with no additional friction at the ocean bottom (see
300{\S\ref{LBC_coast}).
301
302
303% ================================================================
304% Lateral Mixing Operator
305% ================================================================
306\section [Lateral Mixing Operators (\textit{ldftra}, \textit{ldfdyn})]
307        {Lateral Mixing Operators (\mdl{traldf}, \mdl{traldf}) }
308\label{LDF_op}
309
310
311   
312% ================================================================
313% Lateral Mixing Coefficients
314% ================================================================
315\section [Lateral Mixing Coefficient (\textit{ldftra}, \textit{ldfdyn})]
316        {Lateral Mixing Coefficient (\mdl{ldftra}, \mdl{ldfdyn}) }
317\label{LDF_coef}
318
319Introducing a space variation in the lateral eddy mixing coefficients changes
320the model core memory requirement, adding up to four extra three-dimensional
321arrays for the geopotential or isopycnal second order operator applied to
322momentum. Six CPP keys control the space variation of eddy coefficients:
323three for momentum and three for tracer. The three choices allow:
324a space variation in the three space directions (\key{traldf\_c3d}\key{dynldf\_c3d}),
325in the horizontal plane (\key{traldf\_c2d}\key{dynldf\_c2d}),
326or in the vertical only (\key{traldf\_c1d}\key{dynldf\_c1d}).
327The default option is a constant value over the whole ocean on both momentum and tracers.
328   
329The number of additional arrays that have to be defined and the gridpoint
330position at which they are defined depend on both the space variation chosen
331and the type of operator used. The resulting eddy viscosity and diffusivity
332coefficients can be a function of more than one variable. Changes in the
333computer code when switching from one option to another have been
334minimized by introducing the eddy coefficients as statement functions
335(include file \hf{ldftra\_substitute} and \hf{ldfdyn\_substitute}). The functions
336are replaced by their actual meaning during the preprocessing step (CPP).
337The specification of the space variation of the coefficient is made in
338\mdl{ldftra} and \mdl{ldfdyn}, or more precisely in include files
339\textit{traldf\_cNd.h90} and \textit{dynldf\_cNd.h90}, with N=1, 2 or 3.
340The user can modify these include files as he/she wishes. The way the
341mixing coefficient are set in the reference version can be briefly described
342as follows:
343
344\subsubsection{Constant Mixing Coefficients (default option)}
345When none of the \textbf{key\_dynldf\_...} and \textbf{key\_traldf\_...} keys are
346defined, a constant value is used over the whole ocean for momentum and
347tracers, which is specified through the \np{rn\_ahm0} and \np{rn\_aht0} namelist
348parameters.
349
350\subsubsection{Vertically varying Mixing Coefficients (\key{traldf\_c1d} and \key{dynldf\_c1d})} 
351The 1D option is only available when using the $z$-coordinate with full step.
352Indeed in all the other types of vertical coordinate, the depth is a 3D function
353of (\textbf{i},\textbf{j},\textbf{k}) and therefore, introducing depth-dependent
354mixing coefficients will require 3D arrays. In the 1D option, a hyperbolic variation
355of the lateral mixing coefficient is introduced in which the surface value is
356\np{rn\_aht0} (\np{rn\_ahm0}), the bottom value is 1/4 of the surface value,
357and the transition takes place around z=300~m with a width of 300~m
358($i.e.$ both the depth and the width of the inflection point are set to 300~m).
359This profile is hard coded in file \hf{traldf\_c1d}, but can be easily modified by users.
360
361\subsubsection{Horizontally Varying Mixing Coefficients (\key{traldf\_c2d} and \key{dynldf\_c2d})}
362By default the horizontal variation of the eddy coefficient depends on the local mesh
363size and the type of operator used:
364\begin{equation} \label{Eq_title}
365  A_l = \left\{     
366   \begin{aligned}
367         & \frac{\max(e_1,e_2)}{e_{max}} A_o^l           & \text{for laplacian operator } \\
368         & \frac{\max(e_1,e_2)^{3}}{e_{max}^{3}} A_o^l          & \text{for bilaplacian operator } 
369   \end{aligned}    \right.
370\end{equation}
371where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked
372ocean domain, and $A_o^l$ is the \np{rn\_ahm0} (momentum) or \np{rn\_aht0} (tracer)
373namelist parameter. This variation is intended to reflect the lesser need for subgrid
374scale eddy mixing where the grid size is smaller in the domain. It was introduced in
375the context of the DYNAMO modelling project \citep{Willebrand_al_PO01}.
376Note that such a grid scale dependance of mixing coefficients significantly increase
377the range of stability of model configurations presenting large changes in grid pacing
378such as global ocean models. Indeed, in such a case, a constant mixing coefficient
379can lead to a blow up of the model due to large coefficient compare to the smallest
380grid size (see \S\ref{STP_forward_imp}), especially when using a bilaplacian operator.
381
382Other formulations can be introduced by the user for a given configuration.
383For example, in the ORCA2 global ocean model (see Configurations), the laplacian
384viscosity operator uses \np{rn\_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ 
385north and south and decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s
386at the equator \citep{Madec_al_JPO96, Delecluse_Madec_Bk00}. This modification
387can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}.
388Similar modified horizontal variations can be found with the Antarctic or Arctic
389sub-domain options of ORCA2 and ORCA05 (see \&namcfg namelist).
390
391\subsubsection{Space Varying Mixing Coefficients (\key{traldf\_c3d} and \key{dynldf\_c3d})}
392
393The 3D space variation of the mixing coefficient is simply the combination of the
3941D and 2D cases, $i.e.$ a hyperbolic tangent variation with depth associated with
395a grid size dependence of the magnitude of the coefficient.
396
397\subsubsection{Space and Time Varying Mixing Coefficients}
398
399There are no default specifications of space and time varying mixing coefficient.  One
400available case is specific to the ORCA2 and ORCA05 global ocean configurations. It
401provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both
402iso-neutral and eddy induced velocity (ORCA05) that depends on the local growth rate of
403baroclinic instability. This specification is actually used when an ORCA key
404and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined.
405
406\subsubsection{Smagorinsky viscosity (\np{nn\_ahm\_ijk\_t}=32)}
407
408Setting \np{nn\_ahm\_ijk\_t}=32 activates a 3D, time-varying viscosity that depends on the
409resolved motions. Three control parameters are described in the following sections which
410may be set in \textit{\ngn{namdyn\_ldf}}. Following \citep{Smagorinsky_93} and
411\citep{Griffies_Hallberg_MWR00} the viscosity coefficient is set proportional to a local
412deformation rate based on the horizontal shear and tension, namely:
413
414\begin{equation}
415A_{hm} = \left(\frac{{\sf C_{Smag}}}{\pi}\right)^2L^2\vert{D}\vert
416\end{equation}
417where $D$ is the deformation rate and $L$ is the local gridscale given by:
418\begin{equation}
419L = \frac{2{e_1} {e_2}}{\left ( {e_1} + {e_2} \right )}
420\end{equation}
421In orthogonal curvilinear coordinate, we have the tension and shearing strains ($D_T$ and $D_S$) given by
422\begin{equation}       
423\begin{split}   \label{EVP_strain}
424D_T &= \frac{e_2}{e_1} \,\frac{\partial}{\partial i} \left( \frac{u} {e_2} \right)
425     - \frac{e_1}{e_2} \,\frac{\partial}{\partial j} \left( \frac{v} {e_1} \right) \\
426D_S &= \frac{e_1}{e_2} \,\frac{\partial}{\partial j} \left( \frac{u} {e_1} \right)
427     + \frac{e_2}{e_1} \,\frac{\partial}{\partial i} \left( \frac{v} {e_2} \right) \\
428\end{split}   \end{equation}
429and the deformation rate, $D$, is given by:
430\begin{equation}       
431D = \sqrt{  {D_T}^2 + {D_S}^2 }
432\end{equation}
433On an Arakawa C-grid, $D_T$ is naturally defined at $t$-points and $D_S$ at $f$-points.
434Their discretizations are:
435\begin{equation} \label{VP_strain_discrete} \begin{split}
436%
437{D_T} &\equiv \frac{1}{e_{1t} e_{2t}} \left(  e_{2t}^2 \,\delta_i
438                                          \left[ \frac{u}{e_{2u}} \right] 
439                    - e_{1t}^2 \,\delta_j \left[ \frac{v}{e_{1v}} \right]   \right\\
440%
441{D_S} &\equiv \frac{1}{e_{1f} e_{2f}} \left(  e_{1f}^2 \,\delta_{j+1/2} 
442                                         \left[ \frac{u}{e_{1u}} \right] 
443                    - e_{2f}^2 \,\delta_{i+1/2} \left[ \frac{v}{e_{2v}} \right] \right
444\end{split} \end{equation}
445The deformation rate, $D$ needs to be defined at both $t$-points and $f$-points:
446\begin{equation} \begin{split}
447D_t &= \sqrt{  {D_T}^2  +   \overline{ \overline{ {D_S}^2 }}^{\,i,\,j} \,  }   \\
448D_f &= \sqrt{               \overline{ \overline{ {D_T}^2 }}^{\,i,\,j} + {D_S}^2   \, }
449\end{split} \end{equation}
450
451\citep{Griffies_Hallberg_MWR00} suggest values in the range 2.2 to 4.0 of the coefficient
452$\sf C_{Smag}$ for oceanic flows. This value is set via the \np{rn\_csmc} namelist
453parameter.  The same reference also derives the numerical stability criteria that suggest
454the calculated viscosity should be bounded according to the following:
455
456\begin{equation}
457{ \vert U \vert L\over 2} \leq A_{hm} \leq {L^2 \over 4 (2\Delta t)}
458\end{equation}
459which differ from \citep{Griffies_Hallberg_MWR00} only in the explicit acknowledgement of
460the time interval as $2\Delta t$ to reflect the LeapFrog scheme used in NEMO.  These
461bounds are imposed on the coefficients calculated at T and F points with the speed
462calculated using the appropriate average of the squares of the surrounding velocities.
463
464Some flexibility in these bounds may be appropriate in certain situations. For example,
465when modelling the equatorial region the presence of a strong equatorial current may
466impose a higher lower bound than is desirable or strictly necessary. Flexability is
467provided via two extra namelist parameters which can be used to adjust either or both
468bounds by constant factors.  These parameters are: \np{rn\_minfac} and \np{rn\_maxfac}.
469Their default values are 1.0 which provide the exact bounds specified above but values
470other than 1.0 will adjust the bounds according to:
471
472\begin{equation}
473\np{rn\_minfac} \times { \vert U \vert L\over 2} \leq A_{hm} \leq 
474                                   {L^2 \over 8 \Delta t} \times \np{rn\_maxfac}
475\end{equation}
476
477\bigskip 
478
479When $ln\_dynldf\_bilap = .true.$, a biharmonic version of the Smagorinsky viscosity is
480also available which sets a coefficient for the biharmonic viscosity ( following
481\citep{Griffies_Hallberg_MWR00}) as:
482
483\begin{equation}\begin{split}
484B_{hm} &= \left(\frac{{\sf C_{Smag}}}{\pi}\right)^2 {L^4\over 8}\vert{D}\vert \\
485             &= A_{hm} {L^2\over 8}
486\end{split}\end{equation}
487where $A_{hm}$ is the laplacian form of the Smagorinsky viscosity.  This coefficient is
488computed from the laplacian values with bounds already applied so the effective bounds on
489$B_{hm}$ would be:
490\begin{equation}
491\np{rn\_minfac} \times { \vert U \vert L^3\over 16} \leq B_{hm} \leq 
492                                   {L^4 \over 32 \Delta t} \times \np{rn\_maxfac}
493\end{equation}
494but this lower limit does not reflect the behaviour of the 4th order advection schemes
495available in NEMO that have a diffusive component with a biharmonic operator whose eddy
496coefficient is equivalent to:
497\begin{equation}
498{ \vert U \vert L^3\over 12} 
499\end{equation}
500(see equation \eqref{Eq_traadv_ubs2} and surrounding discussion). A safer lower bound for
501the biharmonic Smagorinsky eddy coefficient in NEMO is therefore implemented by
502introducing a hard-coded multiplier of $4/3$ to the lower bound such that the bounds on
503the biharmonic coefficient are:
504\begin{equation}
505\np{rn\_minfac} \times { \vert U \vert L^3\over 12} \leq B_{hm} \leq 
506                                   {L^4 \over 32 \Delta t} \times \np{rn\_maxfac}
507\end{equation} 
508Note that the bilaplacian operator is implemented as a re-entrant laplacian
509which means the actual values of the coefficient returned by \np{ldf\_dyn} in the
510biharmonic case are $\sqrt{B_{hm}}$.
511
512$\ $\newline    % force a new ligne
513
514The following points are relevant when the eddy coefficient varies spatially:
515
516(1) the momentum diffusion operator acting along model level surfaces is
517written in terms of curl and divergent components of the horizontal current
518(see \S\ref{PE_ldf}). Although the eddy coefficient could be set to different values
519in these two terms, this option is not currently available.
520
521(2) with an horizontally varying viscosity, the quadratic integral constraints
522on enstrophy and on the square of the horizontal divergence for operators
523acting along model-surfaces are no longer satisfied
524(Appendix~\ref{Apdx_dynldf_properties}).
525
526(3) for isopycnal diffusion on momentum or tracers, an additional purely
527horizontal background diffusion with uniform coefficient can be added by
528setting a non zero value of \np{rn\_ahmb0} or \np{rn\_ahtb0}, a background horizontal
529eddy viscosity or diffusivity coefficient (namelist parameters whose default
530values are $0$). However, the technique used to compute the isopycnal
531slopes is intended to get rid of such a background diffusion, since it introduces
532spurious diapycnal diffusion (see {\S\ref{LDF_slp}).
533
534(4) when an eddy induced advection term is used (\key{traldf\_eiv}), $A^{eiv}$,
535the eddy induced coefficient has to be defined. Its space variations are controlled
536by the same CPP variable as for the eddy diffusivity coefficient ($i.e.$ 
537\textbf{key\_traldf\_cNd}).
538
539(5) the eddy coefficient associated with a biharmonic operator must be set to a \emph{negative} value.
540
541(6) it is not possible to use both the laplacian and biharmonic operators concurrently.
542
543(7) it is possible to run without explicit lateral diffusion on momentum (\np{ln\_dynldf\_lap} =
544\np{ln\_dynldf\_bilap} = false). This is recommended when using the UBS advection
545scheme on momentum (\np{ln\_dynadv\_ubs} = true, see \ref{DYN_adv_ubs})
546and can be useful for testing purposes.
547
548% ================================================================
549% Eddy Induced Mixing
550% ================================================================
551\section  [Eddy Induced Velocity (\textit{traadv\_eiv}, \textit{ldfeiv})]
552      {Eddy Induced Velocity (\mdl{traadv\_eiv}, \mdl{ldfeiv})}
553\label{LDF_eiv}
554
555%%gm  from Triad appendix  : to be incorporated....
556\gmcomment{
557Values of iso-neutral diffusivity and GM coefficient are set as
558described in \S\ref{LDF_coef}. If none of the keys \key{traldf\_cNd},
559N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and
560GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and
561\np{rn\_aeiv\_0}. If 2D-varying coefficients are set with
562\key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal
563scale factor according to \eqref{Eq_title} \footnote{Except in global ORCA
564  $0.5^{\circ}$ runs with \key{traldf\_eiv}, where
565  $A_l$ is set like $A_e$ but with a minimum vale of
566  $100\;\mathrm{m}^2\;\mathrm{s}^{-1}$}. In idealised setups with
567\key{traldf\_c2d}, $A_e$ is reduced similarly, but if \key{traldf\_eiv}
568is set in the global configurations with \key{traldf\_c2d}, a horizontally varying $A_e$ is
569instead set from the Held-Larichev parameterisation\footnote{In this
570  case, $A_e$ at low latitudes $|\theta|<20^{\circ}$ is further
571  reduced by a factor $|f/f_{20}|$, where $f_{20}$ is the value of $f$
572  at $20^{\circ}$~N} (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored
573unless it is zero.
574}
575
576When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined),
577an eddy induced tracer advection term is added, the formulation of which
578depends on the slopes of iso-neutral surfaces. Contrary to the case of iso-neutral
579mixing, the slopes used here are referenced to the geopotential surfaces, $i.e.$ 
580\eqref{Eq_ldfslp_geo} is used in $z$-coordinates, and the sum \eqref{Eq_ldfslp_geo} 
581+ \eqref{Eq_ldfslp_iso} in $s$-coordinates. The eddy induced velocity is given by:
582\begin{equation} \label{Eq_ldfeiv}
583\begin{split}
584 u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\
585v^* & = \frac{1}{e_{1u}e_{3v}}\; \delta_k \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right]\\
586w^* & = \frac{1}{e_{1w}e_{2w}}\; \left\{ \delta_i \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right] + \delta_j \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right] \right\} \\
587\end{split}
588\end{equation}
589where $A^{eiv}$ is the eddy induced velocity coefficient whose value is set
590through \np{rn\_aeiv}, a \textit{nam\_traldf} namelist parameter.
591The three components of the eddy induced velocity are computed and add
592to the eulerian velocity in \mdl{traadv\_eiv}. This has been preferred to a
593separate computation of the advective trends associated with the eiv velocity,
594since it allows us to take advantage of all the advection schemes offered for
595the tracers (see \S\ref{TRA_adv}) and not just the $2^{nd}$ order advection
596scheme as in previous releases of OPA \citep{Madec1998}. This is particularly
597useful for passive tracers where \emph{positivity} of the advection scheme is
598of paramount importance.
599
600At the surface, lateral and bottom boundaries, the eddy induced velocity,
601and thus the advective eddy fluxes of heat and salt, are set to zero.
602
603
604
605
606
Note: See TracBrowser for help on using the repository browser.