Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 2243)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 2282)
@@ -7,6 +7,8 @@
 
 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN.
-\gmcomment{Steven remark (not taken into account : problem here with turbulent vs turbulence.  I've changed "turbulent closure" to "turbulence closure", "turbulent mixing length" to "turbulence mixing length", but I've left "turbulent kinetic energy" alone - though I think it is an historical abberation!
-Gurvan :  I kept "turbulent closure etc "...}
+
+
+\newpage
+$\ $\newline    % force a new ligne
 
 
@@ -38,7 +40,7 @@
 general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 
 These trends can be computed using either a forward time stepping scheme 
-(namelist parameter \np{np\_zdfexp}=true) or a backward time stepping 
-scheme (\np{np\_zdfexp}=false) depending on the magnitude of the mixing 
-coefficients, and thus of the formulation used (see \S\ref{DOM_nxt}).
+(namelist parameter \np{ln\_zdfexp}=true) or a backward time stepping 
+scheme (\np{ln\_zdfexp}=false) depending on the magnitude of the mixing 
+coefficients, and thus of the formulation used (see \S\ref{STP}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -57,9 +59,8 @@
 \begin{align*} 
 A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1} 	\\
-\\
 A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1}
 \end{align*}
 
-These values are set through the \np{avm0} and \np{avt0} namelist parameters. 
+These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters. 
 In all cases, do not use values smaller that those associated with the molecular 
 viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, 
@@ -74,5 +75,5 @@
 
 %--------------------------------------------namric---------------------------------------------------------
-\namdisplay{namric}
+\namdisplay{namzdf_ric}
 %--------------------------------------------------------------------------------------------------------------
 
@@ -84,10 +85,9 @@
 growth of Kelvin-Helmholtz like instabilities leads to a dependency between the 
 vertical eddy coefficients and the local Richardson number ($i.e.$ the 
-ratio of stratification to vertical shear). Following \citet{PacPhil1981}, the following 
+ratio of stratification to vertical shear). Following \citet{Pacanowski_Philander_JPO81}, the following 
 formulation has been implemented:
 \begin{equation} \label{Eq_zdfric}
    \left\{      \begin{aligned}
          A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT}       \\
-         \\
          A^{vm} &= \frac{A^{vT}        }{\left( 1+ a \;Ri  \right)   } + A_b^{vm}
    \end{aligned}    \right.
@@ -99,5 +99,5 @@
 is the maximum value that can be reached by the coefficient when $Ri\leq 0$, 
 $a=5$ and $n=2$. The last three values can be modified by setting the 
-\np{avmri}, \np{alp} and \np{nric} namelist parameters, respectively.
+\np{rn\_avmri}, \np{rn\_alp} and \np{nn\_ric} namelist parameters, respectively.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -107,24 +107,24 @@
 \label{ZDF_tke}
 
-%--------------------------------------------namtke---------------------------------------------------------
-\namdisplay{namtke}
+%--------------------------------------------namzdf_tke--------------------------------------------------
+\namdisplay{namzdf_tke}
 %--------------------------------------------------------------------------------------------------------------
 
 The vertical eddy viscosity and diffusivity coefficients are computed from a TKE 
-turbulent closure model based on a prognostic equation for $\bar {e}$, the turbulent 
+turbulent closure model based on a prognostic equation for $\bar{e}$, the turbulent 
 kinetic energy, and a closure assumption for the turbulent length scales. This 
 turbulent closure model has been developed by \citet{Bougeault1989} in the 
 atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and 
-embedded in OPA by \citet{Blanke1993} for equatorial Atlantic simulations. Since 
-then, significant modifications have been introduced by \citet{Madec1998} in both 
-the implementation and the formulation of the mixing length scale. The time 
-evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical 
-shear, its destruction through stratification, its vertical diffusion, and its dissipation 
-of \citet{Kolmogorov1942} type:
+embedded in OPA, the ancestor of NEMO, by \citet{Blanke1993} for equatorial Atlantic 
+simulations. Since then, significant modifications have been introduced by 
+\citet{Madec1998} in both the implementation and the formulation of the mixing 
+length scale. The time evolution of $\bar{e}$ is the result of the production of 
+$\bar{e}$ through vertical shear, its destruction through stratification, its vertical 
+diffusion, and its dissipation of \citet{Kolmogorov1942} type:
 \begin{equation} \label{Eq_zdftke_e}
 \frac{\partial \bar{e}}{\partial t} = 
-\frac{A^{vm}}{e_3 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
-							+\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
--A^{vT}\,N^2
+\frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
+				        +\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
+-K_\rho\,N^2
 +\frac{1}{e_3}	\;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 }
 				\;\frac{\partial \bar{e}}{\partial k}} \right]
@@ -133,15 +133,16 @@
 \begin{equation} \label{Eq_zdftke_kz}
    \begin{split}
-         A^{vm} &= C_k\  l_k\  \sqrt {\bar{e}}  	\\
-         A^{vT} &= A^{vm} / P_{rt}
+         K_m &= C_k\  l_k\  \sqrt {\bar{e}\; }  	\\
+         K_\rho &= A^{vm} / P_{rt}
    \end{split}
 \end{equation}
 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 
 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 
-$P_{rt} $ is the Prandtl number. The constants $C_k = \sqrt {2} /2$ and 
-$C_\epsilon = 0.1$ are designed to deal with vertical mixing at any depth 
-\citep{Gaspar1990}. They are set through namelist parameters \np{ediff} 
-and \np{ediss}. $P_{rt} $ can be set to unity or, following \citet{Blanke1993}, 
-be a function of the local Richardson number, $R_i $:
+$P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity 
+and diffusivity coefficients. The constants $C_k =  0.1$ and $C_\epsilon = \sqrt {2} /2$  
+$\approx 0.7$ are designed to deal with vertical mixing at any depth \citep{Gaspar1990}. 
+They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}. 
+$P_{rt}$ can be set to unity or, following \citet{Blanke1993}, be a function 
+of the local Richardson number, $R_i$:
 \begin{align*} \label{Eq_prt}
 P_{rt} = \begin{cases}
@@ -151,23 +152,22 @@
             \end{cases}
 \end{align*}
