Changeset 11422 for NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LDF.tex

Timestamp:

2019-08-08T15:40:47+02:00 (5 years ago)

Author:

jchanut

Message:

#1791, merge with trunk

Location:

NEMO/branches/2019/fix_vvl_ticket1791/doc

Files:

• . (modified) (1 prop)
• latex (modified) (1 prop)
• latex/NEMO (modified) (1 prop)
• latex/NEMO/subfiles/chap_LDF.tex (modified) (18 diffs)

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex

@@ -1 @@
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
-*.ilg
-*.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*
+figures

@@ -1 @@
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
-*.ilg
-*.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*
+figures

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LDF.tex

@@ -23 +23 @@
 These three aspects of the lateral diffusion are set through namelist parameters
 (see the \ngn{nam\_traldf} and \ngn{nam\_dynldf} below).
-Note that this chapter describes the standard implementation of iso-neutral tracer mixing,
-and Griffies's implementation, which is used if \np{traldf\_grif}\forcode{ = .true.},
-is described in Appdx\autoref{apdx:triad}
-
-%-----------------------------------nam_traldf - nam_dynldf--------------------------------------------
+Note that this chapter describes the standard implementation of iso-neutral tracer mixing. 
+Griffies's implementation, which is used if \np{ln\_traldf\_triad}\forcode{ = .true.},
+is described in \autoref{apdx:triad}
+
+%-----------------------------------namtra_ldf - namdyn_ldf--------------------------------------------
 
 \nlst{namtra_ldf} 
@@ -34 +34 @@
 %--------------------------------------------------------------------------------------------------------------
 
+% ================================================================
+% Lateral Mixing Operator
+% ================================================================
+\section[Lateral mixing operators]
+{Lateral mixing operators}
+\label{sec:LDF_op}
+We 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}.
+
+\subsection[No lateral mixing (\forcode{ln_traldf_OFF}, \forcode{ln_dynldf_OFF})]
+{No lateral mixing (\protect\np{ln\_traldf\_OFF}, \protect\np{ln\_dynldf\_OFF})}
+
+It is possible to run without explicit lateral diffusion on tracers (\protect\np{ln\_traldf\_OFF}\forcode{ = .true.}) and/or 
+momentum (\protect\np{ln\_dynldf\_OFF}\forcode{ = .true.}). The latter option is even recommended if using the 
+UBS advection scheme on momentum (\np{ln\_dynadv\_ubs}\forcode{ = .true.},
+see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes.
+
+\subsection[Laplacian mixing (\forcode{ln_traldf_lap}, \forcode{ln_dynldf_lap})]
+{Laplacian mixing (\protect\np{ln\_traldf\_lap}, \protect\np{ln\_dynldf\_lap})}
+Setting \protect\np{ln\_traldf\_lap}\forcode{ = .true.} and/or \protect\np{ln\_dynldf\_lap}\forcode{ = .true.} enables 
+a second order diffusion on tracers and momentum respectively. Note that in \NEMO 4, one can not combine 
+Laplacian and Bilaplacian operators for the same variable.
+
+\subsection[Bilaplacian mixing (\forcode{ln_traldf_blp}, \forcode{ln_dynldf_blp})]
+{Bilaplacian mixing (\protect\np{ln\_traldf\_blp}, \protect\np{ln\_dynldf\_blp})}
+Setting \protect\np{ln\_traldf\_blp}\forcode{ = .true.} and/or \protect\np{ln\_dynldf\_blp}\forcode{ = .true.} enables 
+a fourth order diffusion on tracers and momentum respectively. It is implemented by calling the above Laplacian operator twice. 
+We stress again that from \NEMO 4, the simultaneous use Laplacian and Bilaplacian operators is not allowed.
 
 % ================================================================
 % Direction of lateral Mixing
 % ================================================================
-\section{Direction of lateral mixing (\protect\mdl{ldfslp})}
+\section[Direction of lateral mixing (\textit{ldfslp.F90})]
+{Direction of lateral mixing (\protect\mdl{ldfslp})}
 \label{sec:LDF_slp}
 
@@ -44 +72 @@
 \gmcomment{
   we should emphasize here that the implementation is a rather old one.
-  Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme.
+  Better work can be achieved by using \citet{griffies.gnanadesikan.ea_JPO98, griffies_bk04} iso-neutral scheme.
 }
 
@@ -87 +115 @@
 %gm%  caution I'm not sure the simplification was a good idea! 
 
-These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}\forcode{ = .true.}rue,
+These slopes are computed once in \rou{ldf\_slp\_init} when \np{ln\_sco}\forcode{ = .true.},
 and either \np{ln\_traldf\_hor}\forcode{ = .true.} or \np{ln\_dynldf\_hor}\forcode{ = .true.}. 
 
