Index: trunk/DOC/TexFiles/Chapters/Chap_DIA.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_DIA.tex	(revision 6409)
+++ trunk/DOC/TexFiles/Chapters/Chap_DIA.tex	(revision 6497)
@@ -1409,30 +1409,4 @@
 
 % -------------------------------------------------------------------------------------------------------------
-%       25 hour mean and hourly Surface, Mid and Bed 
-% -------------------------------------------------------------------------------------------------------------
-\section{25 hour mean output for tidal models }
-
-A module is available to compute a crudely detided M2 signal by obtaining a 25 hour mean.
-The 25 hour mean is available for daily runs by summing up the 25 hourly instantananeous hourly values from
-midnight at the start of the day to midight at the day end.
-This diagnostic is actived with the logical  $ln\_dia25h$
-
-%------------------------------------------nam_dia25h------------------------------------------------------
-\namdisplay{nam_dia25h}
-%----------------------------------------------------------------------------------------------------------
-
-\section{Top Middle and Bed hourly output }
-
-A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables. 
-This can be a useful diagnostic when hourly or sub-hourly output is required in high resolution tidal outputs.
-The tidal signal is retained but the overall data usage is cut to just three vertical levels. Also the bottom level 
-is calculated for each cell.
-This diagnostic is actived with the logical  $ln\_diatmb$
-
-%------------------------------------------nam_diatmb-----------------------------------------------------
-\namdisplay{nam_diatmb}
-%----------------------------------------------------------------------------------------------------------
-
-% -------------------------------------------------------------------------------------------------------------
 %       Sections transports
 % -------------------------------------------------------------------------------------------------------------
@@ -1440,4 +1414,8 @@
 \label{DIA_diag_dct}
 
+%------------------------------------------namdct----------------------------------------------------
+\namdisplay{namdct}
+%-------------------------------------------------------------------------------------------------------------
+
 A module is available to compute the transport of volume, heat and salt through sections. 
 This diagnostic is actived with \key{diadct}.
@@ -1459,8 +1437,4 @@
 and the time scales over which they are averaged, as well as the level of output for debugging:
 
-%------------------------------------------namdct----------------------------------------------------
-\namdisplay{namdct}
-%-------------------------------------------------------------------------------------------------------------
-
 \np{nn\_dct}: frequency of instantaneous transports computing
 
@@ -1469,7 +1443,7 @@
 \np{nn\_debug}: debugging of the section
 
-\subsubsection{ To create a binary file containing the pathway of each section }
-
-In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, the file \texttt{ {list\_sections.ascii\_global}}
+\subsubsection{ Creating a binary file containing the pathway of each section }
+
+In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, the file \textit{ {list\_sections.ascii\_global}}
 contains a list of all the sections that are to be computed (this list of sections is based on MERSEA project metrics).
 
@@ -1583,68 +1557,4 @@
 \texttt{=/0, =/ 1000.}   &  diagonal   & eastward  & westward  & postive: eastward  \\ \hline                
 \end{tabular}
-
-
-
-% -------------------------------------------------------------------------------------------------------------
-%       Other Diagnostics
-% -------------------------------------------------------------------------------------------------------------
-\section{Other Diagnostics (\key{diahth}, \key{diaar5})}
-\label{DIA_diag_others}
-
-
-Aside from the standard model variables, other diagnostics can be computed 
-on-line. The available ready-to-add diagnostics routines can be found in directory DIA. 
-Among the available diagnostics the following ones are obtained when defining 
-the \key{diahth} CPP key: 
-
-- the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
-
-- the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth})
-
-- the depth of the 20\deg C isotherm (\mdl{diahth})
-
-- the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth})
-
-The poleward heat and salt transports, their advective and diffusive component, and 
-the meriodional stream function can be computed on-line in \mdl{diaptr} 
-\np{ln\_diaptr} to true (see the \textit{\ngn{namptr} } namelist below).  
-When \np{ln\_subbas}~=~true, transports and stream function are computed 
-for the Atlantic, Indian, Pacific and Indo-Pacific Oceans (defined north of 30\deg S) 
-as well as for the World Ocean. The sub-basin decomposition requires an input file 
-(\ifile{subbasins}) which contains three 2D mask arrays, the Indo-Pacific mask 
-been deduced from the sum of the Indian and Pacific mask (Fig~\ref{Fig_mask_subasins}). 
-
-%------------------------------------------namptr----------------------------------------------------
-\namdisplay{namptr} 
-%-------------------------------------------------------------------------------------------------------------
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] 	  \begin{center}
-\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_mask_subasins.pdf}
-\caption{	\label{Fig_mask_subasins}
-Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute
-the heat and salt transports as well as the meridional stream-function: Atlantic basin (red), 
-Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green). 
-Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay 
-are removed from the sub-basins. Note also that the Arctic Ocean has been split 
-into Atlantic and Pacific basins along the North fold line.  }
-\end{center}   \end{figure}  
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-In addition, a series of diagnostics has been added in the \mdl{diaar5}. 
-They corresponds to outputs that are required for AR5 simulations 
-(see Section \ref{DIA_steric} below for one of them). 
-Activating those outputs requires to define the \key{diaar5} CPP key.
-\\
-\\
-
-\section{Courant numbers}
-Courant numbers provide a theoretical indication of the model's numerical stability. The advective Courant numbers can be calculated according to
-\begin{equation}
-\label{eq:CFL}
-C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}}
-\end{equation}
-in the zonal, meridional and vertical directions respectively. The vertical component is included although it is not strictly valid as the vertical velocity is calculated from the continuity equation rather than as a prognostic variable. Physically this represents the rate at which information is propogated across a grid cell. Values greater than 1 indicate that information is propagated across more than one grid cell in a single time step.
-
-The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist. The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii. In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of where the maximum value occurs. At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along with the coordinates of each. The maximum values from the run are also copied to the ocean.output file. 
 
 
@@ -1802,4 +1712,119 @@
 the \key{diaar5} defined to be called.
 