-Note that a horizontal Shapiro filter can optionally be applied to $R_i$. 
-However it is an obsolescent option that is not recommended.  
-The choice of $P_{rt} $ is controlled by the \np{npdl} namelist parameter.
+The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist parameter.
 
 For computational efficiency, the original formulation of the turbulent length 
 scales proposed by \citet{Gaspar1990} has been simplified. Four formulations 
-are proposed, the choice of which is controlled by the \np{nmxl} namelist 
+are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist 
 parameter. The first two are based on the following first order approximation 
 \citep{Blanke1993}:
 \begin{equation} \label{Eq_tke_mxl0_1}
-l_k = l_\epsilon = \sqrt {2 \bar e} / N
-\end{equation}
-which is valid in a stable stratified region with constant values of the brunt-
+l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N
+\end{equation}
+which is valid in a stable stratified region with constant values of the Brunt-
 Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance 
-to the surface or to the bottom (\np{nmxl}=0) or by the local vertical scale factor (\np{nmxl}=1). \citet{Blanke1993} notice that this simplification has two major 
+to the surface or to the bottom (\np{nn\_mxl} = 0) or by the local vertical scale factor 
+(\np{nn\_mxl} = 1). \citet{Blanke1993} notice that this simplification has two major 
 drawbacks: it makes no sense for locally unstable stratification and the 
 computation no longer uses all the information contained in the vertical density 
 profile. To overcome these drawbacks, \citet{Madec1998} introduces the 
-\np{nmxl}=2 or 3 cases, which add an extra assumption concerning the vertical 
+\np{nn\_mxl} = 2 or 3 cases, which add an extra assumption concerning the vertical 
 gradient of the computed length scale. So, the length scales are first evaluated 
 as in \eqref{Eq_tke_mxl0_1} and then bounded such that:
@@ -185,6 +185,5 @@
 constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, 
 the upward and downward length scales, and evaluate the dissipation and 
-mixing turbulent length scales as (and note that here we use numerical 
-indexing):
+mixing length scales as (and note that here we use numerical indexing):
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t] \label{Fig_mixing_length}  \begin{center}
@@ -196,16 +195,16 @@
 \begin{equation} \label{Eq_tke_mxl2}
 \begin{aligned}
- l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3T}^{(k)}\ \ \ \;  \right)
+ l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \;  \right)
     \quad &\text{ from $k=1$ to $jpk$ }\ \\
- l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3T}^{(k-1)}  \right)   
+ l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)}  \right)   
     \quad &\text{ from $k=jpk$ to $1$ }\ \\
 \end{aligned}
 \end{equation}
 where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1}, 
-$i.e.$ $l^{(k)} = \sqrt {2 \bar e^{(k)} / N^{(k)} }$.
-
-In the \np{nmxl}=2 case, the dissipation and mixing length scales take the same 
+$i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
+
+In the \np{nn\_mxl}=2 case, the dissipation and mixing length scales take the same 
 value: $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the 
-\np{nmxl}=2 case, the dissipation and mixing length scales are give 
+\np{nn\_mxl}=2 case, the dissipation and mixing turbulent length scales are give 
 as in \citet{Gaspar1990}:
 \begin{equation} \label{Eq_tke_mxl_gaspar}
@@ -217,18 +216,220 @@
 
 At the sea surface the value of $\bar{e}$ is prescribed from the wind 
-stress field: $\bar{e}=ebb\;\left| \tau \right|$ ($ebb=60$ by default) 
-with a minimal threshold of $emin0=10^{-4}~m^2.s^{-2}$ (namelist
+stress field: $\bar{e}=rn\_ebb\;\left| \tau \right|$ (\np{rn\_ebb}=60 by default) 
+with a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist
 parameters). Its value at the bottom of the ocean is assumed to be 
 equal to the value of the level just above. The time integration of the 
 $\bar{e}$ equation may formally lead to negative values because the 
 numerical scheme does not ensure its positivity. To overcome this 
-problem, a cut-off in the minimum value of $\bar{e}$ is used. Following 
-\citet{Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. 
-This allows the subsequent formulations to match that of\citet{Gargett1984} 
-for the diffusion in the thermocline and deep ocean :  $(A^{vT} = 10^{-3} / N)$. 
-In addition, a cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical 
+problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} 
+namelist parameter). Following \citet{Gaspar1990}, the cut-off value is set 
+to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations 
+to match that of \citet{Gargett1984} for the diffusion in the thermocline and 
+deep ocean :  $K_\rho = 10^{-3} / N$. 
+In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical 
 instabilities associated with too weak vertical diffusion. They must be 
 specified at least larger than the molecular values, and are set through 
