Index: NEMO/trunk/doc/latex/NEMO/main/bibliography.bib
===================================================================
--- NEMO/trunk/doc/latex/NEMO/main/bibliography.bib	(revision 11335)
+++ NEMO/trunk/doc/latex/NEMO/main/bibliography.bib	(revision 11336)
@@ -928,4 +928,18 @@
 }
 
+@article{         foxkemper.ferrari_JPO08,
+  title         = "Parameterization of Mixed Layer Eddies. Part I: Theory and Diagnosis",
+  pages         = "1145--1165",
+  journal       = "Journal of Physical Oceanography",
+  volume        = "38",
+  number        = "6",
+  author        = "B. Fox-Kemper and R. Ferrari and B. Hallberg",
+  year          = "2008",
+  month         = "jun",
+  publisher     = "American Meteorological Society",
+  issn          = "1520-0485",
+  doi           = "10.1175/2007JPO3792.1"
+}
+
 @article{         galperin.kantha.ea_JAS88,
   title         = "A quasi-equilibrium turbulent energy model for geophysical
@@ -2639,4 +2653,18 @@
 }
 
+@article{         smagorinsky_MW63,
+  title         = "General circulation experiments with the primitive equations: I. The basic experiment ",
+  pages         = "99--164",
+  journal       = "Monthly Weather Review",
+  volume        = "91",
+  number        = "3",
+  author        = "J. Smagorinsky",
+  year          = "1963",
+  month         = "mar",
+  publisher     = "American Meteorological Society",
+  issn          = "1520-0493",
+  doi           = "10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2"
+}
+
 @article{         song.haidvogel_JCP94,
   title         = "A semi-implicit ocean circulation model using a
Index: NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex
===================================================================
--- NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex	(revision 11335)
+++ NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex	(revision 11336)
@@ -23,9 +23,9 @@
 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,4 +34,31 @@
 %--------------------------------------------------------------------------------------------------------------
 
+% ================================================================
+% 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.
 
 % ================================================================
@@ -88,5 +115,5 @@
 %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.}. 
 
@@ -145,8 +172,8 @@
 \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.
@@ -198,10 +225,10 @@
 
 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.gnanadesikan.ea_JPO98}.
-Griffies's scheme is now available in \NEMO if \np{traldf\_grif\_iso} is set true; see Appdx \autoref{apdx:triad}.
+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
@@ -242,8 +269,9 @@
 
 For numerical stability reasons \citep{cox_OM87, griffies_bk04}, the slopes must also be bounded by
-$1/100$ everywhere.
+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.}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -298,63 +326,72 @@
 (see \autoref{sec:LBC_coast}).
 
-
-% ================================================================
-% Lateral Mixing Operator
-% ================================================================
-\section[Lateral mixing operators (\textit{traldf.F90})]
-{Lateral mixing operators (\protect\mdl{traldf}, \protect\mdl{traldf})}
-\label{sec:LDF_op}
-
-
    
 % ================================================================
 % Lateral Mixing Coefficients
 % ================================================================
-\section[Lateral mixing coefficient (\textit{ldftra.F90}, \textit{ldfdyn.F90})]
-{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 (\texttt{\textbf{key\_traldf\_c1d}} and \texttt{\textbf{key\_dynldf\_c1d}})]
-{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 (\texttt{\textbf{key\_traldf\_c2d}} and \texttt{\textbf{key\_dynldf\_c2d}})]
-{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}
@@ -362,42 +399,82 @@
   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{\max(e_1,e_2)}{e_{ref}} A_o^l  			& \text{for laplacian operator } \\
+      & \frac{\max(e_1,e_2)^{3}}{e_{ref}^{3}} A_o^l          & \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 $e_{ref}$ is a reference grid size harcoded to a $1^{\circ}$ grid size (\ie $e_{ref}\approx 111 km$),
+and $A_o^l$ is the user defined mixing coefficient defined according to  \autoref{eq:constantah}.
 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.barnier.ea_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.
+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.delecluse.ea_JPO96, delecluse.madec_icol99}.
-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 (\texttt{\textbf{key\_traldf\_c3d}} and \texttt{\textbf{key\_dynldf\_c3d}})]
-{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:
@@ -412,31 +489,18 @@
 (\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 (\textit{traadv\_eiv.F90}, \textit{ldfeiv.F90})]
-{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....
@@ -460,5 +524,5 @@
 }
 
-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.
@@ -466,5 +530,6 @@
 \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}
@@ -475,8 +540,7 @@
   \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
@@ -488,4 +552,26 @@
 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