+
+
+% -------------------------------------------------------------------------------------------------------------
+%       Other Diagnostics
+% -------------------------------------------------------------------------------------------------------------
+\section{Other Diagnostics (\key{diahth}, \key{diaar5})}
+\label{DIA_diag_others}
+
+
+Aside from the standard model variables, other diagnostics can be computed on-line. 
+The available ready-to-add diagnostics modules can be found in directory DIA. 
+
+\subsection{Depth of various quantities (\mdl{diahth})}
+
+Among the available diagnostics the following ones are obtained when defining 
+the \key{diahth} CPP key: 
+
+- the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
+
+- the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth})
+
+- the depth of the 20\deg C isotherm (\mdl{diahth})
+
+- the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth})
+
+% -----------------------------------------------------------
+%     Poleward heat and salt transports
+% -----------------------------------------------------------
+
+\subsection{Poleward heat and salt transports (\mdl{diaptr})}
+
+%------------------------------------------namptr-----------------------------------------
+\namdisplay{namptr} 
+%-----------------------------------------------------------------------------------------
+
+The poleward heat and salt transports, their advective and diffusive component, and 
+the meriodional stream function can be computed on-line in \mdl{diaptr} 
+\np{ln\_diaptr} to true (see the \textit{\ngn{namptr} } namelist below).  
+When \np{ln\_subbas}~=~true, transports and stream function are computed 
+for the Atlantic, Indian, Pacific and Indo-Pacific Oceans (defined north of 30\deg S) 
+as well as for the World Ocean. The sub-basin decomposition requires an input file 
+(\ifile{subbasins}) which contains three 2D mask arrays, the Indo-Pacific mask 
+been deduced from the sum of the Indian and Pacific mask (Fig~\ref{Fig_mask_subasins}). 
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t] 	  \begin{center}
+\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_mask_subasins.pdf}
+\caption{	\label{Fig_mask_subasins}
+Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute
+the heat and salt transports as well as the meridional stream-function: Atlantic basin (red), 
+Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green). 
+Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay 
+are removed from the sub-basins. Note also that the Arctic Ocean has been split 
+into Atlantic and Pacific basins along the North fold line.  }
+\end{center}   \end{figure}  
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+
+% -----------------------------------------------------------
+%       CMIP specific diagnostics 
+% -----------------------------------------------------------
+\subsection{CMIP specific diagnostics (\mdl{diaar5})}
+
+A series of diagnostics has been added in the \mdl{diaar5}. 
+They corresponds to outputs that are required for AR5 simulations (CMIP5)
+(see also Section \ref{DIA_steric} for one of them). 
+Activating those outputs requires to define the \key{diaar5} CPP key.
+
+
+% -----------------------------------------------------------
+%       25 hour mean and hourly Surface, Mid and Bed 
+% -----------------------------------------------------------
+\subsection{25 hour mean output for tidal models }
+
+%------------------------------------------nam_dia25h-------------------------------------
+\namdisplay{nam_dia25h}
+%-----------------------------------------------------------------------------------------
+
+A module is available to compute a crudely detided M2 signal by obtaining a 25 hour mean.
+The 25 hour mean is available for daily runs by summing up the 25 hourly instantananeous hourly values from
+midnight at the start of the day to midight at the day end.
+This diagnostic is actived with the logical  $ln\_dia25h$
+
+
+% -----------------------------------------------------------
+%     Top Middle and Bed hourly output
+% -----------------------------------------------------------
+\subsection{Top Middle and Bed hourly output }
+
+%------------------------------------------nam_diatmb-----------------------------------------------------
+\namdisplay{nam_diatmb}
+%----------------------------------------------------------------------------------------------------------
+
+A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables. 
+This can be a useful diagnostic when hourly or sub-hourly output is required in high resolution tidal outputs.
+The tidal signal is retained but the overall data usage is cut to just three vertical levels. Also the bottom level 
+is calculated for each cell.
+This diagnostic is actived with the logical  $ln\_diatmb$
+
+
+
+% -----------------------------------------------------------
+%     Courant numbers
+% -----------------------------------------------------------
+\subsection{Courant numbers}
+Courant numbers provide a theoretical indication of the model's numerical stability. The advective Courant numbers can be calculated according to
+\begin{equation}
+\label{eq:CFL}
+C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}}
+\end{equation}
+in the zonal, meridional and vertical directions respectively. The vertical component is included although it is not strictly valid as the vertical velocity is calculated from the continuity equation rather than as a prognostic variable. Physically this represents the rate at which information is propogated across a grid cell. Values greater than 1 indicate that information is propagated across more than one grid cell in a single time step.
+
+The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist. The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii. In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of where the maximum value occurs. At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along with the coordinates of each. The maximum values from the run are also copied to the ocean.output file. 
+
+
 % ================================================================
 