-\textit{avm0} and \textit{avt0} (\textbf{namelist} parameters).
+\np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \S\ref{ZDF_cst}).
+
+% -------------------------------------------------------------------------------------------------------------
+%        TKE Turbulent Closure Scheme : new organization to energetic considerations
+% -------------------------------------------------------------------------------------------------------------
+\subsection{TKE discretization considerations (\key{zdftke})}
+\label{ZDF_tke_ene}
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t] \label{Fig_TKE_time_scheme}  \begin{center}
+\includegraphics[width=1.00\textwidth]{./TexFiles/Figures/Fig_ZDF_TKE_time_scheme.pdf}
+\caption {Illustration of the TKE time integration and its links to the momentum and tracer time integration. }
+\end{center}  
+\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+The production of turbulence by vertical shear (the first term of the right hand side 
+of \eqref{Eq_zdftke_e}) should balance the loss of kinetic energy associated with
+the vertical momentum diffusion (first line in \eqref{Eq_PE_zdf}). To do so a special care 
+have to be taken for both the time and space discretization of the TKE equation \citep{Burchard_OM02}.
+
+Let us first address the time stepping issue. Fig.~\ref{Fig_TKE_time_scheme} shows 
+how the two-level Leap-Frog time stepping of the momentum and tracer equations interplays 
+with the one-level forward time stepping of TKE equation. With this framework, the total loss 
+of kinetic energy (in 1D for the demonstration) due to the vertical momentum diffusion is 
+obtained by multiplying this quantity by $u^t$ and summing the result vertically:   
+\begin{equation} \label{Eq_energ1}
+\begin{split}
+\int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\
+&= \Bigl[  u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta}          
+ - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz }
+\end{split}
+\end{equation}
+Here, the vertical diffusion of momentum is discretized backward in time 
+with a coefficient, $K_m$, known at time $t$ (Fig.~\ref{Fig_TKE_time_scheme}), 
+as it is required when using the TKE scheme (see \S\ref{STP_forward_imp}). 
+The first term of the right hand side of \eqref{Eq_energ1} represents the kinetic energy 
+transfer at the surface (atmospheric forcing) and at the bottom (friction effect). 
+The second term is always negative. It is the dissipation rate of kinetic energy, 
+and thus minus the shear production rate of $\bar{e}$. \eqref{Eq_energ1} 
+implies that, to be energetically consistent, the production rate of $\bar{e}$ 
+used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as 
+${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$ (and not by the more straightforward 
+$K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$).
+
+A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification 
+(second term of the right hand side of \eqref{Eq_zdftke_e}). This term 
+must balance the input of potential energy resulting from vertical mixing. 
+The rate of change of potential energy (in 1D for the demonstration) due vertical 
+mixing is obtained by multiplying vertical density diffusion 
+tendency by $g\,z$ and and summing the result vertically:
+\begin{equation} \label{Eq_energ2}
+\begin{split}
+\int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\
+&= \Bigl[  g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta}  
+   - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz   \\
+&= - \Bigl[  z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta}
++ \int_{-H}^{\eta}{  \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz  }
+\end{split}
+\end{equation}
+where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. 
+The first term of the right hand side of \eqref{Eq_energ2}  is always zero 
+because there is no diffusive flux through the ocean surface and bottom). 
+The second term is minus the destruction rate of  $\bar{e}$ due to stratification. 
+Therefore \eqref{Eq_energ1} implies that, to be energetically consistent, the product 
+${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \eqref{Eq_zdftke_e}, the TKE equation.
+
+Let us now address the space discretization issue. 
+The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity 
+components are in the centre of the side faces of a $t$-box in staggered C-grid 
+(Fig.\ref{Fig_cell}). A space averaging is thus required to obtain the shear TKE production term.
+By redoing the \eqref{Eq_energ1} in the 3D case, it can be shown that the product of 
+eddy coefficient by the shear at $t$ and $t-\rdt$ must be performed prior to the averaging.
+Furthermore, the possible time variation of $e_3$ (\key{vvl} case) have to be taken into 
+account.
+
+The above energetic considerations leads to 
+the following final discrete form for the TKE equation:
+\begin{equation} \label{Eq_zdftke_ene}
+\begin{split}
+\frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv  
+\Biggl\{ \Biggr.
+  &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} } 
+                                                                              \ \frac{\delta_{k+1/2}[u^ t         ]}{{e_3u}^ t          }  \right) }^{\,i} \\
++&\overline{  \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} } 
+                                                                               \ \frac{\delta_{k+1/2}[v^ t         ]}{{e_3v}^ t          }  \right) }^{\,j} 
+\Biggr. \Biggr\}   \\
+%
+- &{K_\rho}^{t-\rdt}\,{(N^2)^t}    \\
+%
++&\frac{1}{{e_3w}^{t+\rdt}}  \;\delta_{k+1/2} \left[   {K_m}^{t-\rdt} \,\frac{\delta_{k}[(\bar{e})^{t+\rdt}]} {{e_3w}^{t+\rdt}}   \right]   \\
+%
+- &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt}
+\end{split}
+\end{equation}
+where the last two terms in \eqref{Eq_zdftke_ene} (vertical diffusion and Kolmogorov dissipation) 
+are time stepped using a backward scheme (see\S\ref{STP_forward_imp}). 
+Note that the Kolmogorov term has been linearized in time in order to render 
+the implicit computation possible. The restart of the TKE scheme 
+requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as they all appear in 
+the right hand side of \eqref{Eq_zdftke_ene}. For the latter, it is in fact 
+the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 
+
+% -------------------------------------------------------------------------------------------------------------
+%        GLS Generic Length Scale Scheme 
+% -------------------------------------------------------------------------------------------------------------
+\subsection{GLS Generic Length Scale (\key{zdfgls})}
+\label{ZDF_gls}
+
+%--------------------------------------------namgls---------------------------------------------------------
+\namdisplay{namgls}
+%--------------------------------------------------------------------------------------------------------------
+
+The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on 
+two prognostic equations: one for the turbulent kinetic energy $\bar {e}$, and another 
+for the generic length scale, $\psi$ \citep{Umlauf_Burchard_JMS03, Umlauf_Burchard_CSR05}. 
+This later variable is defined as : $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 
+where the triplet $(p, m, n)$ value given in Tab.\ref{Tab_GLS} allows to recover 
+a number of well-known turbulent closures ($k$-$kl$ \citep{Mellor_Yamada_1982}, 
+$k$-$\epsilon$ \citep{Rodi_1987}, $k$-$\omega$ \citep{Wilcox_1988} 
+among others \citep{Umlauf_Burchard_JMS03,Kantha_Carniel_CSR05}). 
+The GLS scheme is given by the following set of equations:
+\begin{equation} \label{Eq_zdfgls_e}
+\frac{\partial \bar{e}}{\partial t} = 
+\frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
+                                                   +\left( \frac{\partial v}{\partial k} \right)^2} \right]
+-K_\rho \,N^2
++\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right]
+- \epsilon
+\end{equation}
+
+\begin{equation} \label{Eq_zdfgls_psi}
+   \begin{split}
+\frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{
+\frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
+                                                                   +\left( \frac{\partial v}{\partial k} \right)^2} \right]
+- C_3 \,K_\rho\,N^2   - C_2 \,\epsilon \,Fw   \right\}             \\
+&+\frac{1}{e_3}  \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 }
+                                  \;\frac{\partial \psi}{\partial k}} \right]\;
+   \end{split}
+\end{equation}
+
+\begin{equation} \label{Eq_zdfgls_kz}
+   \begin{split}
+         K_m    &= C_{\mu} \ \sqrt {\bar{e}} \ l         \\
+         K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l
+   \end{split}
+\end{equation}
+
+\begin{equation} \label{Eq_zdfgls_eps}
+{\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \;
+\end{equation}
+where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}) 
+and $\epsilon$ the dissipation rate. 
+The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) 
+depends of the choice of the turbulence model. Four different turbulent models are pre-defined 
+(Tab.\ref{Tab_GLS}). They are made available through th \np{gls} namelist parameter. 
+
+%--------------------------------------------------TABLE--------------------------------------------------
+\begin{table}[htbp]  \label{Tab_GLS}
+\begin{center}
+%\begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c}
+\begin{tabular}{ccccc}
+          	             &   $k-kl$   & $k-\epsilon$ & $k-\omega$ &   generic   \\  
+%                        & \citep{Mellor_Yamada_1982} &  \citep{Rodi_1987}       & \citep{Wilcox_1988} &                 \\  
+\hline  \hline 
+\np{nn\_clo}     & \textbf{0} &   \textbf{1}  &   \textbf{2}   &    \textbf{3}   \\  
+\hline 
+$( p , n , m )$	       &   ( 0 , 1 , 1 )   & ( 3 , 1.5 , -1 )   & ( -1 , 0.5 , -1 )    &  ( 2 , 1 , -0.67 )  \\
+$\sigma_k$      &    2.44         &     1.              &      2.                &      0.8          \\
+$\sigma_\psi$  &    2.44         &     1.3            &      2.                 &       1.07       \\
+$C_1$              &      0.9         &     1.44          &      0.555          &       1.           \\
+$C_2$              &      0.5         &     1.92          &      0.833          &       1.22       \\
+$C_3$              &      1.           &     1.              &      1.                &       1.           \\
+$F_{wall}$        &      Yes        &       --             &     --                  &      --          \\
+\hline
+\hline
+\end{tabular}
+\caption {Set of predefined GLS parameters, or equivalently predefined turbulence models available with \key{gls} and controlled by the \np{nn\_clos} namelist parameter.}
+\end{center}
+\end{table}
+%--------------------------------------------------------------------------------------------------------------
+
+In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force
+the convergence of the mixing length towards $K\,z_b$ ($K$: Kappa and $z_b$: rugosity length) 
+value near physical boundaries (logarithmic boundary layer law). $C_{\mu}$ and $C_{\mu'}$ 
+are calculated from stability function proposed by \citet{Galperin_al_JAS88}, or by \citet{Kantha_Clayson_1994} 
+or one of the two functions suggested by \citet{Canuto_2001}  (\np{nn\_stab\_func} = 0, 1, 2 or 3, resp.}). The value of $C_{0\mu}$ depends of the choice of the stability function.
+
+The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated 
+thanks to Dirichlet or Neumann condition through \np{nn\_tkebc\_surf} and \np{nn\_tkebc\_bot}, resp. 
+The wave effect on the mixing could be also being considered \citep{Craig_Banner_1994}.
+
+The $\psi$ equation is known to fail in stably stratified flows, and for this reason 
+almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy. 
+With this clipping, the maximum permissible length scale is determined by 
+$l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$. A value of $c_{lim} = 0.53$ is often used 
+\citep{Galperin_al_JAS88}. \cite{Umlauf_Burchard_CSR05} show that the value of 
+the clipping factor is of crucial importance for the entrainment depth predicted in 
+stably stratified situations, and that its value has to be chosen in accordance 
+with the algebraic model for the turbulent ßuxes. The clipping is only activated 
+if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{clim\_galp} value.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -239,9 +440,8 @@
 
 %--------------------------------------------namkpp--------------------------------------------------------