@@ -119 +147 @@
 %In practice, \autoref{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 \autoref{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. 
 
-%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{McDougall1987}, therefore in \autoref{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).
+%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{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).
 
 %In the $z$-coordinate, the derivative of the  \autoref{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. 
@@ -135 +163 @@
   thus the $in situ$ density can be used.
   This is not the case for the vertical derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$,
-  where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{McDougall1987}
+  where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{mcdougall_JPO87}
   (see \autoref{subsec:TRA_bn2}). 
 
@@ -144 +172 @@
 \item[$s$- or hybrid $s$-$z$- coordinate: ]
   in the current release of \NEMO, iso-neutral mixing is only employed for $s$-coordinates if
-  the Griffies scheme is used (\np{traldf\_grif}\forcode{ = .true.};
-  see Appdx \autoref{apdx:triad}).
+  the Griffies scheme is used (\np{ln\_traldf\_triad}\forcode{ = .true.};
+  see \autoref{apdx:triad}).
   In other words, iso-neutral mixing will only be accurately represented with a linear equation of state
-  (\np{nn\_eos}\forcode{ = 1..2}).
+  (\np{ln\_seos}\forcode{ = .true.}).
   In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:ldfslp_iso}
   will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes.
@@ -154 +182 @@
   Note: The solution for $s$-coordinate passes trough the use of different (and better) expression for
   the constraint on iso-neutral fluxes.
-  Following \citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral diffusive flux of
+  Following \citet{griffies_bk04}, instead of specifying directly that there is a zero neutral diffusive flux of
   locally referenced potential density, we stay in the $T$-$S$ plane and consider the balance between
   the neutral direction diffusive fluxes of potential temperature and salinity:
@@ -197 +225 @@
 
 This implementation is a rather old one.
-It is similar to the one proposed by Cox [1987], except for the background horizontal diffusion.
-Indeed, the Cox implementation of isopycnal diffusion in GFDL-type models requires
+It is similar to the one proposed by \citet{cox_OM87}, except for the background horizontal diffusion.
+Indeed, the \citet{cox_OM87} implementation of isopycnal diffusion in GFDL-type models requires
 a minimum background horizontal diffusion for numerical stability reasons.
 To overcome this problem, several techniques have been proposed in which the numerical schemes of
-the ocean model are modified \citep{Weaver_Eby_JPO97, Griffies_al_JPO98}.
-Griffies's scheme is now available in \NEMO if \np{traldf\_grif\_iso} is set true; see Appdx \autoref{apdx:triad}.
-Here, another strategy is presented \citep{Lazar_PhD97}:
+the ocean model are modified \citep{weaver.eby_JPO97, griffies.gnanadesikan.ea_JPO98}.
+Griffies's scheme is now available in \NEMO if \np{ln\_traldf\_triad}=\forcode{= .true.}; see \autoref{apdx:triad}.
+Here, another strategy is presented \citep{lazar_phd97}:
 a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of
 grid point noise generated by the iso-neutral diffusion operator (\autoref{fig:LDF_ZDF1}).
@@ -212 +240 @@
 
 Nevertheless, this iso-neutral operator does not ensure that variance cannot increase,
-contrary to the \citet{Griffies_al_JPO98} operator which has that property. 
+contrary to the \citet{griffies.gnanadesikan.ea_JPO98} operator which has that property. 
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht]
   \begin{center}
-    \includegraphics[width=0.70\textwidth]{Fig_LDF_ZDF1}
+    \includegraphics[width=\textwidth]{Fig_LDF_ZDF1}
     \caption {
       \protect\label{fig:LDF_ZDF1}
@@ -235 +263 @@
 
 
-% In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04}, 
+% In addition and also for numerical stability reasons \citep{cox_OM87, griffies_bk04}, 
 % the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly 
 % to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the 
 % surface motivates this flattening of isopycnals near the surface).
 