Index: trunk/DOC/TexFiles/Chapters/Chap_DOM.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 6409)
+++ trunk/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 6497)
@@ -486,5 +486,6 @@
 The last choice in terms of vertical coordinate concerns the presence (or not) in the model domain 
 of ocean cavities beneath ice shelves. Setting \np{ln\_isfcav} to true allows to manage ocean cavities, 
-otherwise they are filled in.
+otherwise they are filled in. This option is currently only available in $z$- or $zps$-coordinate,
+and partial step are also applied at the ocean/ice shelf interface. 
 
 Contrary to the horizontal grid, the vertical grid is computed in the code and no 
@@ -772,7 +773,7 @@
 \end{equation}
 
-where $s_{min}$ is the depth at which the s-coordinate stretching starts and 
-allows a z-coordinate to placed on top of the stretched coordinate, 
-and z is the depth (negative down from the asea surface).
+where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and 
+allows a $z$-coordinate to placed on top of the stretched coordinate, 
+and $z$ is the depth (negative down from the asea surface).
 
 \begin{equation}
@@ -886,8 +887,9 @@
 that do not communicate with another ocean point at the same level are eliminated.
 
-In case of ice shelf cavities, as for the representation of bathymetry, a 2D integer array, misfdep, is created. 
-misfdep defines the level of the first wet $t$-point (ie below the ice-shelf/ocean interface). All the cells between $k=1$ and $misfdep(i,j)-1$ are masked. 
-By default, $misfdep(:,:)=1$ and no cells are masked.
-Modifications of the model bathymetry and ice shelf draft into 
+As for the representation of bathymetry, a 2D integer array, misfdep, is created. 
+misfdep defines the level of the first wet $t$-point. All the cells between $k=1$ and $misfdep(i,j)-1$ are masked. 
+By default, misfdep(:,:)=1 and no cells are masked.
+
+In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into 
 the cavities are performed in the \textit{zgr\_isf} routine. The compatibility between ice shelf draft and bathymetry is checked. 
 All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded ($i.e.$ masked). 
@@ -903,10 +905,12 @@
 vmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j+1,k)	\\
 fmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)	\\
-                   & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\
+             &    \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\
 wmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j,k-1) \text{ with } wmask(i,j,1) = tmask(i,j,1) 
 \end{align*}
 
-Note, wmask is now defined. It allows, in case of ice shelves, 
-to deal with the top boundary (ice shelf/ocean interface) exactly in the same way as for the bottom boundary. 
+Note that, without ice shelves cavities, masks at $t-$ and $w-$points are identical with 
+the numerical indexing used (\S~\ref{DOM_Num_Index}). Nevertheless, $wmask$ are required 
+with oceean cavities to deal with the top boundary (ice shelf/ocean interface) 
+exactly in the same way as for the bottom boundary. 
 
 The specification of closed lateral boundaries requires that at least the first and last 