-\namdisplay{namkpp}
-%--------------------------------------------------------------------------------------------------------------
-
-The K-Profile Parametrization (KKP) developed by \cite{Large_al_RG94} has been 
-implemented in \NEMO by J. Chanut (PhD reference to be added here!).
+\namdisplay{namzdf_kpp}
+%--------------------------------------------------------------------------------------------------------------
+
+The KKP scheme has been implemented by J. Chanut ...
 
 \colorbox{yellow}{Add a description of KPP here.}
@@ -274,8 +474,7 @@
 \label{ZDF_npc}
 
-%--------------------------------------------namnpc--------------------------------------------------------
-\namdisplay{namnpc}
-%--------------------------------------------------------------------------------------------------------------
-
+%--------------------------------------------namzdf--------------------------------------------------------
+\namdisplay{namzdf}
+%--------------------------------------------------------------------------------------------------------------
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -295,8 +494,8 @@
 
 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}=true. 
-It is applied at each \np{nnpc1} time step and mixes downwards instantaneously 
+It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously 
 the statically unstable portion of the water column, but only until the density 
 structure becomes neutrally stable ($i.e.$ until the mixed portion of the water 
-column has \textit{exactly} the density of the water just below) \citep{Madec1991a}. 
+column has \textit{exactly} the density of the water just below) \citep{Madec_al_JPO91}. 
 The associated algorithm is an iterative process used in the following way 
 (Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is 
@@ -320,5 +519,5 @@
 convective algorithm has been proved successful in studies of the deep 
 water formation in the north-western Mediterranean Sea 
-\citep{Madec1991a, Madec1991b, Madec1991c}.
+\citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}.
 
 Note that in the current implementation of this algorithm presents several 
@@ -341,5 +540,5 @@
 after a static adjustment. In that case, we recommend the addition of momentum 
 mixing in a manner that mimics the mixing in temperature and salinity 
-\citep{Speich1992, Speich1996}.
+\citep{Speich_PhD92, Speich_al_JPO96}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -347,5 +546,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Enhanced Vertical Diffusion (\np{ln\_zdfevd})]
-			{Enhanced Vertical Diffusion (\np{ln\_zdfevd}=.true.)}
+			{Enhanced Vertical Diffusion (\np{ln\_zdfevd}=true)}
 \label{ZDF_evd}
 
