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 NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex @ 11561

Last change on this file since 11561 was 11561, checked in by nicolasmartin, 3 years ago

Apply two thirds ratio on figures width as default

File size: 32.3 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4
5% ================================================================
6% Chapter Lateral Ocean Physics (LDF)
7% ================================================================
8\chapter{Lateral Ocean Physics (LDF)}
9\label{chap:LDF}
10
11\chaptertoc
12
13\newpage
14
15The lateral physics terms in the momentum and tracer equations have been described in \autoref{eq:MB_zdf} and
16their discrete formulation in \autoref{sec:TRA_ldf} and \autoref{sec:DYN_ldf}).
17In this section we further discuss each lateral physics option.
18Choosing one lateral physics scheme means for the user defining,
19(1) the type of operator used (laplacian or bilaplacian operators, or no lateral mixing term);
20(2) the direction along which the lateral diffusive fluxes are evaluated
21(model level, geopotential or isopycnal surfaces); and
22(3) the space and time variations of the eddy coefficients.
23These three aspects of the lateral diffusion are set through namelist parameters
24(see the \nam{tra\_ldf} and \nam{dyn\_ldf} below).
25Note that this chapter describes the standard implementation of iso-neutral tracer mixing.
26Griffies's implementation, which is used if \np{ln\_traldf\_triad}\forcode{=.true.},
27is described in \autoref{apdx:TRIADS}
28
29%-----------------------------------namtra_ldf - namdyn_ldf--------------------------------------------
30
31%--------------------------------------------------------------------------------------------------------------
32
33% ================================================================
34% Lateral Mixing Operator
35% ================================================================
36\section[Lateral mixing operators]
37{Lateral mixing operators}
38\label{sec:LDF_op}
39We remind here the different lateral mixing operators that can be used. Further details can be found in \autoref{subsec:TRA_ldf_op} and  \autoref{sec:DYN_ldf}.
40
41\subsection[No lateral mixing (\forcode{ln_traldf_OFF}, \forcode{ln_dynldf_OFF})]
42{No lateral mixing (\protect\np{ln\_traldf\_OFF}, \protect\np{ln\_dynldf\_OFF})}
43
44It is possible to run without explicit lateral diffusion on tracers (\protect\np{ln\_traldf\_OFF}\forcode{=.true.}) and/or
45momentum (\protect\np{ln\_dynldf\_OFF}\forcode{=.true.}). The latter option is even recommended if using the
46UBS advection scheme on momentum (\np{ln\_dynadv\_ubs}\forcode{=.true.},
47see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes.
48
49\subsection[Laplacian mixing (\forcode{ln_traldf_lap}, \forcode{ln_dynldf_lap})]
50{Laplacian mixing (\protect\np{ln\_traldf\_lap}, \protect\np{ln\_dynldf\_lap})}
51Setting \protect\np{ln\_traldf\_lap}\forcode{=.true.} and/or \protect\np{ln\_dynldf\_lap}\forcode{=.true.} enables
52a second order diffusion on tracers and momentum respectively. Note that in \NEMO\ 4, one can not combine
53Laplacian and Bilaplacian operators for the same variable.
54
55\subsection[Bilaplacian mixing (\forcode{ln_traldf_blp}, \forcode{ln_dynldf_blp})]
56{Bilaplacian mixing (\protect\np{ln\_traldf\_blp}, \protect\np{ln\_dynldf\_blp})}
57Setting \protect\np{ln\_traldf\_blp}\forcode{=.true.} and/or \protect\np{ln\_dynldf\_blp}\forcode{=.true.} enables
58a fourth order diffusion on tracers and momentum respectively. It is implemented by calling the above Laplacian operator twice.
59We stress again that from \NEMO\ 4, the simultaneous use Laplacian and Bilaplacian operators is not allowed.
60
61% ================================================================
62% Direction of lateral Mixing
63% ================================================================
64\section[Direction of lateral mixing (\textit{ldfslp.F90})]
65{Direction of lateral mixing (\protect\mdl{ldfslp})}
66\label{sec:LDF_slp}
67
68%%%
69\gmcomment{
70  we should emphasize here that the implementation is a rather old one.
71  Better work can be achieved by using \citet{griffies.gnanadesikan.ea_JPO98, griffies_bk04} iso-neutral scheme.
72}
73
74A direction for lateral mixing has to be defined when the desired operator does not act along the model levels.
75This occurs when $(a)$ horizontal mixing is required on tracer or momentum
76(\np{ln\_traldf\_hor} or \np{ln\_dynldf\_hor}) in $s$- or mixed $s$-$z$- coordinates,
77and $(b)$ isoneutral mixing is required whatever the vertical coordinate is.
78This direction of mixing is defined by its slopes in the \textbf{i}- and \textbf{j}-directions at the face of
79the cell of the quantity to be diffused.
80For a tracer, this leads to the following four slopes:
81$r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \autoref{eq:TRA_ldf_iso}),
82while for momentum the slopes are  $r_{1t}$, $r_{1uw}$, $r_{2f}$, $r_{2uw}$ for $u$ and
83$r_{1f}$, $r_{1vw}$, $r_{2t}$, $r_{2vw}$ for $v$.
84
85%gm% add here afigure of the slope in i-direction
86
87\subsection{Slopes for tracer geopotential mixing in the $s$-coordinate}
88
89In $s$-coordinates, geopotential mixing (\ie\ horizontal mixing) $r_1$ and $r_2$ are the slopes between
90the geopotential and computational surfaces.
91Their discrete formulation is found by locally solving \autoref{eq:TRA_ldf_iso} when
92the diffusive fluxes in the three directions are set to zero and $T$ is assumed to be horizontally uniform,
93\ie\ a linear function of $z_T$, the depth of a $T$-point.
94%gm { Steven : My version is obviously wrong since I'm left with an arbitrary constant which is the local vertical temperature gradient}
95
96\begin{equation}
97  \label{eq:LDF_slp_geo}
98  \begin{aligned}
99    r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)}
100    \;\delta_{i+1/2}[z_t]
101    &\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t] \ \ \ \\
102    r_{2v} &= \frac{e_{3v}}{\left( e_{2v}\;\overline{\overline{e_{3w}}}^{\,j+1/2,\,k} \right)}
103    \;\delta_{j+1/2} [z_t]
104    &\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t] \ \ \ \\
105    r_{1w} &= \frac{1}{e_{1w}}\;\overline{\overline{\delta_{i+1/2}[z_t]}}^{\,i,\,k+1/2}
106    &\approx \frac{1}{e_{1w}}\; \delta_{i+1/2}[z_{uw}\\
107    r_{2w} &= \frac{1}{e_{2w}}\;\overline{\overline{\delta_{j+1/2}[z_t]}}^{\,j,\,k+1/2}
108    &\approx \frac{1}{e_{2w}}\; \delta_{j+1/2}[z_{vw}]
109  \end{aligned}
110\end{equation}
111
112%gm%  caution I'm not sure the simplification was a good idea!
113
114These slopes are computed once in \rou{ldf\_slp\_init} when \np{ln\_sco}\forcode{=.true.},
115and either \np{ln\_traldf\_hor}\forcode{=.true.} or \np{ln\_dynldf\_hor}\forcode{=.true.}.
116
117\subsection{Slopes for tracer iso-neutral mixing}
118\label{subsec:LDF_slp_iso}
119
120In iso-neutral mixing  $r_1$ and $r_2$ are the slopes between the iso-neutral and computational surfaces.
121Their formulation does not depend on the vertical coordinate used.
122Their discrete formulation is found using the fact that the diffusive fluxes of
123locally referenced potential density (\ie\ $in situ$ density) vanish.
124So, substituting $T$ by $\rho$ in \autoref{eq:TRA_ldf_iso} and setting the diffusive fluxes in
125the three directions to zero leads to the following definition for the neutral slopes:
126
127\begin{equation}
128  \label{eq:LDF_slp_iso}
129  \begin{split}
130    r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac{\delta_{i+1/2}[\rho]}
131    {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,i+1/2,\,k}} \\
132    r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac{\delta_{j+1/2}\left[\rho \right]}
133    {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,j+1/2,\,k}} \\
134    r_{1w} &= \frac{e_{3w}}{e_{1w}}\;
135    \frac{\overline{\overline{\delta_{i+1/2}[\rho]}}^{\,i,\,k+1/2}}
136    {\delta_{k+1/2}[\rho]} \\
137    r_{2w} &= \frac{e_{3w}}{e_{2w}}\;
138    \frac{\overline{\overline{\delta_{j+1/2}[\rho]}}^{\,j,\,k+1/2}}
139    {\delta_{k+1/2}[\rho]}
140  \end{split}
141\end{equation}
142
143%gm% rewrite this as the explanation is not very clear !!!
144%In practice, \autoref{eq:LDF_slp_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 \autoref{eq:LDF_slp_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.
145
146%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{eq:LDF_slp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
147
148%In the $z$-coordinate, the derivative of the  \autoref{eq:LDF_slp_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.
149
150As the mixing is performed along neutral surfaces, the gradient of $\rho$ in \autoref{eq:LDF_slp_iso} has to
151be evaluated at the same local pressure
152(which, in decibars, is approximated by the depth in meters in the model).
153Therefore \autoref{eq:LDF_slp_iso} cannot be used as such,
154but further transformation is needed depending on the vertical coordinate used:
155
156\begin{description}
157
158\item[$z$-coordinate with full step: ]
159  in \autoref{eq:LDF_slp_iso} the densities appearing in the $i$ and $j$ derivatives  are taken at the same depth,
160  thus the $in situ$ density can be used.
161  This is not the case for the vertical derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$,
162  where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{mcdougall_JPO87}
163  (see \autoref{subsec:TRA_bn2}).
164
165\item[$z$-coordinate with partial step: ]
166  this case is identical to the full step case except that at partial step level,
167  the \emph{horizontal} density gradient is evaluated as described in \autoref{sec:TRA_zpshde}.
168
169\item[$s$- or hybrid $s$-$z$- coordinate: ]
170  in the current release of \NEMO, iso-neutral mixing is only employed for $s$-coordinates if
171  the Griffies scheme is used (\np{ln\_traldf\_triad}\forcode{=.true.};
172  see \autoref{apdx:TRIADS}).
173  In other words, iso-neutral mixing will only be accurately represented with a linear equation of state
174  (\np{ln\_seos}\forcode{=.true.}).
175  In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:LDF_slp_iso}
176  will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes.
177
178%gm%
179  Note: The solution for $s$-coordinate passes trough the use of different (and better) expression for
180  the constraint on iso-neutral fluxes.
181  Following \citet{griffies_bk04}, instead of specifying directly that there is a zero neutral diffusive flux of
182  locally referenced potential density, we stay in the $T$-$S$ plane and consider the balance between
183  the neutral direction diffusive fluxes of potential temperature and salinity:
184  \[
185    \alpha \ \textbf{F}(T) = \beta \ \textbf{F}(S)
186  \]
187  % gm{  where vector F is ....}
188
189This constraint leads to the following definition for the slopes:
190
191\[
192  % \label{eq:LDF_slp_iso2}
193  \begin{split}
194    r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac
195    {\alpha_u \;\delta_{i+1/2}[T] - \beta_u \;\delta_{i+1/2}[S]}
196    {\alpha_u \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,i+1/2,\,k}
197      -\beta_u  \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,i+1/2,\,k} } \\
198    r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac
199    {\alpha_v \;\delta_{j+1/2}[T] - \beta_v \;\delta_{j+1/2}[S]}
200    {\alpha_v \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,j+1/2,\,k}
201      -\beta_v  \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,j+1/2,\,k} }    \\
202    r_{1w} &= \frac{e_{3w}}{e_{1w}}\; \frac
203    {\alpha_w \;\overline{\overline{\delta_{i+1/2}[T]}}^{\,i,\,k+1/2}
204      -\beta_w  \;\overline{\overline{\delta_{i+1/2}[S]}}^{\,i,\,k+1/2} }
205    {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]} \\
206    r_{2w} &= \frac{e_{3w}}{e_{2w}}\; \frac
207    {\alpha_w \;\overline{\overline{\delta_{j+1/2}[T]}}^{\,j,\,k+1/2}
208      -\beta_w  \;\overline{\overline{\delta_{j+1/2}[S]}}^{\,j,\,k+1/2} }
209    {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]} \\
210  \end{split}
211\]
212where $\alpha$ and $\beta$, the thermal expansion and saline contraction coefficients introduced in
213\autoref{subsec:TRA_bn2}, have to be evaluated at the three velocity points.
214In order to save computation time, they should be approximated by the mean of their values at $T$-points
215(for example in the case of $\alpha$:
216$\alpha_u=\overline{\alpha_T}^{i+1/2}$$\alpha_v=\overline{\alpha_T}^{j+1/2}$ and
217$\alpha_w=\overline{\alpha_T}^{k+1/2}$).
218
219Note that such a formulation could be also used in the $z$-coordinate and $z$-coordinate with partial steps cases.
220
221\end{description}
222
223This implementation is a rather old one.
224It is similar to the one proposed by \citet{cox_OM87}, except for the background horizontal diffusion.
225Indeed, the \citet{cox_OM87} implementation of isopycnal diffusion in GFDL-type models requires
226a minimum background horizontal diffusion for numerical stability reasons.
227To overcome this problem, several techniques have been proposed in which the numerical schemes of
228the ocean model are modified \citep{weaver.eby_JPO97, griffies.gnanadesikan.ea_JPO98}.
229Griffies's scheme is now available in \NEMO\ if \np{ln\_traldf\_triad}\forcode{ = .true.}; see \autoref{apdx:TRIADS}.
230Here, another strategy is presented \citep{lazar_phd97}:
231a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of
232grid point noise generated by the iso-neutral diffusion operator (\autoref{fig:LDF_ZDF1}).
233This allows an iso-neutral diffusion scheme without additional background horizontal mixing.
234This technique can be viewed as a diffusion operator that acts along large-scale
235(2~$\Delta$x) \gmcomment{2deltax doesnt seem very large scale} iso-neutral surfaces.
236The diapycnal diffusion required for numerical stability is thus minimized and its net effect on the flow is quite small when compared to the effect of an horizontal background mixing.
237
238Nevertheless, this iso-neutral operator does not ensure that variance cannot increase,
239contrary to the \citet{griffies.gnanadesikan.ea_JPO98} operator which has that property.
240
241%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
242\begin{figure}[!ht]
243  \centering
244  \includegraphics[width=0.66\textwidth]{Fig_LDF_ZDF1}
245  \caption{Averaging procedure for isopycnal slope computation}
246  \label{fig:LDF_ZDF1}
247\end{figure}
248%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
249
250%There are three additional questions about the slope calculation.
251%First the expression for the rotation tensor has been obtain assuming the "small slope" approximation, so a bound has to be imposed on slopes.
252%Second, numerical stability issues also require a bound on slopes.
253%Third, the question of boundary condition specified on slopes...
254
255%from griffies: chapter 13.1....
256
257
258
259% In addition and also for numerical stability reasons \citep{cox_OM87, griffies_bk04},
260% the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly
261% to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the
262% surface motivates this flattening of isopycnals near the surface).
263
264For numerical stability reasons \citep{cox_OM87, griffies_bk04}, the slopes must also be bounded by
265the namelist scalar \np{rn\_slpmax} (usually $1/100$) everywhere.
266This constraint is applied in a piecewise linear fashion, increasing from zero at the surface to
267$1/100$ at $70$ metres and thereafter decreasing to zero at the bottom of the ocean
268(the fact that the eddies "feel" the surface motivates this flattening of isopycnals near the surface).
269\colorbox{yellow}{The way slopes are tapered has be checked. Not sure that this is still what is actually done.}
270
271%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
272\begin{figure}[!ht]
273  \centering
274  \includegraphics[width=0.66\textwidth]{Fig_eiv_slp}
275  \caption[Vertical profile of the slope used for lateral mixing in the mixed layer]{
276    Vertical profile of the slope used for lateral mixing in the mixed layer:
277    \textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior,
278    which has to be adjusted at the surface boundary
279    \ie\ it must tend to zero at the surface since there is no mixing across the air-sea interface:
280    wall boundary condition).
281    Nevertheless,
282    the profile between the surface zero value and the interior iso-neutral one is unknown,
283    and especially the value at the base of the mixed layer;
284    \textit{(b)} profile of slope using a linear tapering of the slope near the surface and
285    imposing a maximum slope of 1/100;
286    \textit{(c)} profile of slope actually used in \NEMO:
287    a linear decrease of the slope from zero at the surface to
288    its ocean interior value computed just below the mixed layer.
289    Note the huge change in the slope at the base of the mixed layer between
290    \textit{(b)} and \textit{(c)}.}
291  \label{fig:LDF_eiv_slp}
292\end{figure}
293%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
294
295\colorbox{yellow}{add here a discussion about the flattening of the slopes, vs tapering the coefficient.}
296
297\subsection{Slopes for momentum iso-neutral mixing}
298
299The iso-neutral diffusion operator on momentum is the same as the one used on tracers but
300applied to each component of the velocity separately
301(see \autoref{eq:DYN_ldf_iso} in section~\autoref{subsec:DYN_ldf_iso}).
302The slopes between the surface along which the diffusion operator acts and the surface of computation
303($z$- or $s$-surfaces) are defined at $T$-, $f$-, and \textit{uw}- points for the $u$-component, and $T$-, $f$- and
304\textit{vw}- points for the $v$-component.
305They are computed from the slopes used for tracer diffusion,
306\ie\ \autoref{eq:LDF_slp_geo} and \autoref{eq:LDF_slp_iso}:
307
308\[
309  % \label{eq:LDF_slp_dyn}
310  \begin{aligned}
311    &r_{1t}\ \ = \overline{r_{1u}}^{\,i}       &&&    r_{1f}\ \ &= \overline{r_{1u}}^{\,i+1/2} \\
312    &r_{2f} \ \ = \overline{r_{2v}}^{\,j+1/2} &&&  r_{2t}\ &= \overline{r_{2v}}^{\,j} \\
313    &r_{1uw}  = \overline{r_{1w}}^{\,i+1/2} &&\ \ \text{and} \ \ &   r_{1vw}&= \overline{r_{1w}}^{\,j+1/2} \\
314    &r_{2uw}= \overline{r_{2w}}^{\,j+1/2} &&&         r_{2vw}&= \overline{r_{2w}}^{\,j+1/2}\\
315  \end{aligned}
316\]
317
318The major issue remaining is in the specification of the boundary conditions.
319The same boundary conditions are chosen as those used for lateral diffusion along model level surfaces,
320\ie\ using the shear computed along the model levels and with no additional friction at the ocean bottom
321(see \autoref{sec:LBC_coast}).
322
323
324% ================================================================
325% Lateral Mixing Coefficients
326% ================================================================
327\section[Lateral mixing coefficient (\forcode{nn_aht_ijk_t}, \forcode{nn_ahm_ijk_t})]
328{Lateral mixing coefficient (\protect\np{nn\_aht\_ijk\_t}, \protect\np{nn\_ahm\_ijk\_t})}
329\label{sec:LDF_coef}
330
331The specification of the space variation of the coefficient is made in modules \mdl{ldftra} and \mdl{ldfdyn}.
332The way the mixing coefficients are set in the reference version can be described as follows:
333
334\subsection[Mixing coefficients read from file (\forcode{nn_aht_ijk_t=-20, -30}, \forcode{nn_ahm_ijk_t=-20,-30})]
335{ Mixing coefficients read from file (\protect\np{nn\_aht\_ijk\_t}\forcode{=-20, -30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{=-20, -30})}
336
337Mixing coefficients can be read from file if a particular geographical variation is needed. For example, in the ORCA2 global ocean model,
338the laplacian viscosity operator uses $A^l$~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and
339decreases linearly to $A^l$~= 2.10$^3$ m$^2$/s at the equator \citep{madec.delecluse.ea_JPO96, delecluse.madec_icol99}.
340Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of ORCA2 and ORCA05.
341The provided fields can either be 2d (\np{nn\_aht\_ijk\_t}\forcode{=-20}, \np{nn\_ahm\_ijk\_t}\forcode{=-20}) or 3d (\np{nn\_aht\_ijk\_t}\forcode{=-30}\np{nn\_ahm\_ijk\_t}\forcode{=-30}). They must be given at U, V points for tracers and T, F points for momentum (see \autoref{tab:LDF_files}).
342
343%-------------------------------------------------TABLE---------------------------------------------------
344\begin{table}[htb]
345  \centering
346  \begin{tabular}{|l|l|l|l|}
347    \hline
348    Namelist parameter                       & Input filename                               & dimensions & variable names                      \\  \hline
349    \np{nn\_ahm\_ijk\_t}\forcode{=-20}     & \forcode{eddy_viscosity_2D.nc }            &     $(i,j)$         & \forcode{ahmt_2d, ahmf_2d}  \\  \hline
350    \np{nn\_aht\_ijk\_t}\forcode{=-20}           & \forcode{eddy_diffusivity_2D.nc }           &     $(i,j)$           & \forcode{ahtu_2d, ahtv_2d}    \\   \hline
351    \np{nn\_ahm\_ijk\_t}\forcode{=-30}        & \forcode{eddy_viscosity_3D.nc }            &     $(i,j,k)$          & \forcode{ahmt_3d, ahmf_3d}  \\  \hline
352    \np{nn\_aht\_ijk\_t}\forcode{=-30}     & \forcode{eddy_diffusivity_3D.nc }           &     $(i,j,k)$         & \forcode{ahtu_3d, ahtv_3d}    \\   \hline
353  \end{tabular}
354  \caption{Description of expected input files if mixing coefficients are read from NetCDF files}
355  \label{tab:LDF_files}
356\end{table}
357%--------------------------------------------------------------------------------------------------------------
358
359\subsection[Constant mixing coefficients (\forcode{nn_aht_ijk_t=0}, \forcode{nn_ahm_ijk_t=0})]
360{ Constant mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{=0}, \protect\np{nn\_ahm\_ijk\_t}\forcode{=0})}
361
362If constant, mixing coefficients are set thanks to a velocity and a length scales ($U_{scl}$, $L_{scl}$) such that:
363
364\begin{equation}
365  \label{eq:LDF_constantah}
366  A_o^l = \left\{
367    \begin{aligned}
368      & \frac{1}{2} U_{scl} L_{scl}            & \text{for laplacian operator } \\
369      & \frac{1}{12} U_{scl} L_{scl}^3                    & \text{for bilaplacian operator }
370    \end{aligned}
371  \right.
372\end{equation}
373
374 $U_{scl}$ and $L_{scl}$ are given by the namelist parameters \np{rn\_Ud}, \np{rn\_Uv}, \np{rn\_Ld} and \np{rn\_Lv}.
375
376\subsection[Vertically varying mixing coefficients (\forcode{nn_aht_ijk_t=10}, \forcode{nn_ahm_ijk_t=10})]
377{Vertically varying mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{=10}, \protect\np{nn\_ahm\_ijk\_t}\forcode{=10})}
378
379In the vertically varying case, a hyperbolic variation of the lateral mixing coefficient is introduced in which
380the surface value is given by \autoref{eq:LDF_constantah}, the bottom value is 1/4 of the surface value,
381and the transition takes place around z=500~m with a width of 200~m.
382This profile is hard coded in module \mdl{ldfc1d\_c2d}, but can be easily modified by users.
383
384\subsection[Mesh size dependent mixing coefficients (\forcode{nn_aht_ijk_t=20}, \forcode{nn_ahm_ijk_t=20})]
385{Mesh size dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{=20}, \protect\np{nn\_ahm\_ijk\_t}\forcode{=20})}
386
387In that case, the horizontal variation of the eddy coefficient depends on the local mesh size and
388the type of operator used:
389\begin{equation}
390  \label{eq:LDF_title}
391  A_l = \left\{
392    \begin{aligned}
393      & \frac{1}{2} U_{scl}  \max(e_1,e_2)         & \text{for laplacian operator } \\
394      & \frac{1}{12} U_{scl}  \max(e_1,e_2)^{3}             & \text{for bilaplacian operator }
395    \end{aligned}
396  \right.
397\end{equation}
398where $U_{scl}$ is a user defined velocity scale (\np{rn\_Ud}, \np{rn\_Uv}).
399This variation is intended to reflect the lesser need for subgrid scale eddy mixing where
400the grid size is smaller in the domain.
401It was introduced in the context of the DYNAMO modelling project \citep{willebrand.barnier.ea_PO01}.
402Note that such a grid scale dependance of mixing coefficients significantly increases the range of stability of
403model configurations presenting large changes in grid spacing such as global ocean models.
404Indeed, in such a case, a constant mixing coefficient can lead to a blow up of the model due to
405large coefficient compare to the smallest grid size (see \autoref{sec:TD_forward_imp}),
406especially when using a bilaplacian operator.
407
408\colorbox{yellow}{CASE \np{nn\_aht\_ijk\_t} = 21 to be added}
409
410\subsection[Mesh size and depth dependent mixing coefficients (\forcode{nn_aht_ijk_t=30}, \forcode{nn_ahm_ijk_t=30})]
411{Mesh size and depth dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{=30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{=30})}
412
413The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases above,
414\ie\ a hyperbolic tangent variation with depth associated with a grid size dependence of
415the magnitude of the coefficient.
416
417\subsection[Velocity dependent mixing coefficients (\forcode{nn_aht_ijk_t=31}, \forcode{nn_ahm_ijk_t=31})]
418{Flow dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{=31}, \protect\np{nn\_ahm\_ijk\_t}\forcode{=31})}
419In that case, the eddy coefficient is proportional to the local velocity magnitude so that the Reynolds number $Re =  \lvert U \rvert  e / A_l$ is constant (and here hardcoded to $12$):
420\colorbox{yellow}{JC comment: The Reynolds is effectively set to 12 in the code for both operators but shouldn't it be 2 for Laplacian ?}
421
422\begin{equation}
423  \label{eq:LDF_flowah}
424  A_l = \left\{
425    \begin{aligned}
426      & \frac{1}{12} \lvert U \rvert e          & \text{for laplacian operator } \\
427      & \frac{1}{12} \lvert U \rvert e^3             & \text{for bilaplacian operator }
428    \end{aligned}
429  \right.
430\end{equation}
431
432\subsection[Deformation rate dependent viscosities (\forcode{nn_ahm_ijk_t=32})]
433{Deformation rate dependent viscosities (\protect\np{nn\_ahm\_ijk\_t}\forcode{=32})}
434
435This option refers to the \citep{smagorinsky_MW63} scheme which is here implemented for momentum only. Smagorinsky chose as a
436characteristic time scale $T_{smag}$ the deformation rate and for the lengthscale $L_{smag}$ the maximum wavenumber possible on the horizontal grid, e.g.:
437
438\begin{equation}
439  \label{eq:LDF_smag1}
440  \begin{split}
441    T_{smag}^{-1} & = \sqrt{\left( \partial_x u - \partial_y v\right)^2 + \left( \partial_y u + \partial_x v\right)^} \\
442    L_{smag} & = \frac{1}{\pi}\frac{e_1 e_2}{e_1 + e_2}
443  \end{split}
444\end{equation}
445
446Introducing a user defined constant $C$ (given in the namelist as \np{rn\_csmc}), one can deduce the mixing coefficients as follows:
447
448\begin{equation}
449  \label{eq:LDF_smag2}
450  A_{smag} = \left\{
451    \begin{aligned}
452      & C^2  T_{smag}^{-1}  L_{smag}^2       & \text{for laplacian operator } \\
453      & \frac{C^2}{8} T_{smag}^{-1} L_{smag}^4        & \text{for bilaplacian operator }
454    \end{aligned}
455  \right.
456\end{equation}
457
458For stability reasons, upper and lower limits are applied on the resulting coefficient (see \autoref{sec:TD_forward_imp}) so that:
459\begin{equation}
460  \label{eq:LDF_smag3}
461    \begin{aligned}
462      & C_{min} \frac{1}{2}   \lvert U \rvert  e    < A_{smag} < C_{max} \frac{e^2}{   8\rdt}                 & \text{for laplacian operator } \\
463      & C_{min} \frac{1}{12} \lvert U \rvert  e^3 < A_{smag} < C_{max} \frac{e^4}{64 \rdt}                 & \text{for bilaplacian operator }
464    \end{aligned}
465\end{equation}
466
467where $C_{min}$ and $C_{max}$ are adimensional namelist parameters given by \np{rn\_minfac} and \np{rn\_maxfac} respectively.
468
469\subsection{About space and time varying mixing coefficients}
470
471The following points are relevant when the eddy coefficient varies spatially:
472
473(1) the momentum diffusion operator acting along model level surfaces is written in terms of curl and
474divergent components of the horizontal current (see \autoref{subsec:MB_ldf}).
475Although the eddy coefficient could be set to different values in these two terms,
476this option is not currently available.
477
478(2) with an horizontally varying viscosity, the quadratic integral constraints on enstrophy and on the square of
479the horizontal divergence for operators acting along model-surfaces are no longer satisfied
480(\autoref{sec:INVARIANTS_dynldf_properties}).
481
482% ================================================================
483% Eddy Induced Mixing
484% ================================================================
485\section[Eddy induced velocity (\forcode{ln_ldfeiv=.true.})]
486{Eddy induced velocity (\protect\np{ln\_ldfeiv}\forcode{=.true.})}
487
488\label{sec:LDF_eiv}
489
490%--------------------------------------------namtra_eiv---------------------------------------------------
491
492\begin{listing}
493  \nlst{namtra_eiv}
494  \caption{\texttt{namtra\_eiv}}
495  \label{lst:namtra_eiv}
496\end{listing}
497
498%--------------------------------------------------------------------------------------------------------------
499
500
501%%gm  from Triad appendix  : to be incorporated....
502\gmcomment{
503  Values of iso-neutral diffusivity and GM coefficient are set as described in \autoref{sec:LDF_coef}.
504  If none of the keys \key{traldf\_cNd}, N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and
505  GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and \np{rn\_aeiv\_0}.
506  If 2D-varying coefficients are set with \key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal
507  scale factor according to \autoref{eq:title}
508  \footnote{
509    Except in global ORCA  $0.5^{\circ}$ runs with \key{traldf\_eiv},
510    where $A_l$ is set like $A_e$ but with a minimum vale of $100\;\mathrm{m}^2\;\mathrm{s}^{-1}$
511  }.
512  In idealised setups with \key{traldf\_c2d}, $A_e$ is reduced similarly, but if \key{traldf\_eiv} is set in
513  the global configurations with \key{traldf\_c2d}, a horizontally varying $A_e$ is instead set from
514  the Held-Larichev parameterisation
515  \footnote{
516    In this case, $A_e$ at low latitudes $|\theta|<20^{\circ}$ is further reduced by a factor $|f/f_{20}|$,
517    where $f_{20}$ is the value of $f$ at $20^{\circ}$~N
518  } (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored unless it is zero.
519}
520
521When  \citet{gent.mcwilliams_JPO90} diffusion is used (\np{ln\_ldfeiv}\forcode{=.true.}),
522an eddy induced tracer advection term is added,
523the formulation of which depends on the slopes of iso-neutral surfaces.
524Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces,
525\ie\ \autoref{eq:LDF_slp_geo} is used in $z$-coordinates,
526and the sum \autoref{eq:LDF_slp_geo} + \autoref{eq:LDF_slp_iso} in $s$-coordinates.
527
528If isopycnal mixing is used in the standard way, \ie\ \np{ln\_traldf\_triad}\forcode{=.false.}, the eddy induced velocity is given by:
529\begin{equation}
530  \label{eq:LDF_eiv}
531  \begin{split}
532    u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\
533    v^* & = \frac{1}{e_{1u}e_{3v}}\; \delta_k \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right]\\
534    w^* & = \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\} \\
535  \end{split}
536\end{equation}
537where $A^{eiv}$ is the eddy induced velocity coefficient whose value is set through \np{nn\_aei\_ijk\_t} \nam{tra\_eiv} namelist parameter.
538The three components of the eddy induced velocity are computed in \rou{ldf\_eiv\_trp} and
539added to the eulerian velocity in \rou{tra\_adv} where tracer advection is performed.
540This has been preferred to a separate computation of the advective trends associated with the eiv velocity,
541since it allows us to take advantage of all the advection schemes offered for the tracers
542(see \autoref{sec:TRA_adv}) and not just the $2^{nd}$ order advection scheme as in
543previous releases of OPA \citep{madec.delecluse.ea_NPM98}.
544This is particularly useful for passive tracers where \emph{positivity} of the advection scheme is of
545paramount importance.
546
547At the surface, lateral and bottom boundaries, the eddy induced velocity,
548and thus the advective eddy fluxes of heat and salt, are set to zero.
549The value of the eddy induced mixing coefficient and its space variation is controlled in a similar way as for lateral mixing coefficient described in the preceding subsection (\np{nn\_aei\_ijk\_t}, \np{rn\_Ue}, \np{rn\_Le} namelist parameters).
550\colorbox{yellow}{CASE \np{nn\_aei\_ijk\_t} = 21 to be added}
551
552In case of setting \np{ln\_traldf\_triad}\forcode{ = .true.}, a skew form of the eddy induced advective fluxes is used, which is described in \autoref{apdx:TRIADS}.
553
554% ================================================================
555% Mixed layer eddies
556% ================================================================
557\section[Mixed layer eddies (\forcode{ln_mle=.true.})]
558{Mixed layer eddies (\protect\np{ln\_mle}\forcode{=.true.})}
559
560\label{sec:LDF_mle}
561
562%--------------------------------------------namtra_eiv---------------------------------------------------
563
564\begin{listing}
565  \nlst{namtra_mle}
566  \caption{\texttt{namtra\_mle}}
567  \label{lst:namtra_mle}
568\end{listing}
569
570%--------------------------------------------------------------------------------------------------------------
571
572If  \np{ln\_mle}\forcode{=.true.} in \nam{tra\_mle} namelist, a parameterization of the mixing due to unresolved mixed layer instabilities is activated (\citet{foxkemper.ferrari_JPO08}). Additional transport is computed in \rou{ldf\_mle\_trp} and added to the eulerian transport in \rou{tra\_adv} as done for eddy induced advection.
573
574\colorbox{yellow}{TBC}
575
576\biblio
577
578\pindex
579
580\end{document}
Note: See TracBrowser for help on using the repository browser.