@@ -916,7 +920,4 @@
 (and so too the mask arrays) (see \S~\ref{LBC_jperio}).
 
-%%%
-\gmcomment{   \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}.  }
-%%%
 
 % ================================================================
Index: trunk/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 6409)
+++ trunk/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 6497)
@@ -128,7 +128,7 @@
 The ocean model provides, at each time step, to the surface module (\mdl{sbcmod}) 
 the surface currents, temperature and salinity.  
-These variables are averaged over \np{nf\_sbc} time-step (\ref{Tab_ssm}), 
+These variables are averaged over \np{nn\_fsbc} time-step (\ref{Tab_ssm}), 
 and it is these averaged fields which are used to computes the surface fluxes 
-at a frequency of \np{nf\_sbc} time-step.
+at a frequency of \np{nn\_fsbc} time-step.
 
 
@@ -144,6 +144,6 @@
 \caption{  \label{Tab_ssm}   
 Ocean variables provided by the ocean to the surface module (SBC). 
-The variable are averaged over nf{\_}sbc time step, $i.e.$ the frequency of 
-computation of surface fluxes.}
+The variable are averaged over nn{\_}fsbc time step, 
+$i.e.$ the frequency of computation of surface fluxes.}
 \end{center}   \end{table}
 %--------------------------------------------------------------------------------------------------------------
@@ -592,4 +592,20 @@
 or larger than the one of the input atmospheric fields.
 
+The \np{sn\_wndi}, \np{sn\_wndj}, \np{sn\_qsr}, \np{sn\_qlw}, \np{sn\_tair}, \np{sn\_humi},
+\np{sn\_prec}, \np{sn\_snow}, \np{sn\_tdif} parameters describe the fields 
+and the way they have to be used (spatial and temporal interpolations). 
+
+\np{cn\_dir} is the directory of location of bulk files
+\np{ln\_taudif} is the flag to specify if we use Hight Frequency (HF) tau information (.true.) or not (.false.)
+\np{rn\_zqt}: is the height of humidity and temperature measurements (m)
+\np{rn\_zu}: is the height of wind measurements (m)
+
+Three multiplicative factors are availables : 
+\np{rn\_pfac} and \np{rn\_efac} allows to adjust (if necessary) the global freshwater budget 
+by increasing/reducing the precipitations (total and snow) and or evaporation, respectively.
+The third one,\np{rn\_vfac}, control to which extend the ice/ocean velocities are taken into account 
+in the calculation of surface wind stress. Its range should be between zero and one, 
+and it is recommended to set it to 0.
+
 % -------------------------------------------------------------------------------------------------------------
 %        CLIO Bulk formulea
@@ -926,5 +942,6 @@
 \begin{description}
 \item[\np{nn\_isf}~=~1]
-The ice shelf cavity is represented (\np{ln\_isfcav}~=~true needed). The fwf and heat flux are computed. Two different bulk formula are available:
+The ice shelf cavity is represented (\np{ln\_isfcav}~=~true needed). The fwf and heat flux are computed. 
+Two different bulk formula are available:
 	\begin{description}
 	\item[\np{nn\_isfblk}~=~1]
@@ -988,7 +1005,7 @@
 This parameter is only used if \np{nn\_isf}~=~1 or \np{nn\_isf}~=~4
 
-If \np{rn\_hisf\_tbl} = 0.0, the fluxes are put in the top level whatever is its tickness. 
-
-If \np{rn\_hisf\_tbl} $>$ 0.0, the fluxes are spread over the first \np{rn\_hisf\_tbl} m (ie over one or several cells).\\
+If \np{rn\_hisf\_tbl} = 0., the fluxes are put in the top level whatever is its tickness. 
+
+If \np{rn\_hisf\_tbl} $>$ 0., the fluxes are spread over the first \np{rn\_hisf\_tbl} m (ie over one or several cells).\\
 
 The ice shelf melt is implemented as a volume flux with in the same way as for the runoff.
Index: trunk/DOC/TexFiles/Chapters/Chap_STO.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_STO.tex	(revision 6409)
+++ trunk/DOC/TexFiles/Chapters/Chap_STO.tex	(revision 6497)
@@ -5,14 +5,169 @@
 \label{STO}
 
+Authors: P.-A. Bouttier
+
 \minitoc
 
+\newpage
 