@@ -357,16 +556,20 @@
 In this case, the vertical eddy mixing coefficients are assigned very large values 
 (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable 
-($i.e.$ when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar1997, Lazar1999}. 
-This is done either on tracers only (\np{n\_evdm}=0) or on both momentum and 
-tracers (\np{n\_evdm}=1).
+($i.e.$ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) 
+\citep{Lazar_PhD97, Lazar_al_JPO99}. This is done either on tracers only 
+(\np{nn\_evdm}=0) or on both momentum and tracers (\np{nn\_evdm}=1).
 
 In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and 
-if \np{n\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 
-values also, are set equal to the namelist parameter \np{avevd}. A typical value 
-for $avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of 
+if \np{nn\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 
+values also, are set equal to the namelist parameter \np{rn\_avevd}. A typical value 
+for $rn\_avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of 
 convective processes is less time consuming than the convective adjustment 
 algorithm presented above when mixing both tracers and momentum in the 
 case of static instabilities. It requires the use of an implicit time stepping on 
 vertical diffusion terms (i.e. \np{ln\_zdfexp}=false). 
+
+Note that the stability test is performed on both \textit{before} and \textit{now} 
+values of $N^2$. This removes a potential source of divergence of odd and
+even time step in a leapfrog environment \citep{Leclair_PhD2010} (see \S\ref{STP_mLF}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -394,5 +597,5 @@
 The KPP turbulent closure scheme already includes enhanced vertical diffusion 
 in the case of convection, as governed by the variables $bvsqcon$ and $difcon$ 
-found in \mdl{zdfkpp}, therefore \np{np\_zdfevd} should not be used with the KPP 
+found in \mdl{zdfkpp}, therefore \np{ln\_zdfevd}=false should be used with the KPP 
 scheme. %gm%  + one word on non local flux with KPP scheme trakpp.F90 module...
 
@@ -404,6 +607,6 @@
 \label{ZDF_ddm}
 
-%-------------------------------------------namddm--------------------------------------------------------
-\namdisplay{namddm}
+%-------------------------------------------namzdf_ddm-------------------------------------------------
+\namdisplay{namzdf_ddm}
 %--------------------------------------------------------------------------------------------------------------
 
@@ -422,6 +625,5 @@
 \end{align*}
 where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection, 
-and $o$ by processes other than double diffusion. The rates of double-diffusive mixing 
-depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, 
+and $o$ by processes other than double diffusion. The rates of double-diffusive mixing depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, 
 where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline 
 contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt 
@@ -454,6 +656,5 @@
 we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.
 
-To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested 
-by Federov (1988) is used: 
+To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by Federov (1988) is used: 
 \begin{align} 	\label{Eq_zdfddm_d}
 A_d^{vT} &= 	\begin{cases}
@@ -481,6 +682,5 @@
 % Bottom Friction
 % ================================================================
-\section  [Bottom Friction (\textit{zdfbfr})]
-		{Bottom Friction (\mdl{zdfbfr} module)}
+\section  [Bottom Friction (\textit{zdfbfr})]   {Bottom Friction (\mdl{zdfbfr} module)}
 \label{ZDF_bfr}
 
@@ -493,9 +693,9 @@
 diffusive flux. For the bottom boundary layer, one has:
 \begin{equation} \label{Eq_zdfbfr_flux}
-A^{vm} \left( \partial \textbf{U}_h / \partial z \right) = \textbf{F}_h
-\end{equation}
-where $\textbf{F}_h$ is supposed to represent the horizontal momentum flux 
+A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U}
+\end{equation}
+where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum 
 outside the logarithmic turbulent boundary layer (thickness of the order of 
-1~m in the ocean). How $\textbf{F}_h$ influences the interior depends on the 
+1~m in the ocean). How ${\cal F}_h^{\textbf U}$ influences the interior depends on the 
 vertical resolution of the model near the bottom relative to the Ekman layer 
 depth. For example, in order to obtain an Ekman layer depth 
@@ -509,66 +709,78 @@
 this, consider the equation for $u$ at $k$, the last ocean level:
 \begin{equation} \label{Eq_zdfbfr_flux2}
-\frac{\partial u \; (k)}{\partial t} = \frac{1}{e_{3u}} \left[ A^{vm} \; (k) \frac{U \; (k-1) - U \; (k)}{e_{3uw} \; (k-1)} - F_u \right] \approx - \frac{F_u}{e_{3u}}
-\end{equation}
-For example, if the bottom layer thickness is 200~m, the Ekman transport will 
+\frac{\partial u_k}{\partial t} = \frac{1}{e_{3u}} \left[ \frac{A_{uw}^{vm}}{e_{3uw}} \delta_{k+1/2}\;[u] - {\cal F}^u_h \right] \approx - \frac{{\cal F}^u_{h}}{e_{3u}}
+\end{equation}
+If the bottom layer thickness is 200~m, the Ekman transport will 
 be distributed over that depth. On the other hand, if the vertical resolution 
 is high (1~m or less) and a turbulent closure model is used, the turbulent 
 Ekman layer will be represented explicitly by the model. However, the 
 logarithmic layer is never represented in current primitive equation model 
-applications: it is \emph{necessary} to parameterize the flux $\textbf{F}_h $. 
+applications: it is \emph{necessary} to parameterize the flux ${\cal F}^u_h $. 
 Two choices are available in \NEMO: a linear and a quadratic bottom friction. 
 Note that in both cases, the rotation between the interior velocity and the 
 bottom friction is neglected in the present release of \NEMO.
 
+In the code, the bottom friction is imposed by adding the trend due to the bottom 
+friction to the general momentum trend in \mdl{dynbfr}. For the time-split surface 
+pressure gradient algorithm, the momentum trend due to the barotropic component 
+needs to be handled separately. For this purpose it is convenient to compute and 
+store coefficients which can be simply combined with bottom velocities and geometric 
+values to provide the momentum trend due to bottom friction. 
+These coefficients are computed in \mdl{zdfbfr} and generally take the form 
+$c_b^{\textbf U}$ where:
+\begin{equation} \label{Eq_zdfbfr_bdef}
+\frac{\partial {\textbf U_h}}{\partial t} = 
+  - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b
+\end{equation}
+where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity.
+
 % -------------------------------------------------------------------------------------------------------------
 %       Linear Bottom Friction
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Linear Bottom Friction (\np{nbotfr} = 1) }
+\subsection{Linear Bottom Friction (\np{nn\_botfr} = 0 or 1) }
 \label{ZDF_bfr_linear}
 
-The linear bottom friction parameterisation assumes that the bottom friction 
-is proportional to the interior velocity (i.e. the velocity of the last model level):
+The linear bottom friction parameterisation (including the special case 
+of a free-slip condition) assumes that the bottom friction 
+is proportional to the interior velocity (i.e. the velocity of the last 
+model level):
 \begin{equation} \label{Eq_zdfbfr_linear}
-\textbf{F}_h = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b
-\end{equation}
-where $\textbf{U}_h^b$ is the horizontal velocity vector of the bottom ocean 
-layer and $r$ is a friction coefficient expressed in m.s$^{-1}$. This coefficient 
-is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 
+{\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b
+\end{equation}
+where $r$ is a friction coefficient expressed in ms$^{-1}$. 
+This coefficient is generally estimated by setting a typical decay time 
+$\tau$ in the deep ocean, 
 and setting $r = H / \tau$, where $H$ is the ocean depth. Commonly accepted 
-values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly1984}. 
+values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly_JMR84}. 
 A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used 
 in quasi-geostrophic models. One may consider the linear friction as an 
 approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$ (\citet{Gill1982}, 
 Eq. 9.6.6). For example, with a drag coefficient $C_D = 0.002$, a typical speed 