-For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also be bounded by
-$1/100$ everywhere.
+For numerical stability reasons \citep{cox_OM87, griffies_bk04}, the slopes must also be bounded by
+the namelist scalar \np{rn\_slpmax} (usually $1/100$) everywhere.
 This constraint is applied in a piecewise linear fashion, increasing from zero at the surface to
 $1/100$ at $70$ metres and thereafter decreasing to zero at the bottom of the ocean
 (the fact that the eddies "feel" the surface motivates this flattening of isopycnals near the surface).
+\colorbox{yellow}{The way slopes are tapered has be checked. Not sure that this is still what is actually done.}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht]
   \begin{center}
-    \includegraphics[width=0.70\textwidth]{Fig_eiv_slp}
+    \includegraphics[width=\textwidth]{Fig_eiv_slp}
     \caption{
       \protect\label{fig:eiv_slp}
@@ -297 +326 @@
 (see \autoref{sec:LBC_coast}).
 
-
-% ================================================================
-% Lateral Mixing Operator
-% ================================================================
-\section{Lateral mixing operators (\protect\mdl{traldf}, \protect\mdl{traldf}) }
-\label{sec:LDF_op}
-
-
    
 % ================================================================
 % Lateral Mixing Coefficients
 % ================================================================
-\section{Lateral mixing coefficient (\protect\mdl{ldftra}, \protect\mdl{ldfdyn}) }
+\section[Lateral mixing coefficient (\forcode{nn_aht_ijk_t}, \forcode{nn_ahm_ijk_t})]
+{Lateral mixing coefficient (\protect\np{nn\_aht\_ijk\_t}, \protect\np{nn\_ahm\_ijk\_t})}
 \label{sec:LDF_coef}
 
-Introducing a space variation in the lateral eddy mixing coefficients changes the model core memory requirement,
-adding up to four extra three-dimensional arrays for the geopotential or isopycnal second order operator applied to 
-momentum.
-Six CPP keys control the space variation of eddy coefficients: three for momentum and three for tracer.
-The three choices allow:
-a space variation in the three space directions (\key{traldf\_c3d},  \key{dynldf\_c3d}),
-in the horizontal plane (\key{traldf\_c2d},  \key{dynldf\_c2d}),
-or in the vertical only (\key{traldf\_c1d},  \key{dynldf\_c1d}).
-The default option is a constant value over the whole ocean on both momentum and tracers. 
-   
-The number of additional arrays that have to be defined and the gridpoint position at which
-they are defined depend on both the space variation chosen and the type of operator used.
-The resulting eddy viscosity and diffusivity coefficients can be a function of more than one variable.
-Changes in the computer code when switching from one option to another have been minimized by
-introducing the eddy coefficients as statement functions
-(include file \textit{ldftra\_substitute.h90} and \textit{ldfdyn\_substitute.h90}).
-The functions are replaced by their actual meaning during the preprocessing step (CPP).
-The specification of the space variation of the coefficient is made in \mdl{ldftra} and \mdl{ldfdyn},
-or more precisely in include files \textit{traldf\_cNd.h90} and \textit{dynldf\_cNd.h90}, with N=1, 2 or 3.
-The user can modify these include files as he/she wishes.
-The way the mixing coefficient are set in the reference version can be briefly described as follows:
-
-\subsubsection{Constant mixing coefficients (default option)}
-When none of the \key{dynldf\_...} and \key{traldf\_...} keys are defined,
-a constant value is used over the whole ocean for momentum and tracers,
-which is specified through the \np{rn\_ahm0} and \np{rn\_aht0} namelist parameters.
-
-\subsubsection{Vertically varying mixing coefficients (\protect\key{traldf\_c1d} and \key{dynldf\_c1d})} 
-The 1D option is only available when using the $z$-coordinate with full step.
-Indeed in all the other types of vertical coordinate,
-the depth is a 3D function of (\textbf{i},\textbf{j},\textbf{k}) and therefore,
-introducing depth-dependent mixing coefficients will require 3D arrays.
-In the 1D option, a hyperbolic variation of the lateral mixing coefficient is introduced in which
-the surface value is \np{rn\_aht0} (\np{rn\_ahm0}), the bottom value is 1/4 of the surface value,
-and the transition takes place around z=300~m with a width of 300~m
-(\ie both the depth and the width of the inflection point are set to 300~m).
-This profile is hard coded in file \textit{traldf\_c1d.h90}, but can be easily modified by users.
-
-\subsubsection{Horizontally varying mixing coefficients (\protect\key{traldf\_c2d} and \protect\key{dynldf\_c2d})}
-By default the horizontal variation of the eddy coefficient depends on the local mesh size and
+The specification of the space variation of the coefficient is made in modules \mdl{ldftra} and \mdl{ldfdyn}. 
+The way the mixing coefficients are set in the reference version can be described as follows:
+
+\subsection[Mixing coefficients read from file (\forcode{nn_aht_ijk_t = -20, -30}, \forcode{nn_ahm_ijk_t = -20,-30})]
+{ Mixing coefficients read from file (\protect\np{nn\_aht\_ijk\_t}\forcode{ = -20, -30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = -20, -30})}
+
+Mixing coefficients can be read from file if a particular geographical variation is needed. For example, in the ORCA2 global ocean model, 
+the laplacian viscosity operator uses $A^l$~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and
+decreases linearly to $A^l$~= 2.10$^3$ m$^2$/s at the equator \citep{madec.delecluse.ea_JPO96, delecluse.madec_icol99}. 
+Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of ORCA2 and ORCA05. 
+The 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}).
+
+%-------------------------------------------------TABLE---------------------------------------------------
+\begin{table}[htb]
+  \begin{center}
+    \begin{tabular}{|l|l|l|l|}
+      \hline
+      Namelist parameter                        & Input filename                               & dimensions & variable names                      \\  \hline
+      \np{nn\_ahm\_ijk\_t}\forcode{ = -20}       & \forcode{eddy_viscosity_2D.nc }            &     $(i,j)$         & \forcode{ahmt_2d, ahmf_2d}  \\  \hline
+      \np{nn\_aht\_ijk\_t}\forcode{ = -20}           & \forcode{eddy_diffusivity_2D.nc }           &     $(i,j)$          & \forcode{ahtu_2d, ahtv_2d}    \\   \hline
+      \np{nn\_ahm\_ijk\_t}\forcode{ = -30}       & \forcode{eddy_viscosity_3D.nc }            &     $(i,j,k)$          & \forcode{ahmt_3d, ahmf_3d}  \\  \hline
+      \np{nn\_aht\_ijk\_t}\forcode{ = -30}       & \forcode{eddy_diffusivity_3D.nc }           &     $(i,j,k)$         & \forcode{ahtu_3d, ahtv_3d}    \\   \hline
+    \end{tabular}
+    \caption{
+      \protect\label{tab:LDF_files}
+      Description of expected input files if mixing coefficients are read from NetCDF files.
+    }
+  \end{center}
+\end{table}
+%--------------------------------------------------------------------------------------------------------------
+
+\subsection[Constant mixing coefficients (\forcode{nn_aht_ijk_t = 0}, \forcode{nn_ahm_ijk_t = 0})]
+{ Constant mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 0}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 0})}
+
+If constant, mixing coefficients are set thanks to a velocity and a length scales ($U_{scl}$, $L_{scl}$) such that:
+
+\begin{equation}
+  \label{eq:constantah}
+  A_o^l = \left\{
+    \begin{aligned}
+      & \frac{1}{2} U_{scl} L_{scl}            & \text{for laplacian operator } \\
+      & \frac{1}{12} U_{scl} L_{scl}^3                    & \text{for bilaplacian operator }
+    \end{aligned}
+  \right.
+\end{equation}
+
+ $U_{scl}$ and $L_{scl}$ are given by the namelist parameters \np{rn\_Ud}, \np{rn\_Uv}, \np{rn\_Ld} and \np{rn\_Lv}.
+
+\subsection[Vertically varying mixing coefficients (\forcode{nn_aht_ijk_t = 10}, \forcode{nn_ahm_ijk_t = 10})]
+{Vertically varying mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 10}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 10})}
+
+In the vertically varying case, a hyperbolic variation of the lateral mixing coefficient is introduced in which
+the surface value is given by \autoref{eq:constantah}, the bottom value is 1/4 of the surface value,
+and the transition takes place around z=500~m with a width of 200~m.
+This profile is hard coded in module \mdl{ldfc1d\_c2d}, but can be easily modified by users.
+
+\subsection[Mesh size dependent mixing coefficients (\forcode{nn_aht_ijk_t = 20}, \forcode{nn_ahm_ijk_t = 20})]
+{Mesh size dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 20}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 20})}
+
+In that case, the horizontal variation of the eddy coefficient depends on the local mesh size and
 the type of operator used:
 \begin{equation}