-\newpage
-$\ $\newline    % force a new line
+
+The stochastic parametrization module aims to explicitly simulate uncertainties in the model. 
+More particularly, \cite{Brankart_OM2013} has shown that, 
+because of the nonlinearity of the seawater equation of state, unresolved scales represent 
+a major source of uncertainties in the computation of the large scale horizontal density gradient 
+(from T/S large scale fields), and that the impact of these uncertainties can be simulated 
+by random processes representing unresolved T/S fluctuations.
+
+The stochastic formulation of the equation of state can be written as:
+\begin{equation}
+ \label{eq:eos_sto}
+  \rho = \frac{1}{2} \sum_{i=1}^m\{ \rho[T+\Delta T_i,S+\Delta S_i,p_o(z)] + \rho[T-\Delta T_i,S-\Delta S_i,p_o(z)] \}
+\end{equation}
+where $p_o(z)$ is the reference pressure depending on the depth and, 
+$\Delta T_i$ and $\Delta S_i$ are a set of T/S perturbations defined as the scalar product 
+of the respective local T/S gradients with random walks $\mathbf{\xi}$:
+\begin{equation}
+ \label{eq:sto_pert}
+ \Delta T_i = \mathbf{\xi}_i \cdot \nabla T \qquad \hbox{and} \qquad \Delta S_i = \mathbf{\xi}_i \cdot \nabla S
+\end{equation}
+$\mathbf{\xi}_i$ are produced by a first-order autoregressive processes (AR-1) with 
+a parametrized decorrelation time scale, and horizontal and vertical standard deviations $\sigma_s$. 
+$\mathbf{\xi}$ are uncorrelated over the horizontal and fully correlated along the vertical.
+
+
+\section{Stochastic processes}
+\label{STO_the_details}
+
+The starting point of our implementation of stochastic parameterizations
+in NEMO is to observe that many existing parameterizations are based
+on autoregressive processes, which are used as a basic source of randomness
+to transform a deterministic model into a probabilistic model.
+A generic approach is thus to add one single new module in NEMO,
+generating processes with appropriate statistics
+to simulate each kind of uncertainty in the model
+(see \cite{Brankart_al_GMD2015} for more details).
+
+In practice, at every model grid point, independent Gaussian autoregressive
+processes~$\xi^{(i)},\,i=1,\ldots,m$ are first generated
+using the same basic equation:
+
+\begin{equation}
+\label{eq:autoreg}
+\xi^{(i)}_{k+1} = a^{(i)} \xi^{(i)}_k + b^{(i)} w^{(i)} + c^{(i)}
+\end{equation}
+
+\noindent
+where $k$ is the index of the model timestep; and
+$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are parameters defining
+the mean ($\mu^{(i)}$) standard deviation ($\sigma^{(i)}$)
+and correlation timescale ($\tau^{(i)}$) of each process:
+
+\begin{itemize}
+\item for order~1 processes, $w^{(i)}$ is a Gaussian white noise,
+with zero mean and standard deviation equal to~1, and the parameters
+$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by:
+
+\begin{equation}
+\label{eq:ord1}
+\left\{
+\begin{array}{l}
+a^{(i)} = \varphi \\
+b^{(i)} = \sigma^{(i)} \sqrt{ 1 - \varphi^2 } 
+ \qquad\qquad\mbox{with}\qquad\qquad
+\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\
+c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\
+\end{array}
+\right.
+\end{equation}
+
+\item for order~$n>1$ processes, $w^{(i)}$ is an order~$n-1$ autoregressive process,
+with zero mean, standard deviation equal to~$\sigma^{(i)}$; correlation timescale
+equal to~$\tau^{(i)}$; and the parameters $a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by:
+
+\begin{equation}
+\label{eq:ord2}
+\left\{
+\begin{array}{l}
+a^{(i)} = \varphi \\
+b^{(i)} = \frac{n-1}{2(4n-3)} \sqrt{ 1 - \varphi^2 } 
+ \qquad\qquad\mbox{with}\qquad\qquad
+\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\
+c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\
+\end{array}
+\right.
+\end{equation}
+
+\end{itemize}
+
+\noindent
+In this way, higher order processes can be easily generated recursively using 
+the same piece of code implementing Eq.~(\ref{eq:autoreg}), 
+and using succesively processes from order $0$ to~$n-1$ as~$w^{(i)}$.
+The parameters in Eq.~(\ref{eq:ord2}) are computed so that this recursive application
+of Eq.~(\ref{eq:autoreg}) leads to processes with the required standard deviation
+and correlation timescale, with the additional condition that
+the $n-1$ first derivatives of the autocorrelation function
+are equal to zero at~$t=0$, so that the resulting processes
+become smoother and smoother as $n$ is increased.
+
+Overall, this method provides quite a simple and generic way of generating 
+a wide class of stochastic processes. 
+However, this also means that new model parameters are needed to specify each of 
+these stochastic processes. As in any parameterization of lacking physics, 
+a very important issues then to tune these new parameters using either first principles, 
+model simulations, or real-world observations.
+
+\section{Implementation details}
+\label{STO_thech_details}
+
 %---------------------------------------namsbc--------------------------------------------------
 \namdisplay{namsto}
 %--------------------------------------------------------------------------------------------------------------
-$\ $\newline    % force a new ligne
 
+The computer code implementing stochastic parametrisations can be found in the STO directory.
+It involves three modules : 
+\begin{description}
+\item[\mdl{stopar}] : define the Stochastic parameters and their time evolution.
+\item[\mdl{storng}] : a random number generator based on (and includes) the 64-bit KISS 
+                      (Keep It Simple Stupid) random number generator distributed by George Marsaglia 
+                      (see \href{https://groups.google.com/forum/#!searchin/comp.lang.fortran/64-bit$20KISS$20RNGs}{here})
+\item[\mdl{stopts}] : stochastic parametrisation associated with the non-linearity of the equation of seawater, 
+ implementing Eq~\ref{eq:sto_pert} and specific piece of code in the equation of state implementing Eq~\ref{eq:eos_sto}.
+\end{description}
 
-See \cite{Brankart_OM2013} and \cite{Brankart_al_GMD2015} papers for a description of the parameterization.
+The \mdl{stopar} module has 3 public routines to be called by the model (in our case, NEMO):
+
+The first routine (\rou{sto\_par}) is a direct implementation of Eq.~(\ref{eq:autoreg}),
+applied at each model grid point (in 2D or 3D), 
+and called at each model time step ($k$) to update
+every autoregressive process ($i=1,\ldots,m$).
+This routine also includes a filtering operator, applied to $w^{(i)}$,
+to introduce a spatial correlation between the stochastic processes.
+
+The second routine (\rou{sto\_par\_init}) is an initialization routine mainly dedicated
+to the computation of parameters $a^{(i)}, b^{(i)}, c^{(i)}$
+for each autoregressive process, as a function of the statistical properties
+required by the model user (mean, standard deviation, time correlation,
+order of the process,\ldots). 
+
+Parameters for the processes can be specified through the following \ngn{namsto} namelist parameters:
+\begin{description}
+   \item[\np{nn\_sto\_eos}]   : number of independent random walks 
+   \item[\np{rn\_eos\_stdxy}] : random walk horz. standard deviation (in grid points)
+   \item[\np{rn\_eos\_stdz}]  : random walk vert. standard deviation (in grid points)
+   \item[\np{rn\_eos\_tcor}]  : random walk time correlation (in timesteps)
+   \item[\np{nn\_eos\_ord}]   : order of autoregressive processes
+   \item[\np{nn\_eos\_flt}]   : passes of Laplacian filter
+   \item[\np{rn\_eos\_lim}]   : limitation factor (default = 3.0)
+\end{description}
+This routine also includes the initialization (seeding) of the random number generator.
+
+The third routine (\rou{sto\_rst\_write}) writes a restart file (which suffix name is 
+given by \np{cn\_storst\_out} namelist parameter) containing the current value of 
+all autoregressive processes to allow restarting a simulation from where it has been interrupted.
+This file also contains the current state of the random number generator.
+When \np{ln\_rststo} is set to \textit{true}), the restart file (which suffix name is 
+given by \np{cn\_storst\_in} namelist parameter) is read by the initialization routine 
+(\rou{sto\_par\_init}). The simulation will continue exactly as if it was not interrupted
+only  when \np{ln\_rstseed} is set to \textit{true}, $i.e.$ when the state of 
+the random number generator is read in the restart file.
Index: trunk/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 6409)
+++ trunk/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 6497)
@@ -734,6 +734,6 @@
 (see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
 
-$\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, (see \S\ref{SBC_isf} for further details 
-on how the ice shelf melt is computed and applied).
+$\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, 
+(see \S\ref{SBC_isf} for further details on how the ice shelf melt is computed and applied).
 
 The surface boundary condition on temperature and salinity is applied as follows:
@@ -840,11 +840,19 @@
 ($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform 
 chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb} 
-in \mdl{trc\_oce} module). Three types of chlorophyll can be chosen in the RGB formulation:
-(1) a constant 0.05 g.Chl/L value everywhere (\np{nn\_chdta}=0) ; (2) an observed 
-time varying chlorophyll (\np{nn\_chdta}=1) ; (3) simulated time varying chlorophyll
-by TOP biogeochemical model (\np{ln\_qsr\_bio}=true). In the latter case, the RGB 
-formulation is used to calculate both the phytoplankton light limitation in PISCES 
-or LOBSTER and the oceanic heating rate. 
-
+in \mdl{trc\_oce} module). Four types of chlorophyll can be chosen in the RGB formulation:
+\begin{description} 
+\item[\np{nn\_chdta}=0] 
+a constant 0.05 g.Chl/L value everywhere ; 
+\item[\np{nn\_chdta}=1]  
+an observed time varying chlorophyll deduced from satellite surface ocean color measurement 
+spread uniformly in the vertical direction ; 
+\item[\np{nn\_chdta}=2]  
+same as previous case except that a vertical profile of chlorophyl is used. 
+Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value ;
+\item[\np{ln\_qsr\_bio}=true]  
+simulated time varying chlorophyll by TOP biogeochemical model. 
+In this case, the RGB formulation is used to calculate both the phytoplankton 
+light limitation in PISCES or LOBSTER and the oceanic heating rate. 
+\end{description} 
 The trend in \eqref{Eq_tra_qsr} associated with the penetration of the solar radiation 
 is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. 
