Index: branches/2015/nemo_v3_6_STABLE/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- branches/2015/nemo_v3_6_STABLE/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 5550)
+++ branches/2015/nemo_v3_6_STABLE/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 6275)
@@ -33,10 +33,9 @@
 points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These 
 coefficients can be assumed to be either constant, or a function of the local 
-Richardson number, or computed from a turbulent closure model (either 
-TKE or KPP formulation). The computation of these coefficients is initialized 
-in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or 
-\mdl{zdfkpp} modules. The trends due to the vertical momentum and tracer 
-diffusion, including the surface forcing, are computed and added to the 
-general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 
+Richardson number, or computed from a turbulent closure model (TKE, GLS or KPP formulation). 
+The computation of these coefficients is initialized in the \mdl{zdfini} module 
+and performed in the \mdl{zdfric}, \mdl{zdftke}, \mdl{zdfgls} or \mdl{zdfkpp} modules. 
+The trends due to the vertical momentum and tracer diffusion, including the surface forcing, 
+are computed and added to the 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{ln\_zdfexp}=true) or a backward time stepping 
@@ -355,9 +354,46 @@
 %--------------------------------------------------------------%
 
-To be add here a description of "penetration of TKE" and the associated namelist parameters
- \np{nn\_etau}, \np{rn\_efr} and \np{nn\_htau}.
+Vertical mixing parameterizations commonly used in ocean general circulation models 
+tend to produce mixed-layer depths that are too shallow during summer months and windy conditions.
+This bias is particularly acute over the Southern Ocean. 
+To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme  \cite{Rodgers_2014}. 
+The parameterization is an empirical one, $i.e.$ not derived from theoretical considerations, 
+but rather is meant to account for observed processes that affect the density structure of 
+the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme 
+($i.e.$ near-inertial oscillations and ocean swells and waves).
+
+When using this parameterization ($i.e.$ when \np{nn\_etau}~=~1), the TKE input to the ocean ($S$) 
+imposed by the winds in the form of near-inertial oscillations, swell and waves is parameterized 
+by \eqref{ZDF_Esbc} the standard TKE surface boundary condition, plus a depth depend one given by:
+\begin{equation}  \label{ZDF_Ehtau}
+S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau} 
+\end{equation}
+where 
+$z$ is the depth,  
+$e_s$ is TKE surface boundary condition, 
+$f_r$ is the fraction of the surface TKE that penetrate in the ocean, 
+$h_\tau$ is a vertical mixing length scale that controls exponential shape of the penetration, 
+and $f_i$ is the ice concentration (no penetration if $f_i=1$, that is if the ocean is entirely 
+covered by sea-ice).
+The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter. 
+The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}~=~0) 
+or a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m 
+at high latitudes (\np{nn\_etau}~=~1). 
+
+Note that two other option existe, \np{nn\_etau}~=~2, or 3. They correspond to applying 
+\eqref{ZDF_Ehtau} only at the base of the mixed layer, or to using the high frequency part 
+of the stress to evaluate the fraction of TKE that penetrate the ocean. 
+Those two options are obsolescent features introduced for test purposes.
+They will be removed in the next release. 
+
+
 
 % from Burchard et al OM 2008 : 
-% the most critical process not reproduced by statistical turbulence models is the activity of internal waves and their interaction with turbulence. After the Reynolds decomposition, internal waves are in principle included in the RANS equations, but later partially excluded by the hydrostatic assumption and the model resolution. Thus far, the representation of internal wave mixing in ocean models has been relatively crude (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
+% the most critical process not reproduced by statistical turbulence models is the activity of 
+% internal waves and their interaction with turbulence. After the Reynolds decomposition, 
+% internal waves are in principle included in the RANS equations, but later partially 
+% excluded by the hydrostatic assumption and the model resolution. 
+% Thus far, the representation of internal wave mixing in ocean models has been relatively crude 
+% (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
 
 
@@ -586,5 +622,6 @@
 Options are defined through the  \ngn{namzdf\_kpp} namelist variables.
 
-\colorbox{yellow}{Add a description of KPP here.}
+Note that KPP is an obsolescent feature of the \NEMO system. 
+It will be removed in the next release (v3.7 and followings).
 
 
@@ -636,5 +673,5 @@
 
 Options are defined through the  \ngn{namzdf} namelist variables.
-The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}=true. 
+The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}~=~\textit{true}. 
 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 
@@ -644,5 +681,5 @@
 (Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is 
 found. Assume in the following that the instability is located between levels 
-$k$ and $k+1$. The potential temperature and salinity in the two levels are 
+$k$ and $k+1$. The temperature and salinity in the two levels are 
 vertically mixed, conserving the heat and salt contents of the water column. 
 The new density is then computed by a linear approximation. If the new 
@@ -664,24 +701,14 @@
 \citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}.
 
-Note that in the current implementation of this algorithm presents several 
-limitations. First, potential density referenced to the sea surface is used to 
-check whether the density profile is stable or not. This is a strong 
-simplification which leads to large errors for realistic ocean simulations. 
-Indeed, many water masses of the world ocean, especially Antarctic Bottom
-Water, are unstable when represented in surface-referenced potential density. 
-The scheme will erroneously mix them up. Second, the mixing of potential 
-density is assumed to be linear. This assures the convergence of the algorithm 
-even when the equation of state is non-linear. Small static instabilities can thus 
-persist due to cabbeling: they will be treated at the next time step. 
-Third, temperature and salinity, and thus density, are mixed, but the 
-corresponding velocity fields remain unchanged. When using a Richardson 
-Number dependent eddy viscosity, the mixing of momentum is done through 
-the vertical diffusion: after a static adjustment, the Richardson Number is zero 
-and thus the eddy viscosity coefficient is at a maximum. When this convective 
-adjustment algorithm is used with constant vertical eddy viscosity, spurious 
-solutions can occur since the vertical momentum diffusion remains small even 
-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{Speich_PhD92, Speich_al_JPO96}.
+The current implementation has been modified in order to deal with any non linear 
+equation of seawater (L. Brodeau, personnal communication). 
+Two main differences have been introduced compared to the original algorithm: 
+$(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency 
+(not the the difference in potential density) ; 
+$(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients 
+are vertically mixed in the same way their temperature and salinity has been mixed.
+These two modifications allow the algorithm to perform properly and accurately 
+with TEOS10 or EOS-80 without having to recompute the expansion coefficients at each 
+mixing iteration.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -689,5 +716,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}
 