@@ -357 +399 @@
   A_l = \left\{
     \begin{aligned}
-      & \frac{\max(e_1,e_2)}{e_{max}} A_o^l           & \text{for laplacian operator } \\
-      & \frac{\max(e_1,e_2)^{3}}{e_{max}^{3}} A_o^l          & \text{for bilaplacian operator }
+      & \frac{1}{2} U_{scl}  \max(e_1,e_2)         & \text{for laplacian operator } \\
+      & \frac{1}{12} U_{scl}  \max(e_1,e_2)^{3}             & \text{for bilaplacian operator }
     \end{aligned}
   \right.
 \end{equation}
-where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked ocean domain,
-and $A_o^l$ is the \np{rn\_ahm0} (momentum) or \np{rn\_aht0} (tracer) namelist parameter.
+where $U_{scl}$ is a user defined velocity scale (\np{rn\_Ud}, \np{rn\_Uv}).
 This variation is intended to reflect the lesser need for subgrid scale eddy mixing where
 the grid size is smaller in the domain.
-It was introduced in the context of the DYNAMO modelling project \citep{Willebrand_al_PO01}.
-Note that such a grid scale dependance of mixing coefficients significantly increase the range of stability of
-model configurations presenting large changes in grid pacing such as global ocean models.
+It was introduced in the context of the DYNAMO modelling project \citep{willebrand.barnier.ea_PO01}.
+Note that such a grid scale dependance of mixing coefficients significantly increases the range of stability of
+model configurations presenting large changes in grid spacing such as global ocean models.
 Indeed, in such a case, a constant mixing coefficient can lead to a blow up of the model due to
 large coefficient compare to the smallest grid size (see \autoref{sec:STP_forward_imp}),
 especially when using a bilaplacian operator.
 
-Other formulations can be introduced by the user for a given configuration.
-For example, in the ORCA2 global ocean model (see Configurations),
-the laplacian viscosity operator uses \np{rn\_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and
-decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s at the equator \citep{Madec_al_JPO96, Delecluse_Madec_Bk00}.
-This modification can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}.
-Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of
-ORCA2 and ORCA05 (see \&namcfg namelist).
-
-\subsubsection{Space varying mixing coefficients (\protect\key{traldf\_c3d} and \key{dynldf\_c3d})}
-
-The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases,
+\colorbox{yellow}{CASE \np{nn\_aht\_ijk\_t} = 21 to be added}
+
+\subsection[Mesh size and depth dependent mixing coefficients (\forcode{nn_aht_ijk_t = 30}, \forcode{nn_ahm_ijk_t = 30})]
+{Mesh size and depth dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 30})}
+
+The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases above,
 \ie a hyperbolic tangent variation with depth associated with a grid size dependence of
 the magnitude of the coefficient. 
 