@@ -1385,10 +1393,10 @@
                    I've changed "derivative" to "difference" and "mean" to "average"}
 
-With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general, tracers in horizontally 
-adjacent cells live at different depths. Horizontal gradients of tracers are needed 
-for horizontal diffusion (\mdl{traldf} module) and for the hydrostatic pressure 
-gradient (\mdl{dynhpg} module) to be active. The partial cell properties 
-at the top (\np{ln\_isfcav}=true) are computed in the same way as for the bottom. So, only the bottom interpolation is shown.
-\gmcomment{STEVEN from gm : question: not sure of  what -to be active- means}
+With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general, 
+tracers in horizontally adjacent cells live at different depths. 
+Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) 
+and the hydrostatic pressure gradient calculations (\mdl{dynhpg} module). 
+The partial cell properties at the top (\np{ln\_isfcav}=true) are computed in the same way as for the bottom. 
+So, only the bottom interpolation is explained below.
 
 Before taking horizontal gradients between the tracers next to the bottom, a linear 
Index: trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 6409)
+++ trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 6497)
@@ -262,9 +262,9 @@
 \end{equation}
 
-At the ocean surface, a non zero length scale is set through the  \np{rn\_lmin0} namelist 
+At the ocean surface, a non zero length scale is set through the  \np{rn\_mxl0} namelist 
 parameter. Usually the surface scale is given by $l_o = \kappa \,z_o$ 
 where $\kappa = 0.4$ is von Karman's constant and $z_o$ the roughness 
 parameter of the surface. Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94} 