-of tidal currents of $U_{av} =0.1$~m.s$^{-1}$, and assuming an ocean depth 
-$H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m.s$^{-1}$. 
+of tidal currents of $U_{av} =0.1$~m\;s$^{-1}$, and assuming an ocean depth 
+$H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$. 
 This is the default value used in \NEMO. It corresponds to a decay time scale 
-of 115~days. It can be changed by specifying \np{bfric1} (namelist parameter).
-
-In the code, the bottom friction is imposed by updating the value of the 
-vertical eddy coefficient at the bottom level. Indeed, the discrete formulation 
-of (\ref{Eq_zdfbfr_linear}) at the last ocean $T-$level, using the fact that 
-$\textbf {U}_h =0$ below the ocean floor, leads to
-\begin{equation} \label{Eq_zdfbfr_linKz}
+of 115~days. It can be changed by specifying \np{rn\_bfric1} (namelist parameter).
+
+For the linear friction case the coefficients defined in the general 
+expression \eqref{Eq_zdfbfr_bdef} are: 
+\begin{equation} \label{Eq_zdfbfr_linbfr_b}
 \begin{split}
-A_u^{vm} &= r\;e_{3uw}\\
-A_v^{vm} &= r\;e_{3vw}\\
+ c_b^u &= - r\\
+ c_b^v &= - r\\
 \end{split}
 \end{equation}
-
-This update is done in \mdl{zdfbfr} when \np{nbotfr}=1. The value of $r$ 
-used is \np{bfric1}. Setting \np{nbotfr}=3 is equivalent to setting $r=0$ and 
-leads to a free-slip bottom boundary condition. Setting \np{nbotfr}=0 sets 
-$r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the background 
-vertical eddy coefficient, and a no-slip boundary condition is imposed. 
-Note that this latter choice generally leads to an underestimation of the 
-bottom friction: for example with a deepest level thickness of $200~m$ 
-and $A_{vb}^{\rm {\bf U}} =10^{-4}$m$^2$.s$^{-1}$, the friction coefficient 
-is only $r=10^{-6}$m.s$^{-1}$.
+When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfric1}. 
+Setting \np{nn\_botfr}=0 is equivalent to setting $r=0$ and leads to a free-slip 
+bottom boundary condition. These values are assigned in \mdl{zdfbfr}. 
+From v3.2 onwards there is support for local enhancement of these values 
+via an externally defined 2D mask array (\np{ln\_bfr2d}=true) given
+in the \ifile{bfr\_coef} input NetCDF file. The mask values should vary from 0 to 1. 
+Locations with a non-zero mask value will have the friction coefficient increased 
+by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfric1}.
 
 % -------------------------------------------------------------------------------------------------------------
 %       Non-Linear Bottom Friction
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Non-Linear Bottom Friction (\np{nbotfr} = 2)}
+\subsection{Non-Linear Bottom Friction (\np{nn\_botfr} = 2)}
 \label{ZDF_bfr_nonlinear}
 
@@ -576,31 +788,254 @@
 friction is quadratic: 
 \begin{equation} \label{Eq_zdfbfr_nonlinear}
-\textbf {F}_h = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 
+{\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 
 }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b 
 \end{equation}
-
-with $\textbf{U}_h^b = (u_b\;,\;v_b)$ the horizontal interior velocity ($i.e.$ 
-the horizontal velocity of the bottom ocean layer), $C_D$ a drag coefficient, 
-and $e_b $ a bottom turbulent kinetic energy due to tides, internal waves 
-breaking and other short time scale currents. A typical value of the drag 
-coefficient is $C_D = 10^{-3} $. As an example, the CME experiment 
-\citep{Treguier1992} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$, 
-while the FRAM experiment \citep{Killworth1992} uses $e_b =0$ 
-and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been 
-set as default values (\np{bfric2} and \np{bfeb2} namelist parameters).
+where $C_D$ is a drag coefficient, and $e_b $ a bottom turbulent kinetic energy 
+due to tides, internal waves breaking and other short time scale currents. 
+A typical value of the drag coefficient is $C_D = 10^{-3} $. As an example, 
+the CME experiment \citep{Treguier_JGR92} uses $C_D = 10^{-3}$ and 
+$e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992} 
+uses $C_D = 1.4\;10^{-3}$ and $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$. 
+The CME choices have been set as default values (\np{rn\_bfric2} and \np{rn\_bfeb2} 
+namelist parameters).
 
 As for the linear case, the bottom friction is imposed in the code by 
-updating the value of the vertical eddy coefficient at the bottom level:
-\begin{equation} \label{Eq_zdfbfr_nonlinKz}
+adding the trend due to the bottom friction to the general momentum trend 
+in \mdl{dynbfr}.
+For the non-linear friction case the terms
+computed in \mdl{zdfbfr}  are: 
+\begin{equation} \label{Eq_zdfbfr_nonlinbfr}
 \begin{split}