-\subsubsection{Space and time varying mixing coefficients}
-
-There is no default specification of space and time varying mixing coefficient. 
-The only case available is specific to the ORCA2 and ORCA05 global ocean configurations.
-It provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and
-eddy induced velocity (ORCA05) that depends on the local growth rate of baroclinic instability.
-This specification is actually used when an ORCA key and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined.
+\subsection[Velocity dependent mixing coefficients (\forcode{nn_aht_ijk_t = 31}, \forcode{nn_ahm_ijk_t = 31})]
+{Flow dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 31}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 31})}
+In 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$):
+\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 ?}
+
+\begin{equation}
+  \label{eq:flowah}
+  A_l = \left\{
+    \begin{aligned}
+      & \frac{1}{12} \lvert U \rvert e          & \text{for laplacian operator } \\
+      & \frac{1}{12} \lvert U \rvert e^3             & \text{for bilaplacian operator } 
+    \end{aligned}
+  \right.
+\end{equation}
+
+\subsection[Deformation rate dependent viscosities (\forcode{nn_ahm_ijk_t = 32})]
+{Deformation rate dependent viscosities (\protect\np{nn\_ahm\_ijk\_t}\forcode{ = 32})}
+
+This option refers to the \citep{smagorinsky_MW63} scheme which is here implemented for momentum only. Smagorinsky chose as a 
+characteristic time scale $T_{smag}$ the deformation rate and for the lengthscale $L_{smag}$ the maximum wavenumber possible on the horizontal grid, e.g.:
+
+\begin{equation}
+  \label{eq:smag1}
+  \begin{split}
+    T_{smag}^{-1} & = \sqrt{\left( \partial_x u - \partial_y v\right)^2 + \left( \partial_y u + \partial_x v\right)^2  } \\
+    L_{smag} & = \frac{1}{\pi}\frac{e_1 e_2}{e_1 + e_2}
+  \end{split}
+\end{equation}
+
+Introducing a user defined constant $C$ (given in the namelist as \np{rn\_csmc}), one can deduce the mixing coefficients as follows:
+
+\begin{equation}
+  \label{eq:smag2}
+  A_{smag} = \left\{
+    \begin{aligned}
+      & C^2  T_{smag}^{-1}  L_{smag}^2       & \text{for laplacian operator } \\
+      & \frac{C^2}{8} T_{smag}^{-1} L_{smag}^4        & \text{for bilaplacian operator } 
+    \end{aligned}
+  \right.
+\end{equation}
+
+For stability reasons, upper and lower limits are applied on the resulting coefficient (see \autoref{sec:STP_forward_imp}) so that:
+\begin{equation}
+  \label{eq:smag3}
+    \begin{aligned}
+      & C_{min} \frac{1}{2}   \lvert U \rvert  e    < A_{smag} < C_{max} \frac{e^2}{   8\rdt}                 & \text{for laplacian operator } \\
+      & C_{min} \frac{1}{12} \lvert U \rvert  e^3 < A_{smag} < C_{max} \frac{e^4}{64 \rdt}                 & \text{for bilaplacian operator } 
+    \end{aligned}
+\end{equation}
+
+where $C_{min}$ and $C_{max}$ are adimensional namelist parameters given by \np{rn\_minfac} and \np{rn\_maxfac} respectively.
+
+\subsection{About space and time varying mixing coefficients}
 
 The following points are relevant when the eddy coefficient varies spatially:
@@ -406 +488 @@
 (\autoref{sec:dynldf_properties}).
 
-(3) for isopycnal diffusion on momentum or tracers, an additional purely horizontal background diffusion with
-uniform coefficient can be added by setting a non zero value of \np{rn\_ahmb0} or \np{rn\_ahtb0},
-a background horizontal eddy viscosity or diffusivity coefficient
-(namelist parameters whose default values are $0$).
-However, the technique used to compute the isopycnal slopes is intended to get rid of such a background diffusion,
-since it introduces spurious diapycnal diffusion (see \autoref{sec:LDF_slp}).
-
-(4) when an eddy induced advection term is used (\key{traldf\_eiv}),
-$A^{eiv}$, the eddy induced coefficient has to be defined.
-Its space variations are controlled by the same CPP variable as for the eddy diffusivity coefficient
-(\ie \key{traldf\_cNd}). 
-
-(5) the eddy coefficient associated with a biharmonic operator must be set to a \emph{negative} value.
-
-(6) it is possible to use both the laplacian and biharmonic operators concurrently.
-
-(7) it is possible to run without explicit lateral diffusion on momentum
-(\np{ln\_dynldf\_lap}\forcode{ = .?.}\np{ln\_dynldf\_bilap}\forcode{ = .false.}).
-This is recommended when using the UBS advection scheme on momentum (\np{ln\_dynadv\_ubs}\forcode{ = .true.},
-see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes.
-
 % ================================================================
 % Eddy Induced Mixing
 % ================================================================
-\section{Eddy induced velocity (\protect\mdl{traadv\_eiv}, \protect\mdl{ldfeiv})}
+\section[Eddy induced velocity (\forcode{ln_ldfeiv = .true.})]
+{Eddy induced velocity (\protect\np{ln\_ldfeiv}\forcode{ = .true.})}
+
 \label{sec:LDF_eiv}
+
+%--------------------------------------------namtra_eiv---------------------------------------------------
+
+\nlst{namtra_eiv}
+
+%--------------------------------------------------------------------------------------------------------------
+
 
 %%gm  from Triad appendix  : to be incorporated....
@@ -453 +523 @@
 }
 
-When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined),
+When  \citet{gent.mcwilliams_JPO90} diffusion is used (\np{ln\_ldfeiv}\forcode{ = .true.}),
 an eddy induced tracer advection term is added,
 the formulation of which depends on the slopes of iso-neutral surfaces.
@@ -459 +529 @@
 \ie \autoref{eq:ldfslp_geo} is used in $z$-coordinates,
 and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $s$-coordinates.
-The eddy induced velocity is given by: 
+
+If isopycnal mixing is used in the standard way, \ie \np{ln\_traldf\_triad}\forcode{ = .false.}, the eddy induced velocity is given by: 
 \begin{equation}
   \label{eq:ldfeiv}
@@ -468 +539 @@
   \end{split}
 \end{equation}
-where $A^{eiv}$ is the eddy induced velocity coefficient whose value is set through \np{rn\_aeiv},
-a \textit{nam\_traldf} namelist parameter.
-The three components of the eddy induced velocity are computed and
-add to the eulerian velocity in \mdl{traadv\_eiv}.
+where $A^{eiv}$ is the eddy induced velocity coefficient whose value is set through \np{nn\_aei\_ijk\_t} \ngn{namtra\_eiv} namelist parameter. 
+The three components of the eddy induced velocity are computed in \rou{ldf\_eiv\_trp} and
+added to the eulerian velocity in \rou{tra\_adv} where tracer advection is performed.
 This has been preferred to a separate computation of the advective trends associated with the eiv velocity,
 since it allows us to take advantage of all the advection schemes offered for the tracers
 (see \autoref{sec:TRA_adv}) and not just the $2^{nd}$ order advection scheme as in
-previous releases of OPA \citep{Madec1998}.
+previous releases of OPA \citep{madec.delecluse.ea_NPM98}.
 This is particularly useful for passive tracers where \emph{positivity} of the advection scheme is of
 paramount importance. 
@@ -481 +551 @@
 At the surface, lateral and bottom boundaries, the eddy induced velocity,
 and thus the advective eddy fluxes of heat and salt, are set to zero. 
+The 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). 
+\colorbox{yellow}{CASE \np{nn\_aei\_ijk\_t} = 21 to be added}
+
+In 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:triad}.
+
+% ================================================================
+% Mixed layer eddies
+% ================================================================
+\section[Mixed layer eddies (\forcode{ln_mle = .true.})]
+{Mixed layer eddies (\protect\np{ln\_mle}\forcode{ = .true.})}
+
+\label{sec:LDF_mle}
+
+%--------------------------------------------namtra_eiv---------------------------------------------------
+
+\nlst{namtra_mle}
+
+%--------------------------------------------------------------------------------------------------------------
+
+If  \np{ln\_mle}\forcode{ = .true.} in \ngn{namtra\_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.
+
+\colorbox{yellow}{TBC}
 
 \biblio

@@ -1 @@
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
-*.ilg
-*.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*
+figures