-leads to a 0.04~m, the default value of \np{rn\_lsurf}. In the ocean interior 
+leads to a 0.04~m, the default value of \np{rn\_mxl0}. In the ocean interior 
 a minimum length scale is set to recover the molecular viscosity when $\bar{e}$ 
 reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
@@ -295,6 +295,6 @@
 As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, 
 with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}~=~67.83 corresponds 
-to $\alpha_{CB} = 100$. further setting  \np{ln\_lsurf} to true applies \eqref{ZDF_Lsbc} 
-as surface boundary condition on length scale, with $\beta$ hard coded to the Stacet's value.
+to $\alpha_{CB} = 100$. Further setting  \np{ln\_mxl0} to true applies \eqref{ZDF_Lsbc} 
+as surface boundary condition on length scale, with $\beta$ hard coded to the Stacey's value.
 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) 
 is applied on surface $\bar{e}$ value.
@@ -852,5 +852,6 @@
 The bottom friction represents the friction generated by the bathymetry. 
 The top friction represents the friction generated by the ice shelf/ocean interface. 
-As the friction processes at the top and bottom are represented similarly, only the bottom friction is described in detail below.\\
+As the friction processes at the top and bottom are treated in similar way, 
+only the bottom friction is described in detail below.
 
 
@@ -926,5 +927,5 @@
 $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{rn\_bfric1} (namelist parameter).
+of 115~days. It can be changed by specifying \np{rn\_bfri1} (namelist parameter).
 
 For the linear friction case the coefficients defined in the general 
@@ -936,5 +937,5 @@
 \end{split}
 \end{equation}
-When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfric1}. 
+When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfri1}. 
 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}. 
@@ -943,5 +944,5 @@
 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}.
+by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -963,5 +964,5 @@
 $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} 
+The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} 
 namelist parameters).
 
@@ -978,10 +979,41 @@
 \end{equation}
 