-A_u^{vm} &=C_D\; e_{3uw} \left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^
-{1/2}\\
-A_v^{vm} &=C_D\; e_{3uw} \left[  \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\
+ c_b^u &= - \; C_D\;\left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^{1/2}\\
+ c_b^v &= - \; C_D\;\left[  \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\
 \end{split}
 \end{equation}
 
-This update is done in \mdl{zdfbfr}. The coefficients that control the strength of the 
-non-linear bottom friction are initialized as namelist parameters: $C_D$= \np{bfri2}, 
-and $e_b$ =\np{bfeb2}.
-
-% ================================================================
+The coefficients that control the strength of the non-linear bottom friction are 
+initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}. 
+Note for applications which treat tides explicitly a low or even zero value of 
+\np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ 
+is possible via an externally defined 2D mask array (\np{ln\_bfr2d}=true). 
+See previous section for details.
+
+% -------------------------------------------------------------------------------------------------------------
+%       Bottom Friction stability
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Bottom Friction stability considerations}
+\label{ZDF_bfr_stability}
+
+Some care needs to exercised over the choice of parameters to ensure that the
+implementation of bottom friction does not induce numerical instability. For 
+the purposes of stability analysis, an approximation to \eqref{Eq_zdfbfr_flux2}
+is:
+\begin{equation} \label{Eqn_bfrstab}
+\begin{split}
+ \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\
+               &= -\frac{ru}{e_{3u}}\;2\rdt\\
+\end{split}
+\end{equation}
+\noindent where linear bottom friction and a leapfrog timestep have been assumed. 
+To ensure that the bottom friction cannot reverse the direction of flow it is necessary to have:
+\begin{equation}
+ |\Delta u| < \;|u| 
+\end{equation}
+\noindent which, using \eqref{Eqn_bfrstab}, gives:
+\begin{equation}
+r\frac{2\rdt}{e_{3u}} < 1 \qquad  \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\
+\end{equation}
+This same inequality can also be derived in the non-linear bottom friction case 
+if a velocity of 1 m.s$^{-1}$ is assumed. Alternatively, this criterion can be 
+rearranged to suggest a minimum bottom box thickness to ensure stability:
+\begin{equation}
+e_{3u} > 2\;r\;\rdt
+\end{equation}
+\noindent which it may be necessary to impose if partial steps are being used. 
+For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then
+$e_{3u}$ should be greater than 3.6 m. For most applications, with physically
+sensible parameters these restrictions should not be of concern. But 
+caution may be necessary if attempts are made to locally enhance the bottom
+friction parameters. 
+To ensure stability limits are imposed on the bottom friction coefficients both during 
+initialisation and at each time step. Checks at initialisation are made in \mdl{zdfbfr} 
+(assuming a 1 m.s$^{-1}$ velocity in the non-linear case).
+The number of breaches of the stability criterion are reported as well as the minimum 
+and maximum values that have been set. The criterion is also checked at each time step, 
+using the actual velocity, in \mdl{dynbfr}. Values of the bottom friction coefficient are 
+reduced as necessary to ensure stability; these changes are not reported.
+
+% -------------------------------------------------------------------------------------------------------------
+%       Bottom Friction with split-explicit time splitting
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Bottom Friction with split-explicit time splitting}
+\label{ZDF_bfr_ts}
+
+When calculating the momentum trend due to bottom friction in \mdl{dynbfr}, the
+bottom velocity at the before time step is used. This velocity includes both the
+baroclinic and barotropic components which is appropriate when using either the
+explicit or filtered surface pressure gradient algorithms (\key{dynspg\_exp} or 
+{\key{dynspg\_flt}). Extra attention is required, however, when using 
+split-explicit time stepping (\key{dynspg\_ts}). In this case the free surface 
+equation is solved with a small time step \np{nn\_baro}*\np{rn\_rdt}, while the three 
+dimensional prognostic variables are solved with a longer time step that is a 
+multiple of \np{rn\_rdt}. The trend in the barotropic momentum due to bottom 
+friction appropriate to this method is that given by the selected parameterisation 
+($i.e.$ linear or non-linear bottom friction) computed with the evolving velocities 
+at each barotropic timestep. 
+
+In the case of non-linear bottom friction, we have elected to partially linearise 
+the problem by keeping the coefficients fixed throughout the barotropic 
+time-stepping to those computed in \mdl{zdfbfr} using the now timestep. 
+This decision allows an efficient use of the $c_b^{\vect{U}}$ coefficients to:
+
+\begin{enumerate}
+\item On entry to \rou{dyn\_spg\_ts}, remove the contribution of the before
+barotropic velocity to the bottom friction component of the vertically
+integrated momentum trend. Note the same stability check that is carried out 
+on the bottom friction coefficient in \mdl{dynbfr} has to be applied here to
+ensure that the trend removed matches that which was added in \mdl{dynbfr}.
+\item At each barotropic step, compute the contribution of the current barotropic
+velocity to the trend due to bottom friction. Add this contribution to the
+vertically integrated momentum trend. This contribution is handled implicitly which
+eliminates the need to impose a stability criteria on the values of the bottom friction
+coefficient within the barotropic loop. 
+\end{enumerate}
+
+Note that the use of an implicit formulation
+for the bottom friction trend means that any limiting of the bottom friction coefficient 
+in \mdl{dynbfr} does not adversely affect the solution when using split-explicit time 
+splitting. This is because the major contribution to bottom friction is likely to come from 
+the barotropic component which uses the unrestricted value of the coefficient.
+
+The implicit formulation takes the form:
+\begin{equation} \label{Eq_zdfbfr_implicitts}
+ \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ]  
+\end{equation}
+where $\bar U$ is the barotropic velocity, $H_e$ is the full depth (including sea surface height), 
+$c_b^u$ is the bottom friction coefficient as calculated in \rou{zdf\_bfr} and $RHS$ represents 
+all the components to the vertically integrated momentum trend except for that due to bottom friction.
+
+
+
+
+% ================================================================
+% Tidal Mixing
+% ================================================================
+\section{Tidal Mixing}
+\label{ZDF_tmx}
+
+%--------------------------------------------namzdf_tmx--------------------------------------------------
+\namdisplay{namzdf_tmx}
+%--------------------------------------------------------------------------------------------------------------
+
+
+% -------------------------------------------------------------------------------------------------------------
+%        Bottom intensified tidal mixing 
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Bottom intensified tidal mixing}
+\label{ZDF_tmx_bottom}
+
+The parameterization of tidal mixing follows the general formulation for 
+the vertical eddy diffusivity proposed by \citet{St_Laurent_al_GRL02} and 
+first introduced in an OGCM by \citep{Simmons_al_OM04}. 
+In this formulation an additional vertical diffusivity resulting from internal tide breaking, 
+$A^{vT}_{tides}$ is expressed as a function of $E(x,y)$, the energy transfer from barotropic 
+tides to baroclinic tides : 
+\begin{equation} \label{Eq_Ktides}
+A^{vT}_{tides} =  q \,\Gamma \,\frac{ E(x,y) \, F(z) }{ \rho \, N^2 }
+\end{equation}
+where $\Gamma$ is the mixing efficiency, $N$ the Brunt-Vais\"{a}l\"{a} frequency 
+(see \S\ref{TRA_bn2}), $\rho$ the density, $q$ the tidal dissipation efficiency, 
+and $F(z)$ the vertical structure function. 
+
+The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn\_me} namelist parameter)
+and is usually taken to be the canonical value of $\Gamma = 0.2$ (Osborn 1980). 
+The tidal dissipation efficiency is given by the parameter $q$ (\np{rn\_tfe} namelist parameter) 
+represents the part of the internal wave energy flux $E(x, y)$ that is dissipated locally, 
+with the remaining $1-q$ radiating away as low mode internal waves and 
+contributing to the background internal wave field. A value of $q=1/3$ is typically used  
+\citet{St_Laurent_al_GRL02}.
+The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical. 
+It is implemented as a simple exponential decaying upward away from the bottom, 
+with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter, with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04}, 
+\begin{equation} \label{Eq_Fz}
+F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) }
+\end{equation}
+and is normalized so that vertical integral over the water column is unity. 
+
+The associated vertical viscosity is calculated from the vertical 
+diffusivity assuming a Prandtl number of 1, $i.e.$ $A^{vm}_{tides}=A^{vT}_{tides}$. 
+In the limit of $N \rightarrow 0$ (or becoming negative), the vertical diffusivity 
+is capped at $300\,cm^2/s$ and impose a lower limit on $N^2$ of \np{rn\_n2min} 
+usually set to $10^{-8} s^{-2}$. These bounds are usually rarely encountered.
+
+The internal wave energy map, $E(x, y)$ in \eqref{Eq_Ktides}, is derived 
+from a barotropic model of the tides utilizing a parameterization of the 
+conversion of barotropic tidal energy into internal waves. 
+The essential goal of the parameterization is to represent the momentum 
+exchange between the barotropic tides and the unrepresented internal waves 
+induced by the tidal ßow over rough topography in a stratified ocean. 
+In the current version of \NEMO, the map is built from the output of 
+the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}.
+This model provides the dissipation associated with internal wave energy for the M2 and K1 
+tides component (Fig.~\ref{Fig_ZDF_M2_K1_tmx}). The S2 dissipation is simply approximated
+as being $1/4$ of the M2 one. The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$. 
+Its global mean value is $1.1$ TW, in agreement with independent estimates 
+\citep{Egbert_Ray_Nat00, Egbert_Ray_JGR01}. 
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t] \label{Fig_ZDF_M2_K1_tmx}  \begin{center}
+\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_ZDF_M2_K1_tmx.pdf}
+\caption {(a) M2 and (b) K2 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$). }
+\end{center}  
+\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
+ 
+% -------------------------------------------------------------------------------------------------------------
+%        Indonesian area specific treatment 
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Indonesian area specific treatment}
+\label{ZDF_tmx_itf}
+
+When the Indonesian Through Flow (ITF) area is included in the model domain,
+a specific treatment of tidal induced mixing in this area can be used. 
+It is activated through the namelist logical \np{ln\_tmx\_itf}, and the user must provide
+an input NetCDF file, \ifile{mask\_itf}, which contains a mask array defining the ITF area
+where the specific treatment is applied.
+
+When \np{ln\_tmx\_itf}=true, the two key parameters $q$ and $F(z)$ are adjusted following 
+the parameterisation developed by \ref{Koch-Larrouy_al_GRL07}:
+
+First, the Indonesian archipelago is a complex geographic region 
+with a series of large, deep, semi-enclosed basins connected via 
+numerous narrow straits. Once generated, internal tides remain 
+confined within this semi-enclosed area and hardly radiate away. 
+Therefore all the internal tides energy is consumed within this area. 
+So it is assumed that $q = 1$, $i.e.$ all the energy generated is available for mixing.
+Note that for test purposed, the ITF tidal dissipation efficiency is a 
+namelist parameter (\np{rn\_tfe\_itf}). A value of $1$ or close to is
+this recommended for this parameter.
+
+Second, the vertical structure function, $F(z)$, is no more associated
+with a bottom intensification of the mixing, but with a maximum of 
+energy available within the thermocline. \ref{Koch-Larrouy_al_GRL07} 
+have suggested that the vertical distribution of the energy dissipation 
+proportional to $N^2$ below the core of the thermocline and to $N$ above. 
+The resulting $F(z)$ is:
+\begin{equation} \label{Eq_Fz_itf}
+F(i,j,k) \sim     \left\{ \begin{aligned}
+\frac{q\,\Gamma E(i,j) } {\rho N \, \int N     dz}    \qquad \text{when $\partial_z N < 0$} \\
+\frac{q\,\Gamma E(i,j) } {\rho     \, \int N^2 dz}    \qquad \text{when $\partial_z N > 0$}
+                      \end{aligned} \right.
+\end{equation}
+
+Averaged over the ITF area, the resulting tidal mixing coefficient is $1.5\,cm^2/s$, 
+which agrees with the independent estimates inferred from observations. 
+Introduced in a regional OGCM, the parameterization improves the water mass 
+characteristics in the different Indonesian seas, suggesting that the horizontal 
+and vertical distributions of the mixing are adequately prescribed 
+\citep{Koch-Larrouy_al_GRL07, Koch-Larrouy_al_OD08a, Koch-Larrouy_al_OD08b}.
+Note also that such a parameterisation has a sugnificant impact on the behaviour 
+of global coupled GCMs \citep{Koch-Larrouy_al_CD10}.
+
+
+% ================================================================