-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.
+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).  This works in the same way
+as for the linear bottom friction case with non-zero masked locations increased by
+$mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}.
+
+% -------------------------------------------------------------------------------------------------------------
+%       Bottom Friction Log-layer
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Log-layer Bottom Friction enhancement (\np{nn\_botfr} = 2, \np{ln\_loglayer} = .true.)}
+\label{ZDF_bfr_loglayer}
+
+In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally
+enhanced using a "law of the wall" scaling. If  \np{ln\_loglayer} = .true., $C_D$ is no
+longer constant but is related to the thickness of the last wet layer in each column by:
+
+\begin{equation}
+C_D = \left ( {\kappa \over {\rm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2
+\end{equation}
+
+\noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness
+length provided via the namelist.
+
+For stability, the drag coefficient is bounded such that it is kept greater or equal to
+the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional
+namelist parameter: \np{rn\_bfri2\_max}, i.e.:
+
+\begin{equation}
+rn\_bfri2 \leq C_D \leq rn\_bfri2\_max
+\end{equation}
+
+\noindent Note also that a log-layer enhancement can also be applied to the top boundary
+friction if under ice-shelf cavities are in use (\np{ln\_isfcav}=.true.).  In this case, the
+relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2}
+and \np{rn\_tfri2\_max}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1267,2 +1299,58 @@
 
 % ================================================================
+% Internal wave-driven mixing
+% ================================================================
+\section{Internal wave-driven mixing (\key{zdftmx\_new})}
+\label{ZDF_tmx_new}
+
+%--------------------------------------------namzdf_tmx_new------------------------------------------
+\namdisplay{namzdf_tmx_new}
+%--------------------------------------------------------------------------------------------------------------
+
+The parameterization of mixing induced by breaking internal waves is a generalization 
+of the approach originally proposed by \citet{St_Laurent_al_GRL02}. 
+A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed, 
+and the resulting diffusivity is obtained as 
+\begin{equation} \label{Eq_Kwave}
+A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 }
+\end{equation}
+where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution 
+of the energy available for mixing. If the \np{ln\_mevar} namelist parameter is set to false, 
+the mixing efficiency is taken as constant and equal to 1/6 \citep{Osborn_JPO80}. 
+In the opposite (recommended) case, $R_f$ is instead a function of the turbulence intensity parameter 
+$Re_b = \frac{ \epsilon}{\nu \, N^2}$, with $\nu$ the molecular viscosity of seawater, 
+following the model of \cite{Bouffard_Boegman_DAO2013} 
+and the implementation of \cite{de_lavergne_JPO2016_efficiency}.
+Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when the mixing efficiency is constant.
+
+In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary 
+as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice). 
+This parameterization of differential mixing, due to \cite{Jackson_Rehmann_JPO2014}, 
+is implemented as in \cite{de_lavergne_JPO2016_efficiency}.
+
+The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$, is constructed 
+from three static maps of column-integrated internal wave energy dissipation, $E_{cri}(i,j)$, 
+$E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures 
+(de Lavergne et al., in prep):
+\begin{align*}
+F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\
+F_{pyc}(i,j,k) &\propto N^{n\_p}\\
+F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} }
+\end{align*} 
+In the above formula, $h_{ab}$ denotes the height above bottom, 
+$h_{wkb}$ denotes the WKB-stretched height above bottom, defined by
+\begin{equation*}
+h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; ,
+\end{equation*}
+The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist)  controls the stratification-dependence of the pycnocline-intensified dissipation. 
+It can take values of 1 (recommended) or 2.
+Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of 
+the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps. 
+$h_{cri}$ is related to the large-scale topography of the ocean (etopo2) 
+and $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of 
+the abyssal hill topography \citep{Goff_JGR2010} and the latitude.
+
+% ================================================================
+
+
+
Index: trunk/DOC/TexFiles/Namelist/namzdf_tmx_new
===================================================================
--- trunk/DOC/TexFiles/Namelist/namzdf_tmx_new	(revision 6497)
+++ trunk/DOC/TexFiles/Namelist/namzdf_tmx_new	(revision 6497)
@@ -0,0 +1,7 @@
+!-----------------------------------------------------------------------
+&namzdf_tmx_new !   internal wave-driven mixing parameterization        ("key_zdftmx_new" & "key_zdfddm")
+!-----------------------------------------------------------------------
+   nn_zpyc     = 1         !  pycnocline-intensified dissipation scales as N (=1) or N^2 (=2)
+   ln_mevar    = .true.    !  variable (T) or constant (F) mixing efficiency
+   ln_tsdiff   = .true.    !  account for differential T/S mixing (T) or not (F)
+/
