Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/NEMO_book.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/NEMO_book.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/NEMO_book.tex	(revision 6040)
@@ -215,6 +215,6 @@
 %\date{\today}
 \date{
-November 2014  \\
-{\small  -- version 3.6 --} \\
+November 2015  \\
+{\small  -- draft of version 4.0 --} \\
 ~  \\
 \textit{\small Note du P\^ole de mod\'{e}lisation de l'Institut Pierre-Simon Laplace No 27 }\\
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Biblio/Biblio.bib	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Biblio/Biblio.bib	(revision 6040)
@@ -276,7 +276,7 @@
   author = {A. Beckmann and H. Goosse},
   title = {A parameterization of ice shelf-ocean interaction for climate models},
-  journal = OM
-  year = {2003}
-  volume = {5}
+  journal = OM,
+  year = {2003},
+  volume = {5},
   pages = {157--170}
 }
@@ -1124,4 +1124,14 @@
 }
 
+@ARTICLE{Graham_McDougall_JPO13,
+  author = {F.S. Graham and T.J. McDougall},
+  title = {Quantifying the nonconservative production of conservative temperature, potential temperature, and entropy},
+  journal = JPO,
+  year = {2013},
+  volume = {43},
+  pages = {838--862},
+  url ={http://dx.doi.org/10.1175/JPO-D-11-0188.1}
+}
+
 @ARTICLE{Greatbatch_JGR94,
   author = {R. J. Greatbatch},
@@ -1428,4 +1438,14 @@
   publisher = {LA-CC-06-012, Los Alamos National Laboratory, N.M.},
   year = {2008}
+}
+
+@TECHREPORT{TEOS10,
+  author = {IOC and SCOR and IAPSO},
+  title = {The international thermodynamic equation of seawater - 2010: Calculation and use of thermodynamic properties},
+  institution = {Intergovernmental Oceanographic Commission},
+  publisher = {Manuals and Guides No. 56, UNESCO (English)},
+  year = {2010},
+  pages = {196pp},
+  url = {http://www.teos-10.org/pubs/TEOS-10_Manual.pdf}
 }
 
@@ -2488,4 +2508,27 @@
 }
 
+@article{Roquet_OM2015,
+author = "F. Roquet and G. Madec and T.J. McDougall and P.M. Barker",
+title = "Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard ",
+journal = OM,
+volume = "90",
+pages = "29--43",
+year = "2015",
+issn = "1463-5003",
+doi = "10.1016/j.ocemod.2015.04.002",
+url = "http://dx.doi.org/10.1016/j.ocemod.2015.04.002"
+}
+
+@article{Roquet_JPO2015,
+author = "F. Roquet and G. Madec and L. Brodeau and J. Nycander",
+title = "Defining a Simplified Yet Realistic Equation of State for Seawater",
+journal = JPO,
+volume = "45",
+pages = "2564--2579",
+year = "2015",
+doi = "10.1175/JPO-D-15-0080.1",
+url = "http://dx.doi.org/10.1175/JPO-D-15-0080.1"
+}
+
 @ARTICLE{Roullet_Madec_JGR00,
   author = {G. Roullet and G. Madec},
@@ -2842,4 +2885,13 @@
   volume = {27},
   pages = {54--69}
+}
+
+@book{Vallis06,
+author = {Vallis, G. K.},
+title = {Atmospheric and Oceanic Fluid Dynamics},
+publisher = {Cambridge University Press},
+address = {Cambridge, U.K.},
+year = {2006},
+pages = {745},
 }
 
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Abstracts_Foreword.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Abstracts_Foreword.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Abstracts_Foreword.tex	(revision 6040)
@@ -13,13 +13,13 @@
 be a flexible tool for studying the ocean and its interactions with the others components of 
 the earth climate system over a wide range of space and time scales. 
-Prognostic variables are the three-dimensional velocity field, a linear 
-or non-linear sea surface height, the temperature and the salinity. In the horizontal direction, 
-the model uses a curvilinear orthogonal grid and in the vertical direction, a full or partial step 
-$z$-coordinate, or $s$-coordinate, or a mixture of the two. The distribution of variables is a 
-three-dimensional Arakawa C-type grid. Various physical choices are available to describe 
-ocean physics, including TKE, GLS and KPP vertical physics. Within NEMO, the ocean is 
-interfaced with a sea-ice model (LIM v2 and v3), passive tracer and biogeochemical models (TOP) 
-and, via the OASIS coupler, with several atmospheric general circulation models. It also 
-support two-way grid embedding via the AGRIF software.
+Prognostic variables are the three-dimensional velocity field, a non-linear sea surface height, 
+the \textit{Conservative} Temperature and the \textit{Absolut} Salinity. 
+In the horizontal direction, the model uses a curvilinear orthogonal grid and in the vertical direction, 
+a full or partial step $z$-coordinate, or $s$-coordinate, or a mixture of the two. 
+The distribution of variables is a three-dimensional Arakawa C-type grid. 
+Various physical choices are available to describe ocean physics, including TKE, and GLS vertical physics. 
+Within NEMO, the ocean is interfaced with a sea-ice model (LIM or CICE), passive tracer and 
+biogeochemical models (TOP) and, via the OASIS coupler, with several atmospheric general circulation models. 
+It also support two-way grid embedding via the AGRIF software.
 
 % ================================================================
@@ -31,13 +31,13 @@
 interactions avec les autres composantes du syst\`{e}me climatique terrestre. 
 Les variables pronostiques sont le champ tridimensionnel de vitesse, une hauteur de la mer 
-lin\'{e}aire ou non, la temperature et la salinit\'{e}. 
+lin\'{e}aire, la Ttemperature Conservative et la Salinit\'{e} Absolue. 
 La distribution des variables se fait sur une grille C d'Arakawa tridimensionnelle utilisant une 
 coordonn\'{e}e verticale $z$ \`{a} niveaux entiers ou partiels, ou une coordonn\'{e}e s, ou encore 
 une combinaison des deux. Diff\'{e}rents choix sont propos\'{e}s pour d\'{e}crire la physique 
-oc\'{e}anique, incluant notamment des physiques verticales TKE, GLS et KPP. A travers l'infrastructure 
-NEMO, l'oc\'{e}an est interfac\'{e} avec des mod\`{e}les de glace de mer, de biog\'{e}ochimie 
-et de traceurs passifs, et, via le coupleur OASIS, \`{a} plusieurs mod\`{e}les de circulation 
-g\'{e}n\'{e}rale atmosph\'{e}rique. Il supporte \'{e}galement l'embo\^{i}tement interactif de 
-maillages via le logiciel AGRIF.
+oc\'{e}anique, incluant notamment des physiques verticales TKE et GLS. A travers l'infrastructure 
+NEMO, l'oc\'{e}an est interfac\'{e} avec des mod\`{e}les de glace de mer (LIM ou CICE), 
+de biog\'{e}ochimie marine et de traceurs passifs, et, via le coupleur OASIS, \`{a} plusieurs 
+mod\`{e}les de circulation g\'{e}n\'{e}rale atmosph\'{e}rique. 
+Il supporte \'{e}galement l'embo\^{i}tement interactif de maillages via le logiciel AGRIF.
 } 
 
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Annex_C.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Annex_C.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Annex_C.tex	(revision 6040)
@@ -410,5 +410,5 @@
 \end{aligned}   } \right.
 \end{equation} 
-where the indices $i_p$ and $k_p$ take the following value: 
+where the indices $i_p$ and $j_p$ take the following value: 
 $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$,
 and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by: 
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_DOM.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 6040)
@@ -1,4 +1,4 @@
 % ================================================================
-% Chapter 2 � Space and Time Domain (DOM)
+% Chapter 2 ——— Space and Time Domain (DOM)
 % ================================================================
 \chapter{Space Domain (DOM) }
@@ -364,5 +364,5 @@
 For both grids here,  the same $w$-point depth has been chosen but in (a) the 
 $t$-points are set half way between $w$-points while in (b) they are defined from 
-an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$. 
+an analytical function: $z(k)=5\,(k-1/2)^3 - 45\,(k-1/2)^2 + 140\,(k-1/2) - 150$. 
 Note the resulting difference between the value of the grid-size $\Delta_k$ and 
 those of the scale factor $e_k$. }
@@ -475,11 +475,11 @@
 (d) hybrid $s-z$ coordinate, 
 (e) hybrid $s-z$ coordinate with partial step, and 
-(f) same as (e) but with variable volume associated with the non-linear free surface. 
-Note that the variable volume option (\key{vvl}) can be used with any of the 
+(f) same as (e) but in the non-linear free surface (\np{ln\_linssh}=false). 
+Note that the non-linear free surface can be used with any of the 
 5 coordinates (a) to (e).}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-The choice of a vertical coordinate, even if it is made through a namelist parameter, 
+The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters, 
 must be done once of all at the beginning of an experiment. It is not intended as an 
 option which can be enabled or disabled in the middle of an experiment. Three main 
@@ -488,11 +488,12 @@
 (\np{ln\_zps}~=~true), or generalized, $s$-coordinate (\np{ln\_sco}~=~true). 
 Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate 
-(Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). When using the variable 
-volume option \key{vvl} ($i.e.$ non-linear free surface), the coordinate follow the 
-time-variation of the free surface so that the transformation is time dependent: 
-$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step 
-bathymetry or $s$-coordinate (hybrid and partial step coordinates have not 
-yet been tested in NEMO v2.3). If using $z$-coordinate with partial step bathymetry
-(\np{ln\_zps}~=~true), ocean cavity beneath ice shelves can be open (\np{ln\_isfcav}~=~true).
+(Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). By default a non-linear free surface is used:
+the coordinate follow the time-variation of the free surface so that the transformation is time dependent: 
+$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). When a linear free surface is assumed (\np{ln\_linssh}=true), 
+the vertical coordinate are fixed in time, but the seawater can move up and down across the z=0 surface 
+(in other words, the top of the ocean in not a rigid-lid). 
+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.
 
 Contrary to the horizontal grid, the vertical grid is computed in the code and no 
@@ -519,17 +520,12 @@
 %%%
 
-The arrays describing the grid point depths and vertical scale factors 
-are three dimensional arrays $(i,j,k)$ even in the case of $z$-coordinate with 
-full step bottom topography. In non-linear free surface (\key{vvl}), their knowledge
-is required at \textit{before}, \textit{now} and \textit{after} time step, while they 
-do not vary in time in linear free surface case. 
-To improve the code readability while providing this flexibility, the vertical coordinate 
-and scale factors are defined as functions of 
-$(i,j,k)$ with "fs" as prefix (examples: \textit{fse3t\_b, fse3t\_n, fse3t\_a,} 
-for the  \textit{before}, \textit{now} and \textit{after} scale factors at $t$-point) 
-that can be either three different arrays when \key{vvl} is defined, or a single fixed arrays. 
-These functions are defined in the file \hf{domzgr\_substitute} of the DOM directory. 
-They are used throughout the code, and replaced by the corresponding arrays at 
-the time of pre-processing (CPP capability).
+Unless a linear free surface is used (\np{ln\_linssh}=false), the arrays describing 
+the grid point depths and vertical scale factors are three set of three dimensional arrays $(i,j,k)$ 
+defined at \textit{before}, \textit{now} and \textit{after} time step. The time at which they are
+defined is indicated by a suffix:$\_b$, $\_n$, or $\_a$, respectively. They are updated at each model time step
+using a fixed reference coordinate system which computer names have a $\_0$ suffix. 
+When the linear free surface option is used (\np{ln\_linssh}=true), \textit{before}, \textit{now} 
+and \textit{after} arrays are simply set one for all to their reference counterpart. 
+
 
 % -------------------------------------------------------------------------------------------------------------
@@ -840,6 +836,6 @@
 %        z*- or s*-coordinate
 % -------------------------------------------------------------------------------------------------------------
-\subsection{$z^*$- or $s^*$-coordinate (add \key{vvl}) }
-\label{DOM_zgr_vvl}
+\subsection{$z^*$- or $s^*$-coordinate (\np{ln\_linssh}=false) }
+\label{DOM_zgr_star}
 
 This option is described in the Report by Levier \textit{et al.} (2007), available on 
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 6040)
@@ -1,10 +1,8 @@
 % ================================================================
-% Chapter � Ocean Dynamics (DYN)
+% Chapter ——— Ocean Dynamics (DYN)
 % ================================================================
 \chapter{Ocean Dynamics (DYN)}
 \label{DYN}
 \minitoc
-
-% add a figure for  dynvor ens, ene latices
 
 %\vspace{2.cm}
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_LDF.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_LDF.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_LDF.tex	(revision 6040)
@@ -1,5 +1,5 @@
 
 % ================================================================
-% Chapter � Lateral Ocean Physics (LDF)
+% Chapter ———  Lateral Ocean Physics (LDF)
 % ================================================================
 \chapter{Lateral Ocean Physics (LDF)}
@@ -15,12 +15,11 @@
 described in \S\ref{PE_zdf} and their discrete formulation in \S\ref{TRA_ldf} 
 and \S\ref{DYN_ldf}). In this section we further discuss each lateral physics option. 
-Choosing one lateral physics scheme means for the user defining, (1) the space 
-and time variations of the eddy coefficients ; (2) the direction along which the 
-lateral diffusive fluxes are evaluated (model level, geopotential or isopycnal 
-surfaces); and (3) the type of operator used (harmonic, or biharmonic operators, 
-and for tracers only, eddy induced advection on tracers). These three aspects 
-of the lateral diffusion are set through namelist parameters and CPP keys 
-(see the \textit{\ngn{nam\_traldf}} and \textit{\ngn{nam\_dynldf}} below). Note
-that this chapter describes the default implementation of iso-neutral
+Choosing one lateral physics scheme means for the user defining, 
+(1) the type of operator used (laplacian or bilaplacian operators, or no lateral mixing term) ; 
+(2) the direction along which the lateral diffusive fluxes are evaluated (model level, geopotential or isopycnal surfaces) ; and 
+(3) the space and time variations of the eddy coefficients. 
+These three aspects of the lateral diffusion are set through namelist parameters 
+(see the \textit{\ngn{nam\_traldf}} and \textit{\ngn{nam\_dynldf}} below). 
+Note that this chapter describes the standard implementation of iso-neutral
 tracer mixing, and Griffies's implementation, which is used if
 \np{traldf\_grif}=true, is described in Appdx\ref{sec:triad}
@@ -33,4 +32,283 @@
 
 % ================================================================
+% Direction of lateral Mixing
+% ================================================================
+\section  [Direction of Lateral Mixing (\textit{ldfslp})]
+		{Direction of Lateral Mixing (\mdl{ldfslp})}
+\label{LDF_slp}
+
+%%%
+\gmcomment{  we should emphasize here that the implementation is a rather old one. 
+Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme. }
+
+A direction for lateral mixing has to be defined when the desired operator does 
+not act along the model levels. This occurs when $(a)$ horizontal mixing is 
+required on tracer or momentum (\np{ln\_traldf\_hor} or \np{ln\_dynldf\_hor}) 
+in $s$- or mixed $s$-$z$- coordinates, and $(b)$ isoneutral mixing is required 
+whatever the vertical coordinate is. This direction of mixing is defined by its 
+slopes in the \textbf{i}- and \textbf{j}-directions at the face of the cell of the 
+quantity to be diffused. For a tracer, this leads to the following four slopes : 
+$r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \eqref{Eq_tra_ldf_iso}), while 
+for momentum the slopes are  $r_{1t}$, $r_{1uw}$, $r_{2f}$, $r_{2uw}$ for 
+$u$ and  $r_{1f}$, $r_{1vw}$, $r_{2t}$, $r_{2vw}$ for $v$. 
+
+%gm% add here afigure of the slope in i-direction
+
+\subsection{slopes for tracer geopotential mixing in the $s$-coordinate}
+
+In $s$-coordinates, geopotential mixing ($i.e.$ horizontal mixing) $r_1$ and 
+$r_2$ are the slopes between the geopotential and computational surfaces. 
+Their discrete formulation is found by locally solving \eqref{Eq_tra_ldf_iso} 
+when the diffusive fluxes in the three directions are set to zero and $T$ is 
+assumed to be horizontally uniform, $i.e.$ a linear function of $z_T$, the 
+depth of a $T$-point. 
+%gm { Steven : My version is obviously wrong since I'm left with an arbitrary constant which is the local vertical temperature gradient}
+
+\begin{equation} \label{Eq_ldfslp_geo}
+\begin{aligned}
+ r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)}
+ 			  \;\delta_{i+1/2}[z_t] 
+ 		&\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t] 
+\\
+ r_{2v} &= \frac{e_{3v}}{\left( e_{2v}\;\overline{\overline{e_{3w}}}^{\,j+1/2,\,k} \right)} 
+ 			  \;\delta_{j+1/2} [z_t] 
+		&\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t] 
+\\
+ r_{1w} &= \frac{1}{e_{1w}}\;\overline{\overline{\delta_{i+1/2}[z_t]}}^{\,i,\,k+1/2}
+ 		&\approx \frac{1}{e_{1w}}\; \delta_{i+1/2}[z_{uw}] 
+ \\
+ r_{2w} &= \frac{1}{e_{2w}}\;\overline{\overline{\delta_{j+1/2}[z_t]}}^{\,j,\,k+1/2}
+		&\approx \frac{1}{e_{2w}}\; \delta_{j+1/2}[z_{vw}] 
+ \\
+\end{aligned}
+\end{equation}
+
+%gm%  caution I'm not sure the simplification was a good idea! 
+
+These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}=True, 
+and either \np{ln\_traldf\_hor}=True or \np{ln\_dynldf\_hor}=True. 
+
+\subsection{Slopes for tracer iso-neutral mixing}\label{LDF_slp_iso}
+In iso-neutral mixing  $r_1$ and $r_2$ are the slopes between the iso-neutral 
+and computational surfaces. Their formulation does not depend on the vertical 
+coordinate used. Their discrete formulation is found using the fact that the 
+diffusive fluxes of locally referenced potential density ($i.e.$ $in situ$ density) 
+vanish. So, substituting $T$ by $\rho$ in \eqref{Eq_tra_ldf_iso} and setting the 
+diffusive fluxes in the three directions to zero leads to the following definition for 
+the neutral slopes:
+
+\begin{equation} \label{Eq_ldfslp_iso}
+\begin{split}
+ r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac{\delta_{i+1/2}[\rho]}
+ 								{\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,i+1/2,\,k}}
+\\
+ r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac{\delta_{j+1/2}\left[\rho \right]}
+ 								{\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,j+1/2,\,k}}
+\\
+ r_{1w} &= \frac{e_{3w}}{e_{1w}}\; 
+ 			\frac{\overline{\overline{\delta_{i+1/2}[\rho]}}^{\,i,\,k+1/2}}
+				 {\delta_{k+1/2}[\rho]}
+\\
+ r_{2w} &= \frac{e_{3w}}{e_{2w}}\; 
+ 			\frac{\overline{\overline{\delta_{j+1/2}[\rho]}}^{\,j,\,k+1/2}}
+				 {\delta_{k+1/2}[\rho]}
+\\
+\end{split}
+\end{equation}
+
+%gm% rewrite this as the explanation is not very clear !!!
+%In practice, \eqref{Eq_ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \eqref{Eq_ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth. 
+
+%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{McDougall1987}, therefore in \eqref{Eq_ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
+
+%In the $z$-coordinate, the derivative of the  \eqref{Eq_ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation. 
+
+As the mixing is performed along neutral surfaces, the gradient of $\rho$ in 
+\eqref{Eq_ldfslp_iso} has to be evaluated at the same local pressure (which, 
+in decibars, is approximated by the depth in meters in the model). Therefore 
+\eqref{Eq_ldfslp_iso} cannot be used as such, but further transformation is 
+needed depending on the vertical coordinate used:
+
+\begin{description}
+
+\item[$z$-coordinate with full step : ] in \eqref{Eq_ldfslp_iso} the densities 
+appearing in the $i$ and $j$ derivatives  are taken at the same depth, thus 
+the $in situ$ density can be used. This is not the case for the vertical 
+derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$, where $N^2$ 
+is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following 
+\citet{McDougall1987} (see \S\ref{TRA_bn2}). 
+
+\item[$z$-coordinate with partial step : ] this case is identical to the full step 
+case except that at partial step level, the \emph{horizontal} density gradient 
+is evaluated as described in \S\ref{TRA_zpshde}.
+
+\item[$s$- or hybrid $s$-$z$- coordinate : ] in the current release of \NEMO, 
+iso-neutral mixing is only employed for $s$-coordinates if the
+Griffies scheme is used (\np{traldf\_grif}=true; see Appdx \ref{sec:triad}). 
+In other words, iso-neutral mixing will only be accurately represented with a 
+linear equation of state (\np{nn\_eos}=1 or 2). In the case of a "true" equation 
+of state, the evaluation of $i$ and $j$ derivatives in \eqref{Eq_ldfslp_iso} 
+will include a pressure dependent part, leading to the wrong evaluation of 
+the neutral slopes.
+
+%gm% 
+Note: The solution for $s$-coordinate passes trough the use of different 
+(and better) expression for the constraint on iso-neutral fluxes. Following 
+\citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral 
+diffusive flux of locally referenced potential density, we stay in the $T$-$S$ 
+plane and consider the balance between the neutral direction diffusive fluxes 
+of potential temperature and salinity:
+\begin{equation}
+\alpha \ \textbf{F}(T) = \beta \ \textbf{F}(S)
+\end{equation}
+%gm{  where vector F is ....}
+
+This constraint leads to the following definition for the slopes:
+
+\begin{equation} \label{Eq_ldfslp_iso2}
+\begin{split}
+ r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac
+ 		{\alpha_u \;\delta_{i+1/2}[T] - \beta_u \;\delta_{i+1/2}[S]}
+		{\alpha_u \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,i+1/2,\,k}
+		 -\beta_u  \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,i+1/2,\,k} }
+\\
+ r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac
+ 		{\alpha_v \;\delta_{j+1/2}[T] - \beta_v \;\delta_{j+1/2}[S]}
+		{\alpha_v \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,j+1/2,\,k}
+		 -\beta_v  \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,j+1/2,\,k} }
+\\
+ r_{1w} &= \frac{e_{3w}}{e_{1w}}\; \frac
+		{\alpha_w \;\overline{\overline{\delta_{i+1/2}[T]}}^{\,i,\,k+1/2}
+		 -\beta_w  \;\overline{\overline{\delta_{i+1/2}[S]}}^{\,i,\,k+1/2} }
+ 		{\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]}
+\\
+ r_{2w} &= \frac{e_{3w}}{e_{2w}}\; \frac
+		{\alpha_w \;\overline{\overline{\delta_{j+1/2}[T]}}^{\,j,\,k+1/2}
+		 -\beta_w  \;\overline{\overline{\delta_{j+1/2}[S]}}^{\,j,\,k+1/2} }
+ 		{\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]}
+\\
+\end{split}
+\end{equation}
+where $\alpha$ and $\beta$, the thermal expansion and saline contraction 
+coefficients introduced in \S\ref{TRA_bn2}, have to be evaluated at the three 
+velocity points. In order to save computation time, they should be approximated 
+by the mean of their values at $T$-points (for example in the case of $\alpha$:  
+$\alpha_u=\overline{\alpha_T}^{i+1/2}$,  $\alpha_v=\overline{\alpha_T}^{j+1/2}$ 
+and $\alpha_w=\overline{\alpha_T}^{k+1/2}$).
+
+Note that such a formulation could be also used in the $z$-coordinate and 
+$z$-coordinate with partial steps cases.
+
+\end{description}
+
+This implementation is a rather old one. It is similar to the one
+proposed by Cox [1987], except for the background horizontal
+diffusion. Indeed, the Cox implementation of isopycnal diffusion in
+GFDL-type models requires a minimum background horizontal diffusion
+for numerical stability reasons.  To overcome this problem, several
+techniques have been proposed in which the numerical schemes of the
+ocean model are modified \citep{Weaver_Eby_JPO97,
+  Griffies_al_JPO98}. Griffies's scheme is now available in \NEMO if
+\np{traldf\_grif\_iso} is set true; see Appdx \ref{sec:triad}. Here,
+another strategy is presented \citep{Lazar_PhD97}: a local
+filtering of the iso-neutral slopes (made on 9 grid-points) prevents
+the development of grid point noise generated by the iso-neutral
+diffusion operator (Fig.~\ref{Fig_LDF_ZDF1}). This allows an
+iso-neutral diffusion scheme without additional background horizontal
+mixing. This technique can be viewed as a diffusion operator that acts
+along large-scale (2~$\Delta$x) \gmcomment{2deltax doesnt seem very
+  large scale} iso-neutral surfaces. The diapycnal diffusion required
+for numerical stability is thus minimized and its net effect on the
+flow is quite small when compared to the effect of an horizontal
+background mixing.
+
+Nevertheless, this iso-neutral operator does not ensure that variance cannot increase, 
+contrary to the \citet{Griffies_al_JPO98} operator which has that property. 
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!ht]      \begin{center}
+\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_LDF_ZDF1.pdf}
+\caption {    \label{Fig_LDF_ZDF1}
+averaging procedure for isopycnal slope computation.}
+\end{center}    \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+%There are three additional questions about the slope calculation. 
+%First the expression for the rotation tensor has been obtain assuming the "small slope" approximation, so a bound has to be imposed on slopes. 
+%Second, numerical stability issues also require a bound on slopes. 
+%Third, the question of boundary condition specified on slopes...
+
+%from griffies: chapter 13.1....
+
+
+
+% In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04}, 
+% the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly 
+% to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the 
+% surface motivates this flattening of isopycnals near the surface).
+
+For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also 
+be bounded by $1/100$ everywhere. This constraint is applied in a piecewise linear 
+fashion, increasing from zero at the surface to $1/100$ at $70$ metres and thereafter 
+decreasing to zero at the bottom of the ocean. (the fact that the eddies "feel" the 
+surface motivates this flattening of isopycnals near the surface).
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!ht]     \begin{center}
+\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_eiv_slp.pdf}
+\caption {     \label{Fig_eiv_slp}
+Vertical profile of the slope used for lateral mixing in the mixed layer : 
+\textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior, 
+which has to be adjusted at the surface boundary (i.e. it must tend to zero at the 
+surface since there is no mixing across the air-sea interface: wall boundary 
+condition). Nevertheless, the profile between the surface zero value and the interior 
+iso-neutral one is unknown, and especially the value at the base of the mixed layer ; 
+\textit{(b)} profile of slope using a linear tapering of the slope near the surface and 
+imposing a maximum slope of 1/100 ; \textit{(c)} profile of slope actually used in 
+\NEMO: a linear decrease of the slope from zero at the surface to its ocean interior 
+value computed just below the mixed layer. Note the huge change in the slope at the 
+base of the mixed layer between  \textit{(b)}  and \textit{(c)}.}
+\end{center}   \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+\colorbox{yellow}{add here a discussion about the flattening of the slopes, vs  tapering the coefficient.}
+
+\subsection{slopes for momentum iso-neutral mixing}
+
+The iso-neutral diffusion operator on momentum is the same as the one used on 
+tracers but applied to each component of the velocity separately (see 
+\eqref{Eq_dyn_ldf_iso} in section~\ref{DYN_ldf_iso}). The slopes between the 
+surface along which the diffusion operator acts and the surface of computation 
+($z$- or $s$-surfaces) are defined at $T$-, $f$-, and \textit{uw}- points for the 
+$u$-component, and $T$-, $f$- and \textit{vw}- points for the $v$-component. 
+They are computed from the slopes used for tracer diffusion, $i.e.$ 
+\eqref{Eq_ldfslp_geo} and \eqref{Eq_ldfslp_iso} :
+
+\begin{equation} \label{Eq_ldfslp_dyn}
+\begin{aligned}
+&r_{1t}\ \ = \overline{r_{1u}}^{\,i}       &&&    r_{1f}\ \ &= \overline{r_{1u}}^{\,i+1/2} \\
+&r_{2f} \ \ = \overline{r_{2v}}^{\,j+1/2} &&& 	r_{2t}\ &= \overline{r_{2v}}^{\,j} \\
+&r_{1uw}  = \overline{r_{1w}}^{\,i+1/2} &&\ \ \text{and} \ \ &   r_{1vw}&= \overline{r_{1w}}^{\,j+1/2} \\
+&r_{2uw}= \overline{r_{2w}}^{\,j+1/2} &&&         r_{2vw}&= \overline{r_{2w}}^{\,j+1/2}\\
+\end{aligned}
+\end{equation}
+
+The major issue remaining is in the specification of the boundary conditions. 
+The same boundary conditions are chosen as those used for lateral 
+diffusion along model level surfaces, i.e. using the shear computed along 
+the model levels and with no additional friction at the ocean bottom (see 
+{\S\ref{LBC_coast}).
+
+
+% ================================================================
+% Lateral Mixing Operator
+% ================================================================
+\section [Lateral Mixing Operators (\textit{ldftra}, \textit{ldfdyn})] 
+		  {Lateral Mixing Operators (\mdl{traldf}, \mdl{traldf}) }
+\label{LDF_op}
+
+
+   
+% ================================================================
 % Lateral Mixing Coefficients
 % ================================================================
@@ -38,5 +316,4 @@
 		  {Lateral Mixing Coefficient (\mdl{ldftra}, \mdl{ldfdyn}) }
 \label{LDF_coef}
-
 
 Introducing a space variation in the lateral eddy mixing coefficients changes 
@@ -165,274 +442,4 @@
 
 % ================================================================
-% Direction of lateral Mixing
-% ================================================================
-\section  [Direction of Lateral Mixing (\textit{ldfslp})]
-		{Direction of Lateral Mixing (\mdl{ldfslp})}
-\label{LDF_slp}
-
-%%%
-\gmcomment{  we should emphasize here that the implementation is a rather old one. 
-Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme. }
-
-A direction for lateral mixing has to be defined when the desired operator does 
-not act along the model levels. This occurs when $(a)$ horizontal mixing is 
-required on tracer or momentum (\np{ln\_traldf\_hor} or \np{ln\_dynldf\_hor}) 
-in $s$- or mixed $s$-$z$- coordinates, and $(b)$ isoneutral mixing is required 
-whatever the vertical coordinate is. This direction of mixing is defined by its 
-slopes in the \textbf{i}- and \textbf{j}-directions at the face of the cell of the 
-quantity to be diffused. For a tracer, this leads to the following four slopes : 
-$r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \eqref{Eq_tra_ldf_iso}), while 
-for momentum the slopes are  $r_{1t}$, $r_{1uw}$, $r_{2f}$, $r_{2uw}$ for 
-$u$ and  $r_{1f}$, $r_{1vw}$, $r_{2t}$, $r_{2vw}$ for $v$. 
-
-%gm% add here afigure of the slope in i-direction
-
-\subsection{slopes for tracer geopotential mixing in the $s$-coordinate}
-
-In $s$-coordinates, geopotential mixing ($i.e.$ horizontal mixing) $r_1$ and 
-$r_2$ are the slopes between the geopotential and computational surfaces. 
-Their discrete formulation is found by locally solving \eqref{Eq_tra_ldf_iso} 
-when the diffusive fluxes in the three directions are set to zero and $T$ is 
-assumed to be horizontally uniform, $i.e.$ a linear function of $z_T$, the 
-depth of a $T$-point. 
-%gm { Steven : My version is obviously wrong since I'm left with an arbitrary constant which is the local vertical temperature gradient}
-
-\begin{equation} \label{Eq_ldfslp_geo}
-\begin{aligned}
- r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)}
- 			  \;\delta_{i+1/2}[z_t] 
- 		&\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t] 
-\\
- r_{2v} &= \frac{e_{3v}}{\left( e_{2v}\;\overline{\overline{e_{3w}}}^{\,j+1/2,\,k} \right)} 
- 			  \;\delta_{j+1/2} [z_t] 
-		&\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t] 
-\\
- r_{1w} &= \frac{1}{e_{1w}}\;\overline{\overline{\delta_{i+1/2}[z_t]}}^{\,i,\,k+1/2}
- 		&\approx \frac{1}{e_{1w}}\; \delta_{i+1/2}[z_{uw}] 
- \\
- r_{2w} &= \frac{1}{e_{2w}}\;\overline{\overline{\delta_{j+1/2}[z_t]}}^{\,j,\,k+1/2}
-		&\approx \frac{1}{e_{2w}}\; \delta_{j+1/2}[z_{vw}] 
- \\
-\end{aligned}
-\end{equation}
-
-%gm%  caution I'm not sure the simplification was a good idea! 
-
-These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}=True, 
-and either \np{ln\_traldf\_hor}=True or \np{ln\_dynldf\_hor}=True. 
-
-\subsection{Slopes for tracer iso-neutral mixing}\label{LDF_slp_iso}
-In iso-neutral mixing  $r_1$ and $r_2$ are the slopes between the iso-neutral 
-and computational surfaces. Their formulation does not depend on the vertical 
-coordinate used. Their discrete formulation is found using the fact that the 
-diffusive fluxes of locally referenced potential density ($i.e.$ $in situ$ density) 
-vanish. So, substituting $T$ by $\rho$ in \eqref{Eq_tra_ldf_iso} and setting the 
-diffusive fluxes in the three directions to zero leads to the following definition for 
-the neutral slopes:
-
-\begin{equation} \label{Eq_ldfslp_iso}
-\begin{split}
- r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac{\delta_{i+1/2}[\rho]}
- 								{\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,i+1/2,\,k}}
-\\
- r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac{\delta_{j+1/2}\left[\rho \right]}
- 								{\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,j+1/2,\,k}}
-\\
- r_{1w} &= \frac{e_{3w}}{e_{1w}}\; 
- 			\frac{\overline{\overline{\delta_{i+1/2}[\rho]}}^{\,i,\,k+1/2}}
-				 {\delta_{k+1/2}[\rho]}
-\\
- r_{2w} &= \frac{e_{3w}}{e_{2w}}\; 
- 			\frac{\overline{\overline{\delta_{j+1/2}[\rho]}}^{\,j,\,k+1/2}}
-				 {\delta_{k+1/2}[\rho]}
-\\
-\end{split}
-\end{equation}
-
-%gm% rewrite this as the explanation is not very clear !!!
-%In practice, \eqref{Eq_ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \eqref{Eq_ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth. 
-
-%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{McDougall1987}, therefore in \eqref{Eq_ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
-
-%In the $z$-coordinate, the derivative of the  \eqref{Eq_ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation. 
-
-As the mixing is performed along neutral surfaces, the gradient of $\rho$ in 
-\eqref{Eq_ldfslp_iso} has to be evaluated at the same local pressure (which, 
-in decibars, is approximated by the depth in meters in the model). Therefore 
-\eqref{Eq_ldfslp_iso} cannot be used as such, but further transformation is 
-needed depending on the vertical coordinate used:
-
-\begin{description}
-
-\item[$z$-coordinate with full step : ] in \eqref{Eq_ldfslp_iso} the densities 
-appearing in the $i$ and $j$ derivatives  are taken at the same depth, thus 
-the $in situ$ density can be used. This is not the case for the vertical 
-derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$, where $N^2$ 
-is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following 
-\citet{McDougall1987} (see \S\ref{TRA_bn2}). 
-
-\item[$z$-coordinate with partial step : ] this case is identical to the full step 
-case except that at partial step level, the \emph{horizontal} density gradient 
-is evaluated as described in \S\ref{TRA_zpshde}.
-
-\item[$s$- or hybrid $s$-$z$- coordinate : ] in the current release of \NEMO, 
-iso-neutral mixing is only employed for $s$-coordinates if the
-Griffies scheme is used (\np{traldf\_grif}=true; see Appdx \ref{sec:triad}). 
-In other words, iso-neutral mixing will only be accurately represented with a 
-linear equation of state (\np{nn\_eos}=1 or 2). In the case of a "true" equation 
-of state, the evaluation of $i$ and $j$ derivatives in \eqref{Eq_ldfslp_iso} 
-will include a pressure dependent part, leading to the wrong evaluation of 
-the neutral slopes.
-
-%gm% 
-Note: The solution for $s$-coordinate passes trough the use of different 
-(and better) expression for the constraint on iso-neutral fluxes. Following 
-\citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral 
-diffusive flux of locally referenced potential density, we stay in the $T$-$S$ 
-plane and consider the balance between the neutral direction diffusive fluxes 
-of potential temperature and salinity:
-\begin{equation}
-\alpha \ \textbf{F}(T) = \beta \ \textbf{F}(S)
-\end{equation}
-%gm{  where vector F is ....}
-
-This constraint leads to the following definition for the slopes:
-
-\begin{equation} \label{Eq_ldfslp_iso2}
-\begin{split}
- r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac
- 		{\alpha_u \;\delta_{i+1/2}[T] - \beta_u \;\delta_{i+1/2}[S]}
-		{\alpha_u \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,i+1/2,\,k}
-		 -\beta_u  \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,i+1/2,\,k} }
-\\
- r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac
- 		{\alpha_v \;\delta_{j+1/2}[T] - \beta_v \;\delta_{j+1/2}[S]}
-		{\alpha_v \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,j+1/2,\,k}
-		 -\beta_v  \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,j+1/2,\,k} }
-\\
- r_{1w} &= \frac{e_{3w}}{e_{1w}}\; \frac
-		{\alpha_w \;\overline{\overline{\delta_{i+1/2}[T]}}^{\,i,\,k+1/2}
-		 -\beta_w  \;\overline{\overline{\delta_{i+1/2}[S]}}^{\,i,\,k+1/2} }
- 		{\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]}
-\\
- r_{2w} &= \frac{e_{3w}}{e_{2w}}\; \frac
-		{\alpha_w \;\overline{\overline{\delta_{j+1/2}[T]}}^{\,j,\,k+1/2}
-		 -\beta_w  \;\overline{\overline{\delta_{j+1/2}[S]}}^{\,j,\,k+1/2} }
- 		{\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]}
-\\
-\end{split}
-\end{equation}
-where $\alpha$ and $\beta$, the thermal expansion and saline contraction 
-coefficients introduced in \S\ref{TRA_bn2}, have to be evaluated at the three 
-velocity points. In order to save computation time, they should be approximated 
-by the mean of their values at $T$-points (for example in the case of $\alpha$:  
-$\alpha_u=\overline{\alpha_T}^{i+1/2}$,  $\alpha_v=\overline{\alpha_T}^{j+1/2}$ 
-and $\alpha_w=\overline{\alpha_T}^{k+1/2}$).
-
-Note that such a formulation could be also used in the $z$-coordinate and 
-$z$-coordinate with partial steps cases.
-
-\end{description}
-
-This implementation is a rather old one. It is similar to the one
-proposed by Cox [1987], except for the background horizontal
-diffusion. Indeed, the Cox implementation of isopycnal diffusion in
-GFDL-type models requires a minimum background horizontal diffusion
-for numerical stability reasons.  To overcome this problem, several
-techniques have been proposed in which the numerical schemes of the
-ocean model are modified \citep{Weaver_Eby_JPO97,
-  Griffies_al_JPO98}. Griffies's scheme is now available in \NEMO if
-\np{traldf\_grif\_iso} is set true; see Appdx \ref{sec:triad}. Here,
-another strategy is presented \citep{Lazar_PhD97}: a local
-filtering of the iso-neutral slopes (made on 9 grid-points) prevents
-the development of grid point noise generated by the iso-neutral
-diffusion operator (Fig.~\ref{Fig_LDF_ZDF1}). This allows an
-iso-neutral diffusion scheme without additional background horizontal
-mixing. This technique can be viewed as a diffusion operator that acts
-along large-scale (2~$\Delta$x) \gmcomment{2deltax doesnt seem very
-  large scale} iso-neutral surfaces. The diapycnal diffusion required
-for numerical stability is thus minimized and its net effect on the
-flow is quite small when compared to the effect of an horizontal
-background mixing.
-
-Nevertheless, this iso-neutral operator does not ensure that variance cannot increase, 
-contrary to the \citet{Griffies_al_JPO98} operator which has that property. 
-
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht]      \begin{center}
-\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_LDF_ZDF1.pdf}
-\caption {    \label{Fig_LDF_ZDF1}
-averaging procedure for isopycnal slope computation.}
-\end{center}    \end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-%There are three additional questions about the slope calculation. 
-%First the expression for the rotation tensor has been obtain assuming the "small slope" approximation, so a bound has to be imposed on slopes. 
-%Second, numerical stability issues also require a bound on slopes. 
-%Third, the question of boundary condition specified on slopes...
-
-%from griffies: chapter 13.1....
-
-
-
-% In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04}, 
-% the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly 
-% to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the 
-% surface motivates this flattening of isopycnals near the surface).
-
-For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also 
-be bounded by $1/100$ everywhere. This constraint is applied in a piecewise linear 
-fashion, increasing from zero at the surface to $1/100$ at $70$ metres and thereafter 
-decreasing to zero at the bottom of the ocean. (the fact that the eddies "feel" the 
-surface motivates this flattening of isopycnals near the surface).
-
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht]     \begin{center}
-\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_eiv_slp.pdf}
-\caption {     \label{Fig_eiv_slp}
-Vertical profile of the slope used for lateral mixing in the mixed layer : 
-\textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior, 
-which has to be adjusted at the surface boundary (i.e. it must tend to zero at the 
-surface since there is no mixing across the air-sea interface: wall boundary 
-condition). Nevertheless, the profile between the surface zero value and the interior 
-iso-neutral one is unknown, and especially the value at the base of the mixed layer ; 
-\textit{(b)} profile of slope using a linear tapering of the slope near the surface and 
-imposing a maximum slope of 1/100 ; \textit{(c)} profile of slope actually used in 
-\NEMO: a linear decrease of the slope from zero at the surface to its ocean interior 
-value computed just below the mixed layer. Note the huge change in the slope at the 
-base of the mixed layer between  \textit{(b)}  and \textit{(c)}.}
-\end{center}   \end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-\colorbox{yellow}{add here a discussion about the flattening of the slopes, vs  tapering the coefficient.}
-
-\subsection{slopes for momentum iso-neutral mixing}
-
-The iso-neutral diffusion operator on momentum is the same as the one used on 
-tracers but applied to each component of the velocity separately (see 
-\eqref{Eq_dyn_ldf_iso} in section~\ref{DYN_ldf_iso}). The slopes between the 
-surface along which the diffusion operator acts and the surface of computation 
-($z$- or $s$-surfaces) are defined at $T$-, $f$-, and \textit{uw}- points for the 
-$u$-component, and $T$-, $f$- and \textit{vw}- points for the $v$-component. 
-They are computed from the slopes used for tracer diffusion, $i.e.$ 
-\eqref{Eq_ldfslp_geo} and \eqref{Eq_ldfslp_iso} :
-
-\begin{equation} \label{Eq_ldfslp_dyn}
-\begin{aligned}
-&r_{1t}\ \ = \overline{r_{1u}}^{\,i}       &&&    r_{1f}\ \ &= \overline{r_{1u}}^{\,i+1/2} \\
-&r_{2f} \ \ = \overline{r_{2v}}^{\,j+1/2} &&& 	r_{2t}\ &= \overline{r_{2v}}^{\,j} \\
-&r_{1uw}  = \overline{r_{1w}}^{\,i+1/2} &&\ \ \text{and} \ \ &   r_{1vw}&= \overline{r_{1w}}^{\,j+1/2} \\
-&r_{2uw}= \overline{r_{2w}}^{\,j+1/2} &&&         r_{2vw}&= \overline{r_{2w}}^{\,j+1/2}\\
-\end{aligned}
-\end{equation}
-
-The major issue remaining is in the specification of the boundary conditions. 
-The same boundary conditions are chosen as those used for lateral 
-diffusion along model level surfaces, i.e. using the shear computed along 
-the model levels and with no additional friction at the ocean bottom (see 
-{\S\ref{LBC_coast}).
-
-
-% ================================================================
 % Eddy Induced Mixing
 % ================================================================
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_Model_Basics.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_Model_Basics.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_Model_Basics.tex	(revision 6040)
@@ -247,75 +247,15 @@
 sufficient to solve a linearized version of (\ref{Eq_PE_ssh}), which still allows 
 to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}.
+Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost.
 
 The filtering of EGWs in models with a free surface is usually a matter of discretisation 
-of the temporal derivatives, using the time splitting method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92} 
-or the implicit scheme \citep{Dukowicz1994}. In \NEMO, we use a slightly different approach 
-developed by \citet{Roullet_Madec_JGR00}: the damping of EGWs is ensured by introducing an 
-additional force in the momentum equation. \eqref{Eq_PE_dyn} becomes: 
-\begin{equation} \label{Eq_PE_flt}
-\frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}}
-- g \nabla \left( \tilde{\rho} \ \eta \right) 
-- g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) 
-\end{equation}
-where $T_c$, is a parameter with dimensions of time which characterizes the force, 
-$\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$ 
-represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 
-non-linear and viscous terms in \eqref{Eq_PE_dyn}.
-
-The new force can be interpreted as a diffusion of vertically integrated volume flux divergence. 
-The time evolution of $D$ is thus governed by a balance of two terms, $-g$ \textbf{A} $\eta$ 
-and $g \, T_c \,$ \textbf{A} $D$, associated with a propagative regime and a diffusive regime 
-in the temporal spectrum, respectively. In the diffusive regime, the EGWs no longer propagate, 
-$i.e.$ they are stationary and damped. The diffusion regime applies to the modes shorter than 
-$T_c$. For longer ones, the diffusion term vanishes. Hence, the temporally unresolved EGWs 
-can be damped by choosing $T_c > \rdt$. \citet{Roullet_Madec_JGR00} demonstrate that 
-(\ref{Eq_PE_flt}) can be integrated with a leap frog scheme except the additional term which 
-has to be computed implicitly. This is not surprising since the use of a large time step has a 
-necessarily numerical cost. Two gains arise in comparison with the previous formulations. 
-Firstly, the damping of EGWs can be quantified through the magnitude of the additional term. 
-Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as 
-soon as $T_c > \rdt$.
-
-When the variations of free surface elevation are small compared to the thickness of the first 
-model layer, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized 
-by \citet{Roullet_Madec_JGR00} the linearization of (\ref{Eq_PE_ssh}) has consequences on the 
-conservation of salt in the model. With the nonlinear free surface equation, the time evolution 
-of the total salt content is 
-\begin{equation} \label{Eq_PE_salt_content}
-    \frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv} 
-                        =\int\limits_S {S\;(-\frac{\partial \eta }{\partial t}-D+P-E)\;ds}
-\end{equation}
-where $S$ is the salinity, and the total salt is integrated over the whole ocean volume 
-$D_\eta$ bounded by the time-dependent free surface. The right hand side (which is an 
-integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh}) 
-is satisfied, so that the salt is perfectly conserved. When the free surface equation is 
-linearized, \citet{Roullet_Madec_JGR00} show that the total salt content integrated in the fixed 
-volume $D$ (bounded by the surface $z=0$) is no longer conserved:
-\begin{equation} \label{Eq_PE_salt_content_linear}
-         \frac{\partial }{\partial t}\int\limits_D {S\;dv} 
-         		= - \int\limits_S {S\;\frac{\partial \eta }{\partial t}ds} 
-\end{equation}
-
-The right hand side of (\ref{Eq_PE_salt_content_linear}) is small in equilibrium solutions 
-\citep{Roullet_Madec_JGR00}. It can be significant when the freshwater forcing is not balanced and 
-the globally averaged free surface is drifting. An increase in sea surface height \textit{$\eta $} 
-results in a decrease of the salinity in the fixed volume $D$. Even in that case though, 
-the total salt integrated in the variable volume $D_{\eta}$ varies much less, since 
-(\ref{Eq_PE_salt_content_linear}) can be rewritten as 
-\begin{equation} \label{Eq_PE_salt_content_corrected}
-\frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv} 
-=\frac{\partial}{\partial t} \left[ \;{\int\limits_D {S\;dv} +\int\limits_S {S\eta \;ds} } \right]
-=\int\limits_S {\eta \;\frac{\partial S}{\partial t}ds}
-\end{equation}
-
-Although the total salt content is not exactly conserved with the linearized free surface, 
-its variations are driven by correlations of the time variation of surface salinity with the 
-sea surface height, which is a negligible term. This situation contrasts with the case of 
-the rigid lid approximation in which case freshwater forcing is represented by a virtual 
-salt flux, leading to a spurious source of salt at the ocean surface 
-\citep{Huang_JPO93, Roullet_Madec_JGR00}.
-
-\newpage
-$\ $\newline    % force a new ligne
+of the temporal derivatives, using a split-explicit method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92} 
+or the implicit scheme \citep{Dukowicz1994} or the addition of a filtering force in the momentum equation 
+\citep{Roullet_Madec_JGR00}. With the present release, \NEMO offers the choice between 
+an explicit free surface (see \S\ref{DYN_spg_exp}) or a split-explicit scheme strongly 
+inspired the one proposed by \citet{Shchepetkin_McWilliams_OM05} (see \S\ref{DYN_spg_ts}).
+
+%\newpage
+%$\ $\newline    % force a new ligne
 
 % ================================================================
@@ -655,9 +595,5 @@
 the surface pressure, is given by:
 \begin{equation} \label{Eq_PE_spg}
-p_s = \left\{ \begin{split} 
-\rho \,g \,\eta &                                 \qquad  \qquad  \;   \qquad \text{ standard free surface} \\ 
-\rho \,g \,\eta &+ \rho_o \,\mu \,\frac{\partial \eta }{\partial t}      \qquad \text{ filtered     free surface}    
-\end{split} 
-\right.
+p_s =  \rho \,g \,\eta 
 \end{equation}
 with $\eta$ is solution of \eqref{Eq_PE_ssh}
@@ -773,10 +709,10 @@
 \end{equation}
 
-The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinate can be written as follows:
+The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinate can be written as follows (see Appendix~\ref{Apdx_A_momentum}):
 
  \vspace{0.5cm}
-* momentum equation:
+$\bullet$ Vector invariant form of the momentum equation :
 \begin{multline} \label{Eq_PE_sco_u}
-\frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t}=
+\frac{\partial  u   }{\partial t}=
 	+   \left( {\zeta +f} \right)\,v                                    
 	-   \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left(  u^2+v^2   \right) 
@@ -787,5 +723,5 @@
 \end{multline}
 \begin{multline} \label{Eq_PE_sco_v}
-\frac{1}{e_3} \frac{\partial \left(  e_3\,v  \right) }{\partial t}=
+\frac{\partial v }{\partial t}=
 	-   \left( {\zeta +f} \right)\,u   
 	-   \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left(  u^2+v^2  \right)     	
@@ -795,9 +731,38 @@
 	+  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad
 \end{multline}
+
+ \vspace{0.5cm}
+$\bullet$ Vector invariant form of the momentum equation :
+\begin{multline} \label{Eq_PE_sco_u}
+\frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t}=
+   +   \left( { f + \frac{1}{e_1 \; e_2 }
+					\left( 	 v \frac{\partial e_2}{\partial i}
+						-u \frac{\partial e_1}{\partial j} 	\right)} 	\right) \, v    \\
+   - \frac{1}{e_1 \; e_2 \; e_3 } 	\left( 
+					\frac{\partial \left( {e_2 \, e_3 \, u\,u} \right)}{\partial i}
+		+			\frac{\partial \left( {e_1 \, e_3 \, v\,u} \right)}{\partial j}	\right)
+   - \frac{1}{e_3 }\frac{\partial \left( { \omega\,u} \right)}{\partial k}    \\
+	- \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho _o}    \right)    
+	+  g\frac{\rho }{\rho _o}\sigma _1 
+	+   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad
+\end{multline}
+\begin{multline} \label{Eq_PE_sco_v}
+\frac{1}{e_3} \frac{\partial \left(  e_3\,v  \right) }{\partial t}=
+   -   \left( { f + \frac{1}{e_1 \; e_2}
+					\left( 	 v \frac{\partial e_2}{\partial i}
+						-u \frac{\partial e_1}{\partial j} 	\right)} 	\right) \, u   \\
+   - \frac{1}{e_1 \; e_2 \; e_3 } 	\left( 
+					\frac{\partial \left( {e_2 \; e_3  \,u\,v} \right)}{\partial i}
+		+			\frac{\partial \left( {e_1 \; e_3  \,v\,v} \right)}{\partial j}	\right)
+                 - \frac{1}{e_3 } \frac{\partial \left( { \omega\,v} \right)}{\partial k}    \\
+	-   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o}  \right) 
+	 +  g\frac{\rho }{\rho _o }\sigma _2   
+	+  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad
+\end{multline}
+
 where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, and the hydrostatic 
 pressure have the same expressions as in $z$-coordinates although they do not represent 
 exactly the same quantities. $\omega$ is provided by the continuity equation 
 (see Appendix~\ref{Apdx_A}):
-
 \begin{equation} \label{Eq_PE_sco_continuity}
 \frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0   
@@ -809,5 +774,5 @@
 
  \vspace{0.5cm}
-* tracer equations:
+$\bullet$ tracer equations:
 \begin{multline} \label{Eq_PE_sco_t}
 \frac{1}{e_3} \frac{\partial \left(  e_3\,T  \right) }{\partial t}=
@@ -1024,5 +989,7 @@
 
 The $\tilde{z}$-coordinate has been developed by \citet{Leclair_Madec_OM10s}.
-It is not available in the current version of \NEMO.
+It is available in \NEMO since the version 3.4. Nevertheless, it is currently not robust enough 
+to be used in all possible configurations. Its use is therefore not recommended.
+We 
 
 \newpage 
@@ -1180,8 +1147,10 @@
 ocean (see Appendix~\ref{Apdx_B}).
 
+For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. $\Re $ reduces to the identity 
+in the horizontal direction, no rotation is applied. 
+
 For \textit{geopotential} diffusion, $r_1$ and $r_2 $ are the slopes between the 
-geopotential and computational surfaces: in $z$-coordinates they are zero 
-($r_1 = r_2 = 0$) while in $s$-coordinate (including $\textit{z*}$ case) they are 
-equal to $\sigma _1$ and $\sigma _2$, respectively (see \eqref{Eq_PE_sco_slope} ).
+geopotential and computational surfaces: they are equal to $\sigma _1$ and $\sigma _2$, 
+respectively (see \eqref{Eq_PE_sco_slope} ).
 
 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral 
@@ -1235,10 +1204,9 @@
 The lateral fourth order tracer diffusive operator is defined by:
 \begin{equation} \label{Eq_PE_bilapT}
-D^{lT}=\Delta \left( {A^{lT}\;\Delta T} \right) 
-\qquad \text{where} \  D^{lT}=\Delta \left( {A^{lT}\;\Delta T} \right)
+D^{lT}=\Delta \left( \;\Delta T \right) 
+\qquad \text{where} \;\; \Delta \bullet = \nabla \left( {\sqrt{B^{lT}\,}\;\Re \;\nabla \bullet} \right)
  \end{equation}
-
 It is the second order operator given by \eqref{Eq_PE_iso_tensor} applied twice with 
-the eddy diffusion coefficient correctly placed. 
+the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 
 
 
@@ -1262,13 +1230,13 @@
 
 Such a formulation ensures a complete separation between the vorticity and 
-horizontal divergence fields (see Appendix~\ref{Apdx_C}). Unfortunately, it is not 
-available for geopotential diffusion in $s-$coordinates and for isoneutral 
-diffusion in both $z$- and $s$-coordinates ($i.e.$ when a rotation is required). 
-In these two cases, the $u$ and $v-$fields are considered as independent scalar 
-fields, so that the diffusive operator is given by:
+horizontal divergence fields (see Appendix~\ref{Apdx_C}). 
+Unfortunately, it is only available in \textit{iso-level} direction. 
+When a rotation is required ($i.e.$ geopotential diffusion in $s-$coordinates 
+or isoneutral diffusion in both $z$- and $s$-coordinates), the $u$ and $v-$fields 
+are considered as independent scalar fields, so that the diffusive operator is given by:
 \begin{equation} \label{Eq_PE_lapU_iso}
 \begin{split}
- D_u^{l{\rm {\bf U}}} &= \nabla .\left( {\Re \;\nabla u} \right) \\ 
- D_v^{l{\rm {\bf U}}} &= \nabla .\left( {\Re \;\nabla v} \right)
+ D_u^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla u} \right) \\ 
+ D_v^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla v} \right)
  \end{split}
  \end{equation}
@@ -1282,29 +1250,5 @@
 
 As for tracers, the fourth order momentum diffusive operator along $z$ or $s$-surfaces 
-is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU} 
-with the eddy viscosity coefficient correctly placed:
-
-geopotential diffusion in $z$-coordinate:
-\begin{equation} \label{Eq_PE_bilapU}
-\begin{split}
-{\rm {\bf D}}^{l{\rm {\bf U}}} &=\nabla _h \left\{ {\;\nabla _h {\rm {\bf 
-.}}\left[ {A^{lm}\,\nabla _h \left( \chi \right)} \right]\;} 
-\right\}\;   \\
-&+\nabla _h \times \left\{ {\;{\rm {\bf k}}\cdot \nabla \times 
-\left[ {A^{lm}\,\nabla _h \times \left( {\zeta \;{\rm {\bf k}}} \right)} 
-\right]\;} \right\}
-\end{split}
-\end{equation}
-
-\gmcomment{  change the position of the coefficient, both here and in the code}
-
-geopotential diffusion in $s$-coordinate:
-\begin{equation} \label{Eq_bilapU_iso}
-   \left\{   \begin{aligned}
-         D_u^{l{\rm {\bf U}}} =\Delta \left( {A^{lm}\;\Delta u} \right) \\ 
-         D_v^{l{\rm {\bf U}}} =\Delta \left( {A^{lm}\;\Delta v} \right)
-   \end{aligned}    \right.
-   \quad \text{where} \quad 
-   \Delta \left( \bullet \right) = \nabla \cdot \left( \Re \nabla(\bullet) \right) 
-\end{equation}
-
+is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU_iso} 
+with the harmonic eddy diffusion coefficient set to the square root of the biharmonic one.
+
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex	(revision 6040)
@@ -1,4 +1,4 @@
 % ================================================================
-% Chapter 1 � Model Basics
+% Chapter 1 ——— Model Basics
 % ================================================================
 % ================================================================
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_STP.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_STP.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_STP.tex	(revision 6040)
@@ -1,5 +1,5 @@
 
 % ================================================================
-% Chapter 2 � Time Domain (step.F90)
+% Chapter 2 ——— Time Domain (step.F90)
 % ================================================================
 \chapter{Time Domain (STP) }
@@ -21,5 +21,7 @@
 
 Having defined the continuous equations in Chap.~\ref{PE}, we need now to choose 
-a time discretization. In the present chapter, we provide a general description of the \NEMO 
+a time discretization, a key feature of an ocean model as it exerts a strong influence  
+on the structure of the computer code ($i.e.$ on its flowchart). 
+In the present chapter, we provide a general description of the \NEMO 
 time stepping strategy and the consequences for the order in which the equations are
 solved.
@@ -158,4 +160,8 @@
 \end{equation} 
 
+%%gm
+%%gm   UPDATE the next paragraphs with time varying thickness ...
+%%gm
+
 This scheme is rather time consuming since it requires a matrix inversion, 
 but it becomes attractive since a value of 3 or more is needed for N in
@@ -188,8 +194,10 @@
 
 % -------------------------------------------------------------------------------------------------------------
-%        Hydrostatic Pressure gradient
-% -------------------------------------------------------------------------------------------------------------
-\section{Hydrostatic Pressure Gradient --- semi-implicit scheme}
-\label{STP_hpg_imp}
+%        Surface Pressure gradient
+% -------------------------------------------------------------------------------------------------------------
+\section{Surface Pressure Gradient}
+\label{STP_spg_ts}
+
+===>>>>  TO BE written....  :-)
 
 %\gmcomment{ 
@@ -209,34 +217,4 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 %}
-
-The range of stability of the Leap-Frog scheme can be extended by a factor of two
-by introducing a semi-implicit computation of the hydrostatic pressure gradient term
-\citep{Brown_Campana_MWR78}. Instead of evaluating the pressure at $t$, a linear 
-combination of values at $t-\rdt$, $t$ and $t+\rdt$ is used (see \S~\ref{DYN_hpg_imp}).  
-This technique, controlled by the \np{nn\_dynhpg\_rst} namelist parameter, does not 
-introduce a significant additional computational cost when tracers and thus density 
-is time stepped before the dynamics. This time step ordering is used in \NEMO 
-(Fig.\ref{Fig_TimeStep_flowchart}).
-
-
-This technique, used in several GCMs (\NEMO, POP or MOM for instance), 
-makes the Leap-Frog scheme as efficient 
-\footnote{The efficiency is defined as the maximum allowed Courant number of the time 
-stepping scheme divided by the number of computations of the right-hand side per time step.} 
-as the Forward-Backward scheme used in MOM \citep{Griffies_al_OS05} and more 
-efficient than the LF-AM3 scheme (leapfrog time stepping combined with a third order
-Adams-Moulthon interpolation for the predictor phase) used in ROMS 
-\citep{Shchepetkin_McWilliams_OM05}. 
-
-In fact, this technique is efficient when the physical phenomenon that 
-limits the time-step is internal gravity waves (IGWs). Indeed, it is 
-equivalent to applying a time filter to the pressure gradient to eliminate high 
-frequency IGWs. Obviously, the doubling of the time-step is achievable only 
-if no other factors control the time-step, such as the stability limits associated 
-with advection, diffusion or Coriolis terms. For example, it is inefficient in low resolution
-global ocean configurations, since inertial oscillations in the vicinity of the North Pole 
-are the limiting factor for the time step. It is also often inefficient in very high 
-resolution configurations where strong currents and small grid cells exert 
-the strongest constraint on the time step.
 
 % -------------------------------------------------------------------------------------------------------------
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 6040)
@@ -1,4 +1,4 @@
 % ================================================================
-% Chapter 1 � Ocean Tracers (TRA)
+% Chapter 1 ——— Ocean Tracers (TRA)
 % ================================================================
 \chapter{Ocean Tracers (TRA)}
@@ -40,9 +40,7 @@
 described in chapters \S\ref{SBC}, \S\ref{LDF} and  \S\ref{ZDF}, respectively. 
 Note that \mdl{tranpc}, the non-penetrative convection module,  although 
-(temporarily) located in the NEMO/OPA/TRA directory, is described with the 
-model vertical physics (ZDF).
-%%%
-\gmcomment{change the position of eosbn2 in the reference code}
-%%%
+located in the NEMO/OPA/TRA directory as it directly modifies the tracer fields, 
+is described with the model vertical physics (ZDF) together with other available 
+parameterization of convection.
 
 In the present chapter we also describe the diagnostic equations used to compute 
@@ -50,12 +48,12 @@
 freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}).
 
-The different options available to the user are managed by namelist logicals or 
-CPP keys. For each equation term \textit{ttt}, the namelist logicals are \textit{ln\_trattt\_xxx}, 
+The different options available to the user are managed by namelist logicals or CPP keys. 
+For each equation term  \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx}, 
 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 
-The CPP key (when it exists) is \textbf{key\_trattt}. The equivalent code can be 
-found in the \textit{trattt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory.
+The CPP key (when it exists) is \textbf{key\_traTTT}. The equivalent code can be 
+found in the \textit{traTTT} or \textit{traTTT\_xxx} module, in the NEMO/OPA/TRA directory.
 
 The user has the option of extracting each tendency term on the rhs of the tracer 
-equation for output (\key{trdtra} is defined), as described in Chap.~\ref{MISC}.
+equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{MISC}.
 
 $\ $\newline    % force a new ligne
@@ -82,10 +80,10 @@
 implicitly requires the use of the continuity equation. Indeed, it is obtained
 by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 
-which results from the use of the continuity equation, $\nabla \cdot \vect{U}=0$ or 
-$ \partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ in constant volume or variable volume case, respectively. 
+which results from the use of the continuity equation,  $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ 
+(which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, $i.e.$ \np{ln\_linssh}=true). 
 Therefore it is of paramount importance to design the discrete analogue of the 
 advection tendency so that it is consistent with the continuity equation in order to 
 enforce the conservation properties of the continuous equations. In other words, 
-by replacing $\tau$ by the number 1 in (\ref{Eq_tra_adv}) we recover the discrete form of 
+by setting $\tau = 1$ in (\ref{Eq_tra_adv}) we recover the discrete form of 
 the continuity equation which is used to calculate the vertical velocity.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -113,5 +111,5 @@
 boundary condition depends on the type of sea surface chosen: 
 \begin{description}
-\item [linear free surface:] the first level thickness is constant in time: 
+\item [linear free surface:] (\np{ln\_linssh}=true) the first level thickness is constant in time: 
 the vertical boundary condition is applied at the fixed surface $z=0$ 
 rather than on the moving surface $z=\eta$. There is a non-zero advective 
@@ -119,5 +117,5 @@
 $\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $, $i.e.$ 
 the product of surface velocity (at $z=0$) by the first level tracer value.
-\item [non-linear free surface:] (\key{vvl} is defined) 
+\item [non-linear free surface:] (\np{ln\_linssh}=false) 
 convergence/divergence in the first ocean level moves the free surface 
 up/down. There is no tracer advection through it so that the advective 
@@ -125,29 +123,32 @@
 \end{description}
 In all cases, this boundary condition retains local conservation of tracer. 
-Global conservation is obtained in both rigid-lid and non-linear free surface 
-cases, but not in the linear free surface case. Nevertheless, in the latter
-case, it is achieved to a good approximation since the non-conservative 
+Global conservation is obtained in non-linear free surface case, 
+but \textit{not} in the linear free surface case. Nevertheless, in the latter case, 
+it is achieved to a good approximation since the non-conservative 
 term is the product of the time derivative of the tracer and the free surface 
-height, two quantities that are not correlated (see \S\ref{PE_free_surface}, 
-and also \citet{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}).
+height, two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}.
 
 The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco}) 
-is the centred (\textit{now}) \textit{eulerian} ocean velocity (see Chap.~\ref{DYN}). 
-When eddy induced velocity (\textit{eiv}) parameterisation is used it is the \textit{now} 
-\textit{effective} velocity ($i.e.$ the sum of the eulerian and eiv velocities) which is used.
-
-The choice of an advection scheme is made in the \textit{\ngn{nam\_traadv}} namelist, by 
-setting to \textit{true} one and only one of the logicals \textit{ln\_traadv\_xxx}. The 
+is the centred (\textit{now}) \textit{effective} ocean velocity, $i.e.$ the \textit{eulerian} velocity
+(see Chap.~\ref{DYN}) plus the eddy induced velocity (\textit{eiv}) 
+and/or the mixed layer eddy induced velocity (\textit{eiv})
+when those parameterisations are used (see Chap.~\ref{LDF}).
+
+The choice of an advection scheme is made in the \textit{\ngn{namtra\_adv}} namelist, by 
+setting to \textit{true} one of the logicals \textit{ln\_traadv\_xxx}. The 
 corresponding code can be found in the \textit{traadv\_xxx.F90} module, where 
-\textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. Details 
-of the advection schemes are given below. The choice of an advection scheme 
+\textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 
+By default ($i.e.$ in the reference namelist, \ngn{namelist\_ref}), all the logicals 
+are set to \textit{false}. If the user does not select an advection scheme 
+in the configuration namelist (\ngn{namelist\_cfg}), the tracers will not be advected !
+
+Details of the advection schemes are given below. The choice of an advection scheme 
 is a complex matter which depends on the model physics, model resolution, 
 type of tracer, as well as the issue of numerical cost. 
 
 Note that 
-(1) cen2, cen4 and TVD schemes require an explicit diffusion 
-operator while the other schemes are diffusive enough so that they do not 
-require additional diffusion ; 
-(2) cen2, cen4, MUSCL2, and UBS are not \textit{positive} schemes
+(1) CEN and FCT schemes require an explicit diffusion operator 
+while the other schemes are diffusive enough so that they do not necessarily require additional diffusion ; 
+(2) CEN and UBS are not \textit{positive} schemes
 \footnote{negative values can appear in an initially strictly positive tracer field 
 which is advected}
@@ -163,11 +164,13 @@
 
 % -------------------------------------------------------------------------------------------------------------
+%        2nd and 4th order centred schemes
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [$2^{nd}$ and $4^{th}$ order centred schemes (CEN) (\np{ln\_traadv\_cen})]
+			{$2^{nd}$ and $4^{th}$ order centred schemes (CEN) (\np{ln\_traadv\_cen}=true)}
+\label{TRA_adv_cen}
+
 %        2nd order centred scheme  
-% -------------------------------------------------------------------------------------------------------------
-\subsection   [$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2})]
-			{$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2}=true)}
-\label{TRA_adv_cen2}
-
-In the centred second order formulation, the tracer at velocity points is 
+
+In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points is 
 evaluated as the mean of the two neighbouring $T$-point values. 
 For example, in the $i$-direction :
@@ -176,34 +179,21 @@
 \end{equation}
 
-The scheme is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ 
+CEN2 is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ 
 but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously 
 noisy and must be used in conjunction with an explicit diffusion operator to 
 produce a sensible solution. The associated time-stepping is performed using 
 a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in 
-(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value. The centered second 
-order advection is computed in the \mdl{traadv\_cen2} module. In this module,
-it is advantageous to combine the \textit{cen2} scheme with an upstream scheme
-in specific areas which require a strong diffusion in order to avoid the generation 
-of false extrema. These areas are the vicinity of large river mouths, some straits 
-with coarse resolution, and the vicinity of ice cover area ($i.e.$ when the ocean 
-temperature is close to the freezing point).
-This combined scheme has been included for specific grid points in the ORCA2 
-and ORCA4 configurations only. This is an obsolescent feature as the recommended 
-advection scheme for the ORCA configuration is TVD (see  \S\ref{TRA_adv_tvd}).
-
-Note that using the cen2 scheme, the overall tracer advection is of second 
+(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value. 
+CEN2 is computed in the \mdl{traadv\_cen} module.
+
+Note that using the CEN2, the overall tracer advection is of second 
 order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2}) 
 have this order of accuracy. \gmcomment{Note also that ... blah, blah}
 
-% -------------------------------------------------------------------------------------------------------------
 %        4nd order centred scheme  
-% -------------------------------------------------------------------------------------------------------------
-\subsection   [$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4})]
-		     {$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4}=true)}
-\label{TRA_adv_cen4}
-
-In the $4^{th}$ order formulation (to be implemented), tracer values are 
-evaluated at velocity points as a $4^{th}$ order interpolation, and thus depend on 
-the four neighbouring $T$-points. For example, in the $i$-direction:
+
+In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at velocity points as 
+a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points. 
+For example, in the $i$-direction:
 \begin{equation} \label{Eq_tra_adv_cen4}
 \tau _u^{cen4} 
@@ -211,17 +201,18 @@
 \end{equation}
 
-Strictly speaking, the cen4 scheme is not a $4^{th}$ order advection scheme 
+Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme 
 but a $4^{th}$ order evaluation of advective fluxes, since the divergence of 
-advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order. The phrase ``$4^{th}$ 
-order scheme'' used in oceanographic literature is usually associated 
-with the scheme presented here. Introducing a \textit{true} $4^{th}$ order advection 
-scheme is feasible but, for consistency reasons, it requires changes in the 
-discretisation of the tracer advection together with changes in both the 
-continuity equation and the momentum advection terms.  
+advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order. 
+The expression \textit{$4^{th}$ order scheme} used in oceanographic literature 
+is usually associated with the scheme presented here. 
+Introducing a \textit{true} $4^{th}$ order advection scheme is feasible but, 
+for consistency reasons, it requires changes in the discretisation of the tracer 
+advection together with changes in the continuity equation, 
+and the momentum advection and pressure terms.  
 
 A direct consequence of the pseudo-fourth order nature of the scheme is that 
-it is not non-diffusive, i.e. the global variance of a tracer is not preserved using 
-\textit{cen4}. Furthermore, it must be used in conjunction with an explicit 
-diffusion operator to produce a sensible solution. The time-stepping is also 
+it is not non-diffusive, $i.e.$ the global variance of a tracer is not preserved using 
+CEN4. Furthermore, it must be used in conjunction with an explicit 
+diffusion operator to produce a sensible solution. As in CEN2 case, the time-stepping is  
 performed using a leapfrog scheme in conjunction with an Asselin time-filter, 
 so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
@@ -235,14 +226,14 @@
 
 % -------------------------------------------------------------------------------------------------------------
-%        TVD scheme  
-% -------------------------------------------------------------------------------------------------------------
-\subsection   [Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd})]
-			{Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=true)}
+%        FCT scheme  
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [$2^{nd}$ and $4^{th}$ Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct})]
+			{$2^{nd}$ and $4^{th}$ Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct}=true)}
 \label{TRA_adv_tvd}
 
-In the Total Variance Dissipation (TVD) formulation, the tracer at velocity 
+In the Flux Corrected Transport formulation, the tracer at velocity 
 points is evaluated using a combination of an upstream and a centred scheme. 
 For example, in the $i$-direction :
-\begin{equation} \label{Eq_tra_adv_tvd}
+\begin{equation} \label{Eq_tra_adv_fct}
 \begin{split}
 \tau _u^{ups}&= \begin{cases}
@@ -251,5 +242,5 @@
 				  \end{cases}     \\
 \\
-\tau _u^{tvd}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen2} -\tau _u^{ups} } \right)
+\tau _u^{fct}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen} -\tau _u^{ups} } \right)
 \end{split}
 \end{equation}
@@ -260,20 +251,19 @@
 produces a local extremum in the tracer field. The resulting scheme is quite 
 expensive but \emph{positive}. It can be used on both active and passive tracers. 
-This scheme is tested and compared with MUSCL and the MPDATA scheme in 
-\citet{Levy_al_GRL01}; note that in this paper it is referred to as "FCT" (Flux corrected 
-transport) rather than TVD. The TVD scheme is implemented in the \mdl{traadv\_tvd} module.
+This scheme is tested and compared with MUSCL and a MPDATA scheme in \citet{Levy_al_GRL01}. 
+The FCT scheme is implemented in the \mdl{traadv\_fct} module.
 
 For stability reasons (see \S\ref{STP}),
-$\tau _u^{cen2}$ is evaluated  in (\ref{Eq_tra_adv_tvd}) using the \textit{now} tracer while $\tau _u^{ups}$ 
-is evaluated using the \textit{before} tracer. In other words, the advective part of 
-the scheme is time stepped with a leap-frog scheme while a forward scheme is 
-used for the diffusive part. 
+$\tau _u^{cen}$ is evaluated  in (\ref{Eq_tra_adv_fct}) using the \textit{now} tracer 
+while $\tau _u^{ups}$ is evaluated using the \textit{before} tracer. In other words, 
+the advective part of the scheme is time stepped with a leap-frog scheme 
+while a forward scheme is used for the diffusive part. 
 
 % -------------------------------------------------------------------------------------------------------------
 %        MUSCL scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection[MUSCL scheme  (\np{ln\_traadv\_muscl})]
-	{Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_muscl}=T)}
-\label{TRA_adv_muscl}
+\subsection[MUSCL scheme  (\np{ln\_traadv\_mus})]
+	{Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_mus}=T)}
+\label{TRA_adv_mus}
 
 The Monotone Upstream Scheme for Conservative Laws (MUSCL) has been 
@@ -281,5 +271,5 @@
 is evaluated assuming a linear tracer variation between two $T$-points 
 (Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction :
-\begin{equation} \label{Eq_tra_adv_muscl}
+\begin{equation} \label{Eq_tra_adv_mus}
    \tau _u^{mus} = \left\{      \begin{aligned}
          &\tau _i  &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right)
@@ -296,10 +286,6 @@
 
 For an ocean grid point adjacent to land and where the ocean velocity is 
-directed toward land, two choices are available: an upstream flux 
-(\np{ln\_traadv\_muscl}=true) or a second order flux 
-(\np{ln\_traadv\_muscl2}=true). Note that the latter choice does not ensure 
-the \textit{positive} character of the scheme. Only the former can be used 
-on both active and passive tracers. The two MUSCL schemes are implemented 
-in the \mdl{traadv\_tvd} and \mdl{traadv\_tvd2} modules.
+directed toward land, an upstream flux is used. This choice ensure 
+the \textit{positive} character of the scheme. 
 
 % -------------------------------------------------------------------------------------------------------------
@@ -310,7 +296,7 @@
 \label{TRA_adv_ubs}
 
-The UBS advection scheme is an upstream-biased third order scheme based on 
-an upstream-biased parabolic interpolation. It is also known as the Cell 
-Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective 
+The UBS advection scheme (also often called UP3) is an upstream-biased third order 
+scheme based on an upstream-biased parabolic interpolation. It is also known as 
+the Cell Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective 
 Kinematics). For example, in the $i$-direction :
 \begin{equation} \label{Eq_tra_adv_ubs}
@@ -324,21 +310,21 @@
 
 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation 
-error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of the advection 
-scheme is similar to that reported in \cite{Farrow1995}. 
+error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of
+ the advection scheme is similar to that reported in \cite{Farrow1995}. 
 It is a relatively good compromise between accuracy and smoothness. 
 It is not a \emph{positive} scheme, meaning that false extrema are permitted, 
 but the amplitude of such are significantly reduced over the centred second 
-order method. Nevertheless it is not recommended that it should be applied 
-to a passive tracer that requires positivity. 
+or fourth order method. Nevertheless it is not recommended that it should be 
+applied to a passive tracer that requires positivity. 
 
 The intrinsic diffusion of UBS makes its use risky in the vertical direction 
 where the control of artificial diapycnal fluxes is of paramount importance. 
-Therefore the vertical flux is evaluated using the TVD scheme when 
-\np{ln\_traadv\_ubs}=true.
+Therefore the vertical flux is evaluated using either a 2nd order FCT scheme 
+or a 4th order COMPACT scheme (\np{nn\_cen\_v}=2 or 4).
 
 For stability reasons  (see \S\ref{STP}),
-the first term  in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order centred scheme) 
-is evaluated using the \textit{now} tracer (centred in time) while the 
-second term (which is the diffusive part of the scheme), is 
+the first term  in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order 
+centred scheme) is evaluated using the \textit{now} tracer (centred in time) 
+while the second term (which is the diffusive part of the scheme), is 
 evaluated using the \textit{before} tracer (forward in time). 
 This choice is discussed by \citet{Webb_al_JAOT98} in the context of the 
@@ -350,12 +336,16 @@
 substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme.
 
+???
+
 Four different options are possible for the vertical 
 component used in the UBS scheme. $\tau _w^{ubs}$ can be evaluated 
 using either \textit{(a)} a centred $2^{nd}$ order scheme, or  \textit{(b)} 
-a TVD scheme, or  \textit{(c)} an interpolation based on conservative 
+a FCT scheme, or  \textit{(c)} an interpolation based on conservative 
 parabolic splines following the \citet{Shchepetkin_McWilliams_OM05} 
 implementation of UBS in ROMS, or  \textit{(d)} a UBS. The $3^{rd}$ case 
 has dispersion properties similar to an eighth-order accurate conventional scheme.
-The current reference version uses method b)
+The current reference version uses method (b).
+
+???
 
 Note that :
@@ -390,5 +380,5 @@
 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy 
 coefficient which is simply proportional to the velocity:
- $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that NEMO v3.4 still uses 
+ $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note the current version of NEMO still uses 
  \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}.
  %%%
@@ -416,18 +406,4 @@
 direction (as for the UBS case) should be implemented to restore this property.
 
-
-% -------------------------------------------------------------------------------------------------------------
-%        PPM scheme  
-% -------------------------------------------------------------------------------------------------------------
-\subsection   [Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm})]
-			{Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm}=true)}
-\label{TRA_adv_ppm}
-
-The Piecewise Parabolic Method (PPM) proposed by Colella and Woodward (1984) 
-\sgacomment{reference?}
-is based on a quadradic piecewise construction. Like the QCK scheme, it is associated 
-with the ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented 
-in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the reference 
-version 3.3.
 
 % ================================================================
@@ -1167,61 +1143,116 @@
 %        Equation of State
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Equation of State (\np{nn\_eos} = 0, 1 or 2)}
+\subsection{Equation Of Seawater (\np{nn\_eos} = -1, 0, or 1)}
 \label{TRA_eos}
 
-It is necessary to know the equation of state for the ocean very accurately 
-to determine stability properties (especially the Brunt-Vais\"{a}l\"{a} frequency), 
-particularly in the deep ocean. The ocean seawater volumic mass, $\rho$, 
-abusively called density, is a non linear empirical function of \textit{in situ} 
-temperature, salinity and pressure. The reference equation of state is that 
-defined by the Joint Panel on Oceanographic Tables and Standards 
-\citep{UNESCO1983}. It was the standard equation of state used in early 
-releases of OPA. However, even though this computation is fully vectorised, 
-it is quite time consuming ($15$ to $20${\%} of the total CPU time) since 
-it requires the prior computation of the \textit{in situ} temperature from the 
-model \textit{potential} temperature using the \citep{Bryden1973} polynomial 
-for adiabatic lapse rate and a $4^th$ order Runge-Kutta integration scheme. 
-Since OPA6, we have used the \citet{JackMcD1995} equation of state for 
-seawater instead. It allows the computation of the \textit{in situ} ocean density 
-directly as a function of \textit{potential} temperature relative to the surface 
-(an \NEMO variable), the practical salinity (another \NEMO variable) and the 
-pressure (assuming no pressure variation along geopotential surfaces, $i.e.$ 
-the pressure in decibars is approximated by the depth in meters). 
-Both the \citet{UNESCO1983} and \citet{JackMcD1995} equations of state 
-have exactly the same except that the values of the various coefficients have 
-been adjusted by \citet{JackMcD1995} in order to directly use the \textit{potential} 
-temperature instead of the \textit{in situ} one. This reduces the CPU time of the 
-\textit{in situ} density computation to about $3${\%} of the total CPU time, 
-while maintaining a quite accurate equation of state.
-
-In the computer code, a \textit{true} density anomaly, $d_a= \rho / \rho_o - 1$, 
-is computed, with $\rho_o$ a reference volumic mass. Called \textit{rau0} 
-in the code, $\rho_o$ is defined in \mdl{phycst}, and a value of $1,035~Kg/m^3$. 
+The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship 
+linking seawater density, $\rho$, to a number of state variables, 
+most typically temperature, salinity and pressure. 
+Because density gradients control the pressure gradient force through the hydrostatic balance, 
+the equation of state provides a fundamental bridge between the distribution of active tracers 
+and the fluid dynamics. Nonlinearities of the EOS are of major importance, in particular 
+influencing the circulation through determination of the static stability below the mixed layer, 
+thus controlling rates of exchange between the atmosphere  and the ocean interior \citep{Roquet_JPO2015}. 
+Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983}) 
+or TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real 
+ocean circulation is attempted \citep{Roquet_JPO2015}. 
+The use of TEOS-10 is highly recommended because 
+\textit{(i)} it is the new official EOS, 
+\textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and 
+\textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature 
+and practical salinity for EOS-980, both variables being more suitable for use as model variables 
+\citep{TEOS10, Graham_McDougall_JPO13}. 
+EOS-80 is an obsolescent feature of the NEMO system, kept only for backward compatibility.
+For process studies, it is often convenient to use an approximation of the EOS. To that purposed, 
+a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available.
+
+In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$, 
+is computed, with $\rho_o$ a reference density. Called \textit{rau0} 
+in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$. 
 This is a sensible choice for the reference density used in a Boussinesq ocean 
 climate model, as, with the exception of only a small percentage of the ocean, 
-density in the World Ocean varies by no more than 2$\%$ from $1,035~kg/m^3$ 
-\citep{Gill1982}.
-
-Options are defined through the  \ngn{nameos} namelist variables.
-The default option (namelist parameter \np{nn\_eos}=0) is the \citet{JackMcD1995} 
-equation of state. Its use is highly recommended. However, for process studies, 
-it is often convenient to use a linear approximation of the density.
+density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}.
+
+Options are defined through the  \ngn{nameos} namelist variables, and in particular \np{nn\_eos} 
+which controls the EOS used (=-1 for TEOS10 ; =0 for EOS-80 ; =1 for S-EOS).
+\begin{description}
+
+\item[\np{nn\_eos}$=-1$] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used.  
+The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, 
+but it is optimized for a boussinesq fluid and the polynomial expressions have simpler 
+and more computationally efficient expressions for their derived quantities 
+which make them more adapted for use in ocean models. 
+Note that a slightly higher precision polynomial form is now used replacement of the TEOS-10 
+rational function approximation for hydrographic data analysis  \citep{TEOS10}. 
+A key point is that conservative state variables are used: 
+Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: $\degres C$, notation: $\Theta$).
+The pressure in decibars is approximated by the depth in meters. 
+With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. It is set to 
+$C_p=3991.86795711963~J\,Kg^{-1}\,\degres K^{-1}$, according to \citet{TEOS10}.
+
+Choosing polyTEOS10-bsq implies that the state variables used by the model are 
+$\Theta$ and $S_A$. In particular, the initial state deined by the user have to be given as 
+\textit{Conservative} Temperature and \textit{Absolute} Salinity. 
+In addition, setting \np{ln\_useCT} to \textit{true} convert the Conservative SST to potential SST 
+prior to either computing the air-sea and ice-sea fluxes (forced mode) 
+or sending the SST field to the atmosphere (coupled mode).
+
+\item[\np{nn\_eos}$=0$] the polyEOS80-bsq equation of seawater is used.
+It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized 
+to accurately fit EOS80 (Roquet, personal comm.). The state variables used in both the EOS80 
+and the ocean model are: 
+the Practical Salinity ((unit: psu, notation: $S_p$)) and Potential Temperature (unit: $\degres C$, notation: $\theta$).
+The pressure in decibars is approximated by the depth in meters.  
+With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature, 
+salinity and pressure \citep{UNESCO1983}. Nevertheless, a severe assumption is made in order to 
+have a heat content ($C_p T_p$) which is conserved by the model: $C_p$ is set to a constant 
+value, the TEOS10 value. 
+ 
+\item[\np{nn\_eos}$=1$] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen, 
+the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.) 
+(see also \citet{Roquet_JPO2015}). It provides a simplistic linear representation of both 
+cabbeling and thermobaricity effects which is enough for a proper treatment of the EOS 
+in theoretical studies \citep{Roquet_JPO2015}.
 With such an equation of state there is no longer a distinction between 
-\textit{in situ} and \textit{potential} density and both cabbeling and thermobaric
-effects are removed.
-Two linear formulations are available: a function of $T$ only (\np{nn\_eos}=1) 
-and a function of both $T$ and $S$ (\np{nn\_eos}=2):
-\begin{equation} \label{Eq_tra_eos_linear}
+\textit{conservative} and \textit{potential} temperature, as well as between \textit{absolute} 
+and \textit{practical} salinity.
+S-EOS takes the following expression:
+\begin{equation} \label{Eq_tra_S-EOS}
 \begin{split}
-  d_a(T)       &=  \rho (T)      /  \rho_o   - 1     =  \  0.0285         -  \alpha   \;T     \\ 
-  d_a(T,S)    &=  \rho (T,S)   /  \rho_o   - 1     =  \  \beta \; S       -  \alpha   \;T    
+  d_a(T,S,z)  =  ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a  \\
+                                & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a  \\
+                                & - \nu \; T_a \; S_a \;  ) \; / \; \rho_o                     \\
+  with \ \  T_a = T-10  \; ;  & \;  S_a = S-35  \; ;\;  \rho_o = 1026~Kg/m^3
 \end{split}
 \end{equation} 
-where $\alpha$ and $\beta$ are the thermal and haline expansion 
-coefficients, and $\rho_o$, the reference volumic mass, $rau0$. 
-($\alpha$ and $\beta$ can be modified through the \np{rn\_alpha} and 
-\np{rn\_beta} namelist variables). Note that when $d_a$ is a function 
-of $T$ only (\np{nn\_eos}=1), the salinity is a passive tracer and can be 
-used as such.
+where the computer name of the coefficients as well as their standard value are given in \ref{Tab_SEOS}.
+In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing 
+the associated coefficients. 
+Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS.
+setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS.
+Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S.
+
+\end{description}
+
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{table}[!tb]
+\begin{center} \begin{tabular}{|p{26pt}|p{72pt}|p{56pt}|p{136pt}|}
+\hline
+coeff.	& computer name   & S-EOS		&  description		                  \\ \hline
+$a_0$       & \np{nn\_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff. 	\\ \hline
+$b_0$	      & \np{nn\_b0} 	   & 7.6554 $10^{-1}$ &  linear haline  expansion coeff. 	\\ \hline
+$\lambda_1$	& \np{nn\_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$ 	      \\ \hline
+$\lambda_2$	& \np{nn\_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$	 	   \\ \hline
+$\nu$       & \np{nn\_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$ 	   \\ \hline
+$\mu_1$     & \np{nn\_mu1} 	& 1.4970 $10^{-4}$ &  thermobaric coeff. in T    	   \\ \hline
+$\mu_2$     & \np{nn\_mu2} 	& 1.1090 $10^{-5}$ &  thermobaric coeff. in S   	      \\ \hline
+\end{tabular}
+\caption{ \label{Tab_SEOS}
+Standard value of S-EOS coefficients. }
+\end{center}
+\end{table}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1232,50 +1263,18 @@
 
 An accurate computation of the ocean stability (i.e. of $N$, the brunt-Vais\"{a}l\"{a}
- frequency) is of paramount importance as it is used in several ocean 
- parameterisations (namely TKE, KPP, Richardson number dependent 
- vertical diffusion, enhanced vertical diffusion, non-penetrative convection, 
- iso-neutral diffusion). In particular, one must be aware that $N^2$ has to 
- be computed with an \textit{in situ} reference. The expression for $N^2$ 
- depends on the type of equation of state used (\np{nn\_eos} namelist parameter).
-
-For \np{nn\_eos}=0 (\citet{JackMcD1995} equation of state), the \citet{McDougall1987} 
-polynomial expression is used (with the pressure in decibar approximated by 
-the depth in meters): 
+ frequency) is of paramount importance as determine the ocean stratification and 
+ is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent 
+ vertical diffusion, enhanced vertical diffusion, non-penetrative convection, tidal mixing 
+ parameterisation, iso-neutral diffusion). In particular, $N^2$ has to be computed at the local pressure 
+ (pressure in decibar being approximated by the depth in meters). The expression for $N^2$ 
+ is given by: 
 \begin{equation} \label{Eq_tra_bn2}
-N^2 = \frac{g}{e_{3w}} \; \beta   \ 
-	   \left(  \alpha / \beta \ \delta_{k+1/2}[T]     - \delta_{k+1/2}[S]   \right) 
-\end{equation} 
-where $\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 
-They are a function of  $\overline{T}^{\,k+1/2},\widetilde{S}=\overline{S}^{\,k+1/2} - 35.$, 
-and  $z_w$, with $T$ the \textit{potential} temperature and $\widetilde{S}$ a salinity anomaly. 
-Note that both $\alpha$ and $\beta$ depend on \textit{potential} 
-temperature and salinity which are averaged at $w$-points prior 
-to the computation instead of being computed at $T$-points and 
-then averaged to $w$-points.
-
-When a linear equation of state is used (\np{nn\_eos}=1 or 2, 
-\eqref{Eq_tra_bn2} reduces to:
-\begin{equation} \label{Eq_tra_bn2_linear}
 N^2 =	\frac{g}{e_{3w}} \left(   \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T]   \right)
 \end{equation} 
-where $\alpha$ and $\beta $ are the constant coefficients used to 
-defined the linear equation of state \eqref{Eq_tra_eos_linear}.
-
-% -------------------------------------------------------------------------------------------------------------
-%        Specific Heat
-% -------------------------------------------------------------------------------------------------------------
-\subsection   	[Specific Heat (\textit{phycst})]
-			{Specific Heat (\mdl{phycst})}
-\label{TRA_adv_ldf}
-
-The specific heat of sea water, $C_p$, is a function of temperature, salinity 
-and pressure \citep{UNESCO1983}. It is only used in the model to convert 
-surface heat fluxes into surface temperature increase and so the pressure 
-dependence is neglected. The dependence on $T$ and $S$ is weak. 
-For example, with $S=35~psu$, $C_p$ increases from $3989$ to $4002$ 
-when $T$ varies from -2~\degres C to 31~\degres C. Therefore, $C_p$ has 
-been chosen as a constant: $C_p=4.10^3~J\,Kg^{-1}\,\degres K^{-1}$. 
-Its value is set in \mdl{phycst} module. 
-
+where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS, 
+and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 
+The coefficients are a polynomial function of temperature, salinity and depth which expression 
+depends on the chosen EOS. They are computed through \textit{eos\_rab}, a \textsc{Fortran} 
+function that can be found in \mdl{eosbn2}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1298,6 +1297,17 @@
 sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent 
 terms in \eqref{Eq_tra_eos_fzp} (last term) have been dropped. The freezing
-point is computed through \textit{tfreez}, a \textsc{Fortran} function that can be found 
+point is computed through \textit{eos\_fzp}, a \textsc{Fortran} function that can be found 
 in \mdl{eosbn2}.  
+
+
+% -------------------------------------------------------------------------------------------------------------
+%        Potential Energy     
+% -------------------------------------------------------------------------------------------------------------
+%\subsection{Potential Energy anomalies}
+%\label{TRA_bn2}
+
+%    =====>>>>> TO BE written
+%
+
 
 % ================================================================
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Introduction.tex
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Introduction.tex	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Introduction.tex	(revision 6040)
@@ -24,15 +24,14 @@
 release 8.2, described in \citet{Madec1998}. This model has been used for a wide 
 range of applications, both regional or global, as a forced ocean model and as a 
-model coupled with the atmosphere. A complete list of references is found on the 
-\NEMO web site. 
+model coupled with the sea-ice and/or the atmosphere.  
 
 This manual is organised in as follows. Chapter~\ref{PE} presents the model basics, 
 $i.e.$ the equations and their assumptions, the vertical coordinates used, and the 
 subgrid scale physics. This part deals with the continuous equations of the model 
-(primitive equations, with potential temperature, salinity and an equation of state). 
+(primitive equations, with temperature, salinity and an equation of seawater). 
 The equations are written in a curvilinear coordinate system, with a choice of vertical 
-coordinates ($z$ or $s$, with the rescaled height coordinate formulation \textit{z*}, or  
-\textit{s*}). Momentum equations are formulated in the vector invariant form or in the 
-flux form. Dimensional units in the meter, kilogram, second (MKS) international system 
+coordinates ($z$, $s$, \textit{z*}, \textit{s*}, $\tilde{z}$, $\tilde{s}$, and a mixture of them). 
+Momentum equations are formulated in the vector invariant form or in the flux form. 
+Dimensional units in the meter, kilogram, second (MKS) international system 
 are used throughout.
 
@@ -79,9 +78,10 @@
 space and time variable coefficient \citet{Treguier1997}. The model has vertical harmonic 
 viscosity and diffusion with a space and time variable coefficient, with options to compute 
-the coefficients with \citet{Blanke1993}, \citet{Large_al_RG94}, \citet{Pacanowski_Philander_JPO81}, 
+the coefficients with \citet{Blanke1993}, \citet{Pacanowski_Philander_JPO81}, 
 or \citet{Umlauf_Burchard_JMS03} mixing schemes.
  \vspace{1cm}
  
- 
+%%gm    To be put somewhere else ....
+
 \noindent CPP keys and namelists are used for inputs to the code.  \newline
 
@@ -112,4 +112,5 @@
  \vspace{1cm}
 
+%%gm  end
 
 Model outputs management and specific online diagnostics are described in chapters~\ref{DIA}.
@@ -249,2 +250,16 @@
 
 
+ \vspace{1cm}
+$\bullet$ The main modifications from NEMO/OPA v3.4 and  v3.6 are :\\
+\begin{enumerate}
+\item ... ; 
+\end{enumerate}
+
+
+ \vspace{1cm}
+$\bullet$ The main modifications from NEMO/OPA v3.6 and  v4.0 are :\\
+\begin{enumerate}
+\item ... ; 
+\end{enumerate}
+
+
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nam_tide
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nam_tide	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nam_tide	(revision 6040)
@@ -1,16 +1,8 @@
 !-----------------------------------------------------------------------
-&nam_tide      !   tide parameters (#ifdef key_tide)
+&nam_tide      !   tide parameters                                      ("key_tide")
 !-----------------------------------------------------------------------
    ln_tide_pot   = .true.   !  use tidal potential forcing
-   clname(1)     =   'M2'   !  name of constituent
-   clname(2)     =   'S2'
-   clname(3)     =   'N2'
-   clname(4)     =   'K1'
-   clname(5)     =   'O1'
-   clname(6)     =   'Q1'
-   clname(7)     =   'M4'
-   clname(8)     =   'K2'
-   clname(9)     =   'P1'
-   clname(10)    =   'Mf'
-   clname(11)    =   'Mm'
+   ln_tide_ramp  = .false.  !
+   rdttideramp   =    0.    !
+   clname(1)     = 'DUMMY'  !  name of constituent - all tidal components must be set in namelist_cfg
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namagrif
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namagrif	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namagrif	(revision 6040)
@@ -6,3 +6,4 @@
    rn_sponge_tra = 2880.   !  coefficient for tracer   sponge layer [m2/s]
    rn_sponge_dyn = 2880.   !  coefficient for dynamics sponge layer [m2/s]
+   ln_chk_bathy  = .FALSE. !
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambbl
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambbl	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambbl	(revision 6040)
@@ -1,4 +1,4 @@
 !-----------------------------------------------------------------------
-&nambbl        !   bottom boundary layer scheme
+&nambbl        !   bottom boundary layer scheme                         ("key_trabbl")
 !-----------------------------------------------------------------------
    nn_bbl_ldf  =  1      !  diffusive bbl (=1)   or not (=0)
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambdy
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambdy	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambdy	(revision 6040)
@@ -2,22 +2,33 @@
 &nambdy        !  unstructured open boundaries                          ("key_bdy")
 !-----------------------------------------------------------------------
-    nb_bdy = 1                            !  number of open boundary sets
+    nb_bdy         = 0                    !  number of open boundary sets
     ln_coords_file = .true.               !  =T : read bdy coordinates from file
     cn_coords_file = 'coordinates.bdy.nc' !  bdy coordinates files
-    ln_mask_file = .false.                !  =T : read mask from file
-    cn_mask_file = ''                     !  name of mask file (if ln_mask_file=.TRUE.)
-    nn_dyn2d      =  2                    !  boundary conditions for barotropic fields
-    nn_dyn2d_dta  =  3                    !  = 0, bdy data are equal to the initial state
+    ln_mask_file   = .false.              !  =T : read mask from file
+    cn_mask_file   = ''                   !  name of mask file (if ln_mask_file=.TRUE.)
+    cn_dyn2d       = 'none'               !
+    nn_dyn2d_dta   =  0                   !  = 0, bdy data are equal to the initial state
                                           !  = 1, bdy data are read in 'bdydata   .nc' files
                                           !  = 2, use tidal harmonic forcing data from files
                                           !  = 3, use external data AND tidal harmonic forcing
-    nn_dyn3d      =  0                    !  boundary conditions for baroclinic velocities
+    cn_dyn3d      =  'none'               !
     nn_dyn3d_dta  =  0                    !  = 0, bdy data are equal to the initial state
-                           !  = 1, bdy data are read in 'bdydata   .nc' files
-    nn_tra        =  1                    !  boundary conditions for T and S
-    nn_tra_dta    =  1                    !  = 0, bdy data are equal to the initial state
-                           !  = 1, bdy data are read in 'bdydata   .nc' files
-    nn_rimwidth  = 10                      !  width of the relaxation zone
-    ln_vol     = .false.                  !  total volume correction (see nn_volctl parameter)
-    nn_volctl  = 1                        !  = 0, the total water flux across open boundaries is zero
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+    cn_tra        =  'none'               !
+    nn_tra_dta    =  0                    !  = 0, bdy data are equal to the initial state
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+    cn_ice_lim      =  'none'             !
+    nn_ice_lim_dta  =  0                  !  = 0, bdy data are equal to the initial state
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+    rn_ice_tem      = 270.                !  lim3 only: arbitrary temperature of incoming sea ice
+    rn_ice_sal      = 10.                 !  lim3 only:      --   salinity           --
+    rn_ice_age      = 30.                 !  lim3 only:      --   age                --
+
+    ln_tra_dmp    =.false.                !  open boudaries conditions for tracers
+    ln_dyn3d_dmp  =.false.                !  open boundary condition for baroclinic velocities
+    rn_time_dmp   =  1.                   ! Damping time scale in days
+    rn_time_dmp_out =  1.                 ! Outflow damping time scale
+    nn_rimwidth   = 10                    !  width of the relaxation zone
+    ln_vol        = .false.               !  total volume correction (see nn_volctl parameter)
+    nn_volctl     = 1                     !  = 0, the total water flux across open boundaries is zero
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambdy_dta
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambdy_dta	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambdy_dta	(revision 6040)
@@ -1,15 +1,24 @@
 !-----------------------------------------------------------------------
-&nambdy_dta      !  open boundaries - external data           ("key_bdy")
+&nambdy_dta    !  open boundaries - external data                       ("key_bdy")
 !-----------------------------------------------------------------------
-!              !   file name    ! frequency (hours) !  variable  ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
-!              !                !  (if <0  months)  !    name    !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      !
-   bn_ssh =     'amm12_bdyT_u2d' ,         24        , 'sossheig' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''    , ''
-   bn_u2d =     'amm12_bdyU_u2d' ,         24        , 'vobtcrtx' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''    , ''
-   bn_v2d =     'amm12_bdyV_u2d' ,         24        , 'vobtcrty' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''    , ''
-   bn_u3d  =    'amm12_bdyU_u3d' ,         24        , 'vozocrtx' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''    , ''
-   bn_v3d  =    'amm12_bdyV_u3d' ,         24        , 'vomecrty' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''    , ''
-   bn_tem  =    'amm12_bdyT_tra' ,         24        , 'votemper' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''    , ''
-   bn_sal  =    'amm12_bdyT_tra' ,         24        , 'vosaline' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''    , ''
-   cn_dir  =    'bdydta/'
-   ln_full_vel = .false.
+!              !  file name      ! frequency (hours) ! variable  ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !                 !  (if <0  months)  !   name    !  (logical)   !  (T/F ) ! 'monthly' ! filename ! pairing  ! filename      !
+   bn_ssh =     'amm12_bdyT_u2d' ,         24        , 'sossheig',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+   bn_u2d =     'amm12_bdyU_u2d' ,         24        , 'vobtcrtx',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+   bn_v2d =     'amm12_bdyV_u2d' ,         24        , 'vobtcrty',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+   bn_u3d  =    'amm12_bdyU_u3d' ,         24        , 'vozocrtx',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+   bn_v3d  =    'amm12_bdyV_u3d' ,         24        , 'vomecrty',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+   bn_tem  =    'amm12_bdyT_tra' ,         24        , 'votemper',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+   bn_sal  =    'amm12_bdyT_tra' ,         24        , 'vosaline',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+! for lim2
+!   bn_frld  =   'amm12_bdyT_ice' ,         24        , 'ileadfra',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+!   bn_hicif =   'amm12_bdyT_ice' ,         24        , 'iicethic',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+!   bn_hsnif =   'amm12_bdyT_ice' ,         24        , 'isnowthi',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+! for lim3
+!   bn_a_i  =    'amm12_bdyT_ice' ,         24        , 'ileadfra',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+!   bn_ht_i =    'amm12_bdyT_ice' ,         24        , 'iicethic',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+!   bn_ht_s =    'amm12_bdyT_ice' ,         24        , 'isnowthi',     .true.   , .false. ,  'daily'  ,    ''    ,   ''     , ''
+
+   cn_dir      =    'bdydta/'   !  root directory for the location of the bulk files
+   ln_full_vel = .false.        !  
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambdy_tide
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambdy_tide	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambdy_tide	(revision 6040)
@@ -2,6 +2,6 @@
 &nambdy_tide     ! tidal forcing at open boundaries
 !-----------------------------------------------------------------------
-   filtide          = 'bdydta/amm12_bdytide_'         !  file name root of tidal forcing files
-   ln_bdytide_2ddta = .false.
-   ln_bdytide_conj  = .false.
+   filtide          = 'bdydta/amm12_bdytide_'   !  file name root of tidal forcing files
+   ln_bdytide_2ddta = .false.   !
+   ln_bdytide_conj  = .false.   ! 
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namberg
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namberg	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namberg	(revision 6040)
@@ -1,6 +1,6 @@
 !-----------------------------------------------------------------------
-&namberg       !   iceberg parameters
+&namberg       !   iceberg parameters                                   (default: No iceberg)
 !-----------------------------------------------------------------------
-      ln_icebergs              = .false.
+      ln_icebergs              = .false.              ! iceberg floats or not
       ln_bergdia               = .true.               ! Calculate budgets
       nn_verbose_level         = 1                    ! Turn on more verbose output if level > 0
@@ -9,8 +9,8 @@
                                                       ! Initial mass required for an iceberg of each class
       rn_initial_mass          = 8.8e7, 4.1e8, 3.3e9, 1.8e10, 3.8e10, 7.5e10, 1.2e11, 2.2e11, 3.9e11, 7.4e11
-                                                      ! Proportion of calving mass to apportion to each class  
+                                                      ! Proportion of calving mass to apportion to each class
       rn_distribution          = 0.24, 0.12, 0.15, 0.18, 0.12, 0.07, 0.03, 0.03, 0.03, 0.02
                                                       ! Ratio between effective and real iceberg mass (non-dim)
-                                                      ! i.e. number of icebergs represented at a point         
+                                                      ! i.e. number of icebergs represented at a point
       rn_mass_scaling          = 2000, 200, 50, 20, 10, 5, 2, 1, 1, 1
                                                       ! thickness of newly calved bergs (m)
@@ -21,14 +21,14 @@
       rn_bits_erosion_fraction = 0.                   ! Fraction of erosion melt flux to divert to bergy bits
       rn_sicn_shift            = 0.                   ! Shift of sea-ice concn in erosion flux (0<sicn_shift<1)
-      ln_passive_mode          = .false.              ! iceberg - ocean decoupling   
+      ln_passive_mode          = .false.              ! iceberg - ocean decoupling
       nn_test_icebergs         =  10                  ! Create test icebergs of this class (-1 = no)
                                                       ! Put a test iceberg at each gridpoint in box (lon1,lon2,lat1,lat2)
       rn_test_box              = 108.0,  116.0, -66.0, -58.0
-      rn_speed_limit           = 0.                   ! CFL speed limit for a berg   
+      rn_speed_limit           = 0.                   ! CFL speed limit for a berg
 
-               ! filename ! freq (hours) ! variable ! time interp. ! clim  !'yearly' or ! weights  ! rotation ! land/sea mask !
-               !          ! (<0  months) !   name   !  (logical)   ! (T/F) ! 'monthly'  ! filename ! pairing  ! filename      !
-      sn_icb =  'calving' ,     -1       , 'calvingmask',  .true.      , .true., 'yearly'   , ' '      , ' '  , ''
-   
-      cn_dir = './' 
+!              ! file name ! frequency (hours) !   variable   ! time interp.   !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !           !  (if <0  months)  !     name     !   (logical)    !  (T/F ) ! 'monthly' ! filename ! pairing  ! filename      !
+      sn_icb =  'calving' ,       -1           , 'calvingmask',  .true.        , .true.  , 'yearly'  , ''       , ''       , ''
+
+      cn_dir = './'
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambfr
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambfr	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nambfr	(revision 6040)
@@ -1,4 +1,4 @@
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
+&nambfr        !   bottom friction                                      (default: linear)
 !-----------------------------------------------------------------------
    nn_bfr      =    1      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdct
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdct	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdct	(revision 6040)
@@ -1,4 +1,4 @@
 !-----------------------------------------------------------------------
-&namdct        ! transports through sections
+&namdct        ! transports through some sections
 !-----------------------------------------------------------------------
     nn_dct      = 15       !  time step frequency for transports computing
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdom
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdom	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdom	(revision 6040)
@@ -3,6 +3,7 @@
 !-----------------------------------------------------------------------
    nn_bathy    =    1      !  compute (=0) or read (=1) the bathymetry file
-   nn_closea    =   0      !  remove (=0) or keep (=1) closed seas and lakes (ORCA)
-   nn_msh      =    0      !  create (=1) a mesh file or not (=0)
+   rn_bathy    =    0.     !  value of the bathymetry. if (=0) bottom flat at jpkm1
+   nn_closea   =    0      !  remove (=0) or keep (=1) closed seas and lakes (ORCA)
+   nn_msh      =    1      !  create (=1) a mesh file or not (=0)
    rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0)
    rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of
@@ -16,4 +17,5 @@
    rn_rdtmax   = 28800.          !  maximum time step on tracers (used if nn_acc=1)
    rn_rdth     =  800.           !  depth variation of tracer time step  (used if nn_acc=1)
+   ln_crs      = .false.      !  Logical switch for coarsening module
    jphgr_msh   =       0               !  type of horizontal mesh
                                        !  = 0 curvilinear coordinate on the sphere read in coordinate.nc
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_adv
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_adv	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_adv	(revision 6040)
@@ -1,7 +1,9 @@
 !-----------------------------------------------------------------------
-&namdyn_adv    !   formulation of the momentum advection
+&namdyn_adv    !   formulation of the momentum advection                (default: vector form)
 !-----------------------------------------------------------------------
-   ln_dynadv_vec = .true.  !  vector form (T) or flux form (F)  
+   ln_dynadv_vec = .true.  !  vector form (T) or flux form (F)
+   nn_dynkeg     = 0       ! scheme for grad(KE): =0   C2  ;  =1   Hollingsworth correction
    ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme
-   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme 
-/  
+   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme
+   ln_dynzad_zts = .false. !  Use (T) sub timestepping for vertical momentum advection
+/
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_hpg
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_hpg	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_hpg	(revision 6040)
@@ -1,12 +1,10 @@
 !-----------------------------------------------------------------------
-&namdyn_hpg    !   Hydrostatic pressure gradient option
+&namdyn_hpg    !   Hydrostatic pressure gradient option                 (default: zps)
 !-----------------------------------------------------------------------
    ln_hpg_zco  = .false.   !  z-coordinate - full steps
    ln_hpg_zps  = .true.    !  z-coordinate - partial steps (interpolation)
    ln_hpg_sco  = .false.   !  s-coordinate (standard jacobian formulation)
-   ln_hpg_isf  = .false.   !  s-coordinate (sco ) adapted to ice shelf cavity
+   ln_hpg_isf  = .false.   !  s-coordinate (sco ) adapted to isf
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
    ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
-   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T)
-                                 !           centered      time scheme  (F)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_ldf
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_ldf	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_ldf	(revision 6040)
@@ -1,20 +1,26 @@
 !-----------------------------------------------------------------------
-&namdyn_ldf    !   lateral diffusion on momentum
+&namdyn_ldf    !   lateral diffusion on momentum                        (default: NO)
 !-----------------------------------------------------------------------
    !                       !  Type of the operator :
-   ln_dynldf_lap    =  .true.   !  laplacian operator
-   ln_dynldf_bilap  =  .false.  !  bilaplacian operator
+   !                           !  no diffusion: set ln_dynldf_lap=..._blp=F 
+   ln_dynldf_lap =  .false.    !    laplacian operator
+   ln_dynldf_blp =  .false.    !  bilaplacian operator
    !                       !  Direction of action  :
-   ln_dynldf_level  =  .false.  !  iso-level
-   ln_dynldf_hor    =  .true.   !  horizontal (geopotential)            (require "key_ldfslp" in s-coord.)
-   ln_dynldf_iso    =  .false.  !  iso-neutral                          (require "key_ldfslp")
+   ln_dynldf_lev =  .false.    !  iso-level
+   ln_dynldf_hor =  .false.    !  horizontal (geopotential)
+   ln_dynldf_iso =  .false.    !  iso-neutral
    !                       !  Coefficient
-   rn_ahm_0_lap     = 40000.    !  horizontal laplacian eddy viscosity   [m2/s]
-   rn_ahmb_0        =     0.    !  background eddy viscosity for ldf_iso [m2/s]
-   rn_ahm_0_blp     =     0.    !  horizontal bilaplacian eddy viscosity [m4/s]
-   rn_cmsmag_1      =     3.    !  constant in laplacian Smagorinsky viscosity
-   rn_cmsmag_2      =     3     !  constant in bilaplacian Smagorinsky viscosity
-   rn_cmsh          =     1.    !  1 or 0 , if 0 -use only shear for Smagorinsky viscosity
-   rn_ahm_m_blp     =    -1.e12 !  upper limit for bilap  abs(ahm) < min( dx^4/128rdt, rn_ahm_m_blp)
-   rn_ahm_m_lap     = 40000.    !  upper limit for lap  ahm < min(dx^2/16rdt, rn_ahm_m_lap)
+   nn_ahm_ijk_t  = 0           !  space/time variation of eddy coef
+   !                                !  =-30  read in eddy_viscosity_3D.nc file
+   !                                !  =-20  read in eddy_viscosity_2D.nc file
+   !                                !  =  0  constant 
+   !                                !  = 10  F(k)=c1d
+   !                                !  = 20  F(i,j)=F(grid spacing)=c2d
+   !                                !  = 30  F(i,j,k)=c2d*c1d
+   !                                !  = 31  F(i,j,k)=F(grid spacing and local velocity)
+   rn_ahm_0      =  40000.     !  horizontal laplacian eddy viscosity   [m2/s]
+   rn_ahm_b      =      0.     !  background eddy viscosity for ldf_iso [m2/s]
+   rn_bhm_0      = 1.e+12      !  horizontal bilaplacian eddy viscosity [m4/s]
+   !
+   ! Caution in 20 and 30 cases the coefficient have to be given for a 1 degree grid (~111km)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_spg
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_spg	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_spg	(revision 6040)
@@ -1,6 +1,14 @@
 !-----------------------------------------------------------------------
-!namdyn_spg    !   surface pressure gradient   (CPP key only)
+&namdyn_spg    !   surface pressure gradient                            (default: NO)
 !-----------------------------------------------------------------------
-!                          !  explicit free surface                     ("key_dynspg_exp")
-!                          !  filtered free surface                     ("key_dynspg_flt")
-!                          !  split-explicit free surface               ("key_dynspg_ts")
+   ln_dynspg_exp  = .false.   ! explicit free surface
+   ln_dynspg_ts   = .false.   ! split-explicit free surface
+      ln_bt_fw      = .true.     ! Forward integration of barotropic Eqs.
+      ln_bt_av      = .true.     ! Time filtering of barotropic variables
+         nn_bt_flt     = 1          ! Time filter choice  = 0 None
+         !                          !                     = 1 Boxcar over   nn_baro sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+      ln_bt_auto    = .true.     ! Number of sub-step defined from:
+         rn_bt_cmax   =  0.8        ! =T : the Maximum Courant Number allowed
+         nn_baro      = 30          ! =F : the number of sub-step in rn_rdt seconds
+/
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_vor
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_vor	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namdyn_vor	(revision 6040)
@@ -1,8 +1,10 @@
 !-----------------------------------------------------------------------
-&namdyn_vor    !   option of physics/algorithm (not control by CPP keys)
+&namdyn_vor    !   option of physics/algorithm                          (default: NO)
 !-----------------------------------------------------------------------
    ln_dynvor_ene = .false. !  enstrophy conserving scheme
    ln_dynvor_ens = .false. !  energy conserving scheme
    ln_dynvor_mix = .false. !  mixed scheme
-   ln_dynvor_een = .true.  !  energy & enstrophy scheme
+   ln_dynvor_een = .false. !  energy & enstrophy scheme
+      nn_een_e3f = 1          ! e3f = masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
+   ln_dynvor_msk = .false. !  vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)  ! PLEASE DO NOT USE
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nameos
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nameos	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/nameos	(revision 6040)
@@ -2,9 +2,18 @@
 &nameos        !   ocean physical parameters
 !-----------------------------------------------------------------------
-   nn_eos      =   0       !  type of equation of state and Brunt-Vaisala frequency
-                           !     = 0, UNESCO (formulation of Jackett and McDougall (1994) and of McDougall (1987) )
-                           !     = 1, linear: rho(T)   = rau0 * ( 1.028 - ralpha * T )
-                           !     = 2, linear: rho(T,S) = rau0 * ( rbeta * S - ralpha * T )
-   rn_alpha    =   2.0e-4  !  thermal expension coefficient (nn_eos= 1 or 2)
-   rn_beta     =   7.7e-4  !  saline  expension coefficient (nn_eos= 2)
+   nn_eos      =  -1     !  type of equation of state and Brunt-Vaisala frequency
+                                 !  =-1, TEOS-10
+                                 !  = 0, EOS-80
+                                 !  = 1, S-EOS   (simplified eos)
+   ln_useCT    = .true.  ! use of Conservative Temp. ==> surface CT converted in Pot. Temp. in sbcssm
+                                 !
+   !                     ! S-EOS coefficients :
+                                 !  rd(T,S,Z)*rau0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
+   rn_a0       =  1.6550e-1      !  thermal expension coefficient (nn_eos= 1)
+   rn_b0       =  7.6554e-1      !  saline  expension coefficient (nn_eos= 1)
+   rn_lambda1  =  5.9520e-2      !  cabbeling coeff in T^2  (=0 for linear eos)
+   rn_lambda2  =  7.4914e-4      !  cabbeling coeff in S^2  (=0 for linear eos)
+   rn_mu1      =  1.4970e-4      !  thermobaric coeff. in T (=0 for linear eos)
+   rn_mu2      =  1.1090e-5      !  thermobaric coeff. in S (=0 for linear eos)
+   rn_nu       =  2.4341e-3      !  cabbeling coeff in T*S  (=0 for linear eos)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namhsb
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namhsb	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namhsb	(revision 6040)
@@ -1,4 +1,4 @@
 !-----------------------------------------------------------------------
-&namhsb       !  Heat and salt budgets 
+&namhsb       !  Heat and salt budgets                                  (default F)
 !-----------------------------------------------------------------------
    ln_diahsb  = .false.    !  check the heat and salt budgets (T) or not (F)
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namlbc
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namlbc	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namlbc	(revision 6040)
@@ -2,6 +2,6 @@
 &namlbc        !   lateral momentum boundary condition
 !-----------------------------------------------------------------------
+   !                       !  free slip  !   partial slip  !   no slip   ! strong slip
    rn_shlat    =    2.     !  shlat = 0  !  0 < shlat < 2  !  shlat = 2  !  2 < shlat
-                           !  free slip  !   partial slip  !   no slip   ! strong slip
-   ln_vorlat   = .false.   !  consistency of vorticity boundary condition with analytical eqs.
+   ln_vorlat   = .false.   !  consistency of vorticity boundary condition with analytical Eqs.
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namobs
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namobs	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namobs	(revision 6040)
@@ -5,59 +5,51 @@
    ln_s3d     = .false.    ! Logical switch for S profile observations
    ln_ena     = .false.    ! Logical switch for ENACT insitu data set
-   !                       !     ln_cor                  Logical switch for Coriolis insitu data set
+   ln_cor     = .false.    ! Logical switch for Coriolis insitu data set
    ln_profb   = .false.    ! Logical switch for feedback insitu data set
    ln_sla     = .false.    ! Logical switch for SLA observations
-
    ln_sladt   = .false.    ! Logical switch for AVISO SLA data
-
    ln_slafb   = .false.    ! Logical switch for feedback SLA data
-                           !     ln_ssh                  Logical switch for SSH observations
-
-   ln_sst     = .true.     ! Logical switch for SST observations
-   ln_reysst  = .true.     !     ln_reysst               Logical switch for Reynolds observations
-   ln_ghrsst  = .false.    !     ln_ghrsst               Logical switch for GHRSST observations      
-
+   ln_ssh     = .false.    ! Logical switch for SSH observations
+   ln_sst     = .false.    ! Logical switch for SST observations
+   ln_reysst  = .false.    ! Logical switch for Reynolds observations
+   ln_ghrsst  = .false.    ! Logical switch for GHRSST observations
    ln_sstfb   = .false.    ! Logical switch for feedback SST data
-                           !     ln_sss                  Logical switch for SSS observations
-                           !     ln_seaice               Logical switch for Sea Ice observations
-                           !     ln_vel3d                Logical switch for velocity observations
-                           !     ln_velavcur             Logical switch for velocity daily av. cur.
-                           !     ln_velhrcur             Logical switch for velocity high freq. cur.
-                           !     ln_velavadcp            Logical switch for velocity daily av. ADCP
-                           !     ln_velhradcp            Logical switch for velocity high freq. ADCP
-                           !     ln_velfb                Logical switch for feedback velocity data
-                           !     ln_grid_global          Global distribtion of observations
-                           !     ln_grid_search_lookup   Logical switch for obs grid search w/lookup table
-                           !     grid_search_file        Grid search lookup file header
-                           !     enactfiles              ENACT input observation file names
-                           !     coriofiles              Coriolis input observation file name
-   !                       ! profbfiles: Profile feedback input observation file name
-   profbfiles = 'profiles_01.nc'
-                           !     ln_profb_enatim         Enact feedback input time setting switch
-                           !     slafilesact             Active SLA input observation file name
-                           !     slafilespas             Passive SLA input observation file name
-   !                       ! slafbfiles: Feedback SLA input observation file name
-   slafbfiles = 'sla_01.nc'
-                           !     sstfiles                GHRSST input observation file name
-   !                       ! sstfbfiles: Feedback SST input observation file name
-   sstfbfiles = 'sst_01.nc' 'sst_02.nc' 'sst_03.nc' 'sst_04.nc' 'sst_05.nc'
-                           !     seaicefiles             Sea Ice input observation file name
-                           !     velavcurfiles           Vel. cur. daily av. input file name
-                           !     velhvcurfiles           Vel. cur. high freq. input file name
-                           !     velavadcpfiles          Vel. ADCP daily av. input file name
-                           !     velhvadcpfiles          Vel. ADCP high freq. input file name
-                           !     velfbfiles              Vel. feedback input observation file name
-                           !     dobsini                 Initial date in window YYYYMMDD.HHMMSS
-                           !     dobsend                 Final date in window YYYYMMDD.HHMMSS
-                           !     n1dint                  Type of vertical interpolation method
-                           !     n2dint                  Type of horizontal interpolation method
-                           !     ln_nea                  Rejection of observations near land switch
-   nmsshc     = 0          ! MSSH correction scheme
-                           !     mdtcorr                 MDT  correction
-                           !     mdtcutoff               MDT cutoff for computed correction
+   ln_sss     = .false.    ! Logical switch for SSS observations
+   ln_seaice  = .false.    ! Logical switch for Sea Ice observations
+   ln_vel3d   = .false.    ! Logical switch for velocity observations
+   ln_velavcur= .false     ! Logical switch for velocity daily av. cur.
+   ln_velhrcur= .false     ! Logical switch for velocity high freq. cur.
+   ln_velavadcp=.false.    ! Logical switch for velocity daily av. ADCP
+   ln_velhradcp=.false.    ! Logical switch for velocity high freq. ADCP
+   ln_velfb   = .false.    ! Logical switch for feedback velocity data
+   ln_grid_global=.false.  ! Global distribtion of observations
+   ln_grid_search_lookup = .false. !  Logical switch for obs grid search w/lookup table
+   grid_search_file = 'grid_search'  !  Grid search lookup file header
+! All of the *files* variables below are arrays. Use namelist_cfg to add more files
+   enactfiles = 'enact.nc' !  ENACT input observation file names (specify full array in namelist_cfg)
+   coriofiles = 'corio.nc' !  Coriolis input observation file name
+   profbfiles = 'profiles_01.nc' ! Profile feedback input observation file name
+   ln_profb_enatim = .false !        Enact feedback input time setting switch
+   slafilesact = 'sla_act.nc' !  Active SLA input observation file names
+   slafilespas = 'sla_pass.nc' ! Passive SLA input observation file names
+   slafbfiles = 'sla_01.nc' ! slafbfiles: Feedback SLA input observation file names
+   sstfiles = 'ghrsst.nc'   ! GHRSST input observation file names
+   sstfbfiles = 'sst_01.nc' ! Feedback SST input observation file names
+   seaicefiles = 'seaice_01.nc' ! Sea Ice input observation file names
+   velavcurfiles  = 'velavcurfile.nc'  ! Vel. cur. daily av. input file name
+   velhrcurfiles  = 'velhrcurfile.nc'  ! Vel. cur. high freq. input file name
+   velavadcpfiles = 'velavadcpfile.nc' ! Vel. ADCP daily av. input file name
+   velhradcpfiles = 'velhradcpfile.nc' ! Vel. ADCP high freq. input file name
+   velfbfiles = 'velfbfile.nc' ! Vel. feedback input observation file name
+   dobsini = 20000101.000000  !  Initial date in window YYYYMMDD.HHMMSS
+   dobsend = 20010101.000000  !  Final date in window YYYYMMDD.HHMMSS
+   n1dint = 0  !               Type of vertical interpolation method
+   n2dint = 0  !               Type of horizontal interpolation method
+   ln_nea = .false.   !        Rejection of observations near land switch
+   nmsshc     = 0     !        MSSH correction scheme
+   mdtcorr = 1.61     !        MDT  correction
+   mdtcutoff = 65.0   !        MDT cutoff for computed correction
    ln_altbias = .false.    ! Logical switch for alt bias
    ln_ignmis  = .true.     ! Logical switch for ignoring missing files
-                           !     endailyavtypes   ENACT daily average types
-   ln_grid_global = .true.
-   ln_grid_search_lookup = .false.
+   endailyavtypes = 820    ! ENACT daily average types - array (use namelist_cfg to set more values)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namptr
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namptr	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namptr	(revision 6040)
@@ -1,11 +1,6 @@
 !-----------------------------------------------------------------------
-&namptr       !   Poleward Transport Diagnostic
+&namptr       !   Poleward Transport Diagnostic                         (default F)
 !-----------------------------------------------------------------------
    ln_diaptr  = .false.    !  Poleward heat and salt transport (T) or not (F)
-   ln_diaznl  = .true.     !  Add zonal means and meridional stream functions
-   ln_subbas  = .true.     !  Atlantic/Pacific/Indian basins computation (T) or not
-                           !  (orca configuration only, need input basins mask file named "subbasins.nc"
-   ln_ptrcomp = .true.     !  Add decomposition : overturning
-   nn_fptr    =  1         !  Frequency of ptr computation [time step]
-   nn_fwri    =  15        !  Frequency of ptr outputs [time step]
+   ln_subbas  = .false.     !  Atlantic/Pacific/Indian basins computation (T) or not
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namrun
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namrun	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namrun	(revision 6040)
@@ -9,16 +9,22 @@
    nn_leapy    =       0   !  Leap year calendar (1) or not (0)
    ln_rstart   = .false.   !  start from rest (F) or from a restart file (T)
-   nn_rstctl   =       0   !  restart control => activated only if ln_rstart = T
-                           !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
-                           !    = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart
-                           !    = 2 nn_date0 read in restart  ; nn_it000 : check consistancy between namelist and restart
-   cn_ocerst_in  = "restart"   !  suffix of ocean restart name (input)
-   cn_ocerst_out = "restart"   !  suffix of ocean restart name (output)
+      nn_euler    =    1         !  = 0 : start with forward time step if ln_rstart=T
+      nn_rstctl   =    0         !  restart control ==> activated only if ln_rstart=T
+                                 !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
+                                 !    = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart
+                                 !    = 2 nn_date0 read in restart  ; nn_it000 : check consistancy between namelist and restart
+      cn_ocerst_in  = "restart"  !  suffix of ocean restart name (input)
+      cn_ocerst_indir = "."      !  directory from which to read input ocean restarts
+      cn_ocerst_out = "restart"  !  suffix of ocean restart name (output)
+      cn_ocerst_outdir = "."     !  directory in which to write output ocean restarts
    nn_istate   =       0   !  output the initial state (1) or not (0)
+   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F)
    nn_stock    =    5475   !  frequency of creation of a restart file (modulo referenced to 1)
+   nn_stocklist = 0,0,0,0,0,0,0,0,0,0 ! List of timesteps when a restart file is to be written
    nn_write    =    5475   !  frequency of write in the output file   (modulo referenced to nn_it000)
    ln_dimgnnn  = .false.   !  DIMG file format: 1 file for all processors (F) or by processor (T)
    ln_mskland  = .false.   !  mask land points in NetCDF outputs (costly: + ~15%)
-   ln_clobber  = .false.   !  clobber (overwrite) an existing file
+   ln_cfmeta   = .false.   !  output additional data to netCDF files required for compliance with the CF metadata standard
+   ln_clobber  = .true.    !  clobber (overwrite) an existing file
    nn_chunksz  =       0   !  chunksize (bytes) for NetCDF file (works only with iom_nf90 routines)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc	(revision 6040)
@@ -3,5 +3,6 @@
 !-----------------------------------------------------------------------
    nn_fsbc     = 5         !  frequency of surface boundary condition computation
-                           !     (also = the frequency of sea-ice model call)
+                           !     (also = the frequency of sea-ice & iceberg model call)
+                     ! Type of air-sea fluxes 
    ln_ana      = .false.   !  analytical formulation                    (T => fill namsbc_ana )
    ln_flx      = .false.   !  flux formulation                          (T => fill namsbc_flx )
@@ -9,14 +10,32 @@
    ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core)
    ln_blk_mfs  = .false.   !  MFS bulk formulation                      (T => fill namsbc_mfs )
-   ln_cpl      = .false.   !  Coupled formulation                       (T => fill namsbc_cpl )
-   ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
+                     ! Type of coupling (Ocean/Ice/Atmosphere) :
+   ln_cpl      = .false.   !  atmosphere coupled   formulation          ( requires key_oasis3 )
+   ln_mixcpl   = .false.   !  forced-coupled mixed formulation          ( requires key_oasis3 )
+   nn_components = 0       !  configuration of the opa-sas OASIS coupling
+                           !  =0 no opa-sas OASIS coupling: default single executable configuration
+                           !  =1 opa-sas OASIS coupling: multi executable configuration, OPA component
+                           !  =2 opa-sas OASIS coupling: multi executable configuration, SAS component 
+   nn_limflx = -1          !  LIM3 Multi-category heat flux formulation (use -1 if LIM3 is not used)
+                           !  =-1  Use per-category fluxes, bypass redistributor, forced mode only, not yet implemented coupled
+                           !  = 0  Average per-category fluxes (forced and coupled mode)
+                           !  = 1  Average and redistribute per-category fluxes, forced mode only, not yet implemented coupled
+                           !  = 2  Redistribute a single flux over categories (coupled mode only)
+                     ! Sea-ice :
    nn_ice      = 2         !  =0 no ice boundary condition   ,
                            !  =1 use observed ice-cover      ,
-                           !  =2 ice-model used                         ("key_lim3" or "key_lim2)
-   nn_ice_embd = 0         !  =0 levitating ice (no mass exchange, concentration/dilution effect)
+                           !  =2 ice-model used                         ("key_lim3", "key_lim2", "key_cice")
+   nn_ice_embd = 1         !  =0 levitating ice (no mass exchange, concentration/dilution effect)
                            !  =1 levitating ice with mass and salt exchange but no presure effect
                            !  =2 embedded sea-ice (full salt and mass exchanges and pressure)
-   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave
+                     ! Misc. options of sbc : 
+   ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr )
+      ln_dm2dc = .false.      !  daily mean to diurnal cycle on short wave
    ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf)
+   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr)
+   nn_fwb      = 2         !  FreshWater Budget: =0 unchecked
+                           !     =1 global mean of e-p-r set to zero at each time step
+                           !     =2 annual global mean of e-p-r set to zero
+   ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
    nn_isf      = 0         !  ice shelf melting/freezing                (/=0 => fill namsbc_isf)
                            !  0 =no isf                  1 = presence of ISF
@@ -24,14 +43,6 @@
                            !  4 = ISF fwf specified
                            !  option 1 and 4 need ln_isfcav = .true. (domzgr)
-   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr)
-   nn_fwb      = 3         !  FreshWater Budget: =0 unchecked
-                           !     =1 global mean of e-p-r set to zero at each time step
-                           !     =2 annual global mean of e-p-r set to zero
-                           !     =3 global emp set to zero and spread out over erp area
-   ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave)
-   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
-   ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
-   nn_lsm  = 0             !  =0   land/sea mask for input fields is not applied (the field land/sea mask filename
-                           !       is left empty in namelist)  ,
-                           !  =1:n number of iterations of land/sea mask application for input fields
+   ln_wave = .false.       !  coupling with surface wave                    (T => fill namsbc_wave)
+   nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) ,
+                           !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_apr
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_apr	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_apr	(revision 6040)
@@ -1,11 +1,11 @@
 !-----------------------------------------------------------------------
-&namsbc_apr    !   Atmospheric pressure used as ocean forcing or in bulk
+&namsbc_apr    !   Atmospheric pressure forcing (in ocean or bulk)      (ln_apr_dyn=T)
 !-----------------------------------------------------------------------
-!              ! file name ! frequency (hours) ! variable ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
-!              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_apr      = 'patm'    ,         -1        ,'somslpre',    .true.      , .true.  , 'yearly'  ,  ''      ,   ''     , ''
+!              !  file name ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !            !  (if <0  months)  !   name    !  (logical)   !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_apr      = 'patm'     ,         -1        ,'somslpre',    .true.     , .true. , 'yearly'  ,  ''      ,   ''     , ''
 
    cn_dir      = './'       !  root directory for the location of the bulk files
-   rn_pref     = 101000._wp !  reference atmospheric pressure   [N/m2]/
+   rn_pref     = 101000.    !  reference atmospheric pressure   [N/m2]/
    ln_ref_apr  = .false.    !  ref. pressure: global mean Patm (T) or a constant (F)
    ln_apr_obc  = .false.    !  inverse barometer added to OBC ssh data
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_core
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_core	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_core	(revision 6040)
@@ -2,19 +2,23 @@
 &namsbc_core   !   namsbc_core  CORE bulk formulae
 !-----------------------------------------------------------------------
-!              !  file name                    ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
-!              !                               !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !  filename     !
-   sn_wndi     = 'u_10.15JUNE2009_orca2'       ,         6         , 'U_10_MOD',   .false.    , .true. , 'yearly'  , ''       , 'Uwnd'   , ''
-   sn_wndj     = 'v_10.15JUNE2009_orca2'       ,         6         , 'V_10_MOD',   .false.    , .true. , 'yearly'  , ''       , 'Vwnd'   , ''
-   sn_qsr      = 'ncar_rad.15JUNE2009_orca2'   ,        24         , 'SWDN_MOD',   .false.    , .true. , 'yearly'  , ''       , ''       , ''
-   sn_qlw      = 'ncar_rad.15JUNE2009_orca2'   ,        24         , 'LWDN_MOD',   .false.    , .true. , 'yearly'  , ''       , ''       , ''
-   sn_tair     = 't_10.15JUNE2009_orca2'       ,         6         , 'T_10_MOD',   .false.    , .true. , 'yearly'  , ''       , ''       , ''
-   sn_humi     = 'q_10.15JUNE2009_orca2'       ,         6         , 'Q_10_MOD',   .false.    , .true. , 'yearly'  , ''       , ''       , ''
-   sn_prec     = 'ncar_precip.15JUNE2009_orca2',        -1         , 'PRC_MOD1',   .false.    , .true. , 'yearly'  , ''       , ''       , ''
-   sn_snow     = 'ncar_precip.15JUNE2009_orca2',        -1         , 'SNOW'    ,   .false.    , .true. , 'yearly'  , ''       , ''       , ''
-   sn_tdif     = 'taudif_core'                 ,        24         , 'taudif'  ,   .false.    , .true. , 'yearly'  , ''       , ''       , ''
-!
+!              !  file name                   ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights                               ! rotation ! land/sea mask !
+!              !                              !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename                              ! pairing  ! filename      !
+   sn_wndi     = 'u_10.15JUNE2009_fill'       ,         6         , 'U_10_MOD',   .false.    , .true. , 'yearly'  , 'weights_core_orca2_bicubic_noc.nc'   , 'Uwnd'   , ''
+   sn_wndj     = 'v_10.15JUNE2009_fill'       ,         6         , 'V_10_MOD',   .false.    , .true. , 'yearly'  , 'weights_core_orca2_bicubic_noc.nc'   , 'Vwnd'   , ''
+   sn_qsr      = 'ncar_rad.15JUNE2009_fill'   ,        24         , 'SWDN_MOD',   .false.    , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc'  , ''       , ''
+   sn_qlw      = 'ncar_rad.15JUNE2009_fill'   ,        24         , 'LWDN_MOD',   .false.    , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc'  , ''       , ''
+   sn_tair     = 't_10.15JUNE2009_fill'       ,         6         , 'T_10_MOD',   .false.    , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc'  , ''       , ''
+   sn_humi     = 'q_10.15JUNE2009_fill'       ,         6         , 'Q_10_MOD',   .false.    , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc'  , ''       , ''
+   sn_prec     = 'ncar_precip.15JUNE2009_fill',        -1         , 'PRC_MOD1',   .false.    , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc'  , ''       , ''
+   sn_snow     = 'ncar_precip.15JUNE2009_fill',        -1         , 'SNOW'    ,   .false.    , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc'  , ''       , ''
+   sn_tdif     = 'taudif_core'                ,        24         , 'taudif'  ,   .false.    , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc'  , ''       , ''
+
    cn_dir      = './'      !  root directory for the location of the bulk files
-   ln_2m       = .false.   !  air temperature and humidity referenced at 2m (T) instead 10m (F)
    ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data
+   rn_zqt      = 10.       !  Air temperature and humidity reference height (m)
+   rn_zu       = 10.       !  Wind vector reference height (m)
    rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
+   rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.)
+   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity
+                           !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_cpl
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_cpl	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_cpl	(revision 6040)
@@ -1,23 +1,27 @@
 !-----------------------------------------------------------------------
-&namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_coupled")
+&namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_oasis3")
 !-----------------------------------------------------------------------
-!                    !     description       !  multiple  !    vector   !      vector          ! vector !
-!                    !                       ! categories !  reference  !    orientation       ! grids  !
+!                    !     description      !  multiple  !    vector   !      vector          ! vector !
+!                    !                      ! categories !  reference  !    orientation       ! grids  !
 ! send
-sn_snd_temp   =       'weighted oce and ice' ,    'no'    ,     ''      ,         ''           ,   ''
-sn_snd_alb    =       'weighted ice'         ,    'no'    ,     ''      ,         ''           ,   ''
-sn_snd_thick  =       'none'                 ,    'no'   ,     ''      ,         ''           ,   ''
-sn_snd_crt    =       'none'                 ,    'no'    , 'spherical' , 'eastward-northward' ,  'T'
-sn_snd_co2    =       'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_temp   =   'weighted oce and ice' ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_alb    =   'weighted ice'         ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_thick  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_crt    =   'none'                 ,    'no'    , 'spherical' , 'eastward-northward' ,  'T'
+   sn_snd_co2    =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
 ! receive
-sn_rcv_w10m   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-sn_rcv_taumod =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
-sn_rcv_tau    =       'oce only'             ,    'no'    , 'cartesian' , 'eastward-northward',  'U,V'
-sn_rcv_dqnsdt =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
-sn_rcv_qsr    =       'oce and ice'          ,    'no'    ,     ''      ,         ''          ,   ''
-sn_rcv_qns    =       'oce and ice'          ,    'no'    ,     ''      ,         ''          ,   ''
-sn_rcv_emp    =       'conservative'         ,    'no'    ,     ''      ,         ''          ,   ''
-sn_rcv_rnf    =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
-sn_rcv_cal    =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
-sn_rcv_co2    =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_w10m   =   'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_taumod =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_tau    =   'oce only'             ,    'no'    , 'cartesian' , 'eastward-northward',  'U,V'
+   sn_rcv_dqnsdt =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_qsr    =   'oce and ice'          ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_qns    =   'oce and ice'          ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_emp    =   'conservative'         ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_rnf    =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_cal    =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_co2    =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
+!
+   nn_cplmodel   =     1     !  Maximum number of models to/from which NEMO is potentialy sending/receiving data
+   ln_usecplmask = .false.   !  use a coupling mask file to merge data received from several models
+                             !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_isf
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_isf	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_isf	(revision 6040)
@@ -1,25 +1,26 @@
 !-----------------------------------------------------------------------
-&namsbc_isf    !  Top boundary layer (ISF)
+&namsbc_isf    !  Top boundary layer (ISF)                              (nn_isf >0)
 !-----------------------------------------------------------------------
-!              ! file name ! frequency (hours) ! variable ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation !
-!              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  !
+!              ! file name ! frequency (hours) ! variable ! time interp. ! clim   ! 'yearly'/ ! weights  ! rotation !
+!              !           !  (if <0  months)  !   name   !   (logical)  ! (T/F)  ! 'monthly' ! filename ! pairing  !
 ! nn_isf == 4
-   sn_qisf      = 'rnfisf' ,         -12      ,'sohflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
-   sn_fwfisf    = 'rnfisf' ,         -12      ,'sowflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_qisf      = 'rnfisf' ,         -12       ,'sohflisf',    .false.   , .true. , 'yearly'  ,  ''      ,   ''
+   sn_fwfisf    = 'rnfisf' ,         -12       ,'sowflisf',    .false.   , .true. , 'yearly'  ,  ''      ,   ''
 ! nn_isf == 3
-   sn_rnfisf    = 'runoffs' ,         -12      ,'sofwfisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_rnfisf    = 'runoffs',         -12       ,'sofwfisf',    .false.   , .true. , 'yearly'  ,  ''      ,   ''
 ! nn_isf == 2 and 3
-   sn_depmax_isf = 'runoffs' ,       -12        ,'sozisfmax' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
-   sn_depmin_isf = 'runoffs' ,       -12        ,'sozisfmin' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_depmax_isf = 'runoffs',        -12       ,'sozisfmax',   .false.   , .true. , 'yearly'  ,  ''      ,   ''
+   sn_depmin_isf = 'runoffs',        -12       ,'sozisfmin',   .false.   , .true. , 'yearly'  ,  ''      ,   ''
 ! nn_isf == 2
-   sn_Leff_isf = 'rnfisf' ,       0          ,'Leff'         ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_Leff_isf = 'rnfisf'  ,           0       ,'Leff'    ,    .false.   , .true. , 'yearly'  ,  ''      ,   ''
+
 ! for all case
-   ln_divisf   = .true.  ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?)
+   ln_divisf   = .true.   ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?)
 ! only for nn_isf = 1 or 2
-   rn_gammat0  = 1.0e-4   ! gammat coefficient used in blk formula
-   rn_gammas0  = 1.0e-4   ! gammas coefficient used in blk formula
+   rn_gammat0  = 1.e-4    ! gammat coefficient used in blk formula
+   rn_gammas0  = 1.e-4    ! gammas coefficient used in blk formula
 ! only for nn_isf = 1
    nn_isfblk   =  1       ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??)
-   rn_hisf_tbl =  30.      ! thickness of the top boundary layer           (Losh et al. 2008)
+   rn_hisf_tbl =  30.     ! thickness of the top boundary layer           (Losh et al. 2008)
                           ! 0 => thickness of the tbl = thickness of the first wet cell
    ln_conserve = .true.   ! conservative case (take into account meltwater advection)
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_mfs
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_mfs	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_mfs	(revision 6040)
@@ -2,13 +2,14 @@
 &namsbc_mfs   !   namsbc_mfs  MFS bulk formulae
 !-----------------------------------------------------------------------
-!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
-!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_wndi     =   'ecmwf'   ,        6          , 'u10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''    , ''
-   sn_wndj     =   'ecmwf'   ,        6          , 'v10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''    , ''
-   sn_clc      =   'ecmwf'   ,        6          , 'clc'     ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''    , ''
-   sn_msl      =   'ecmwf'   ,        6          , 'msl'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''    , ''
-   sn_tair     =   'ecmwf'   ,        6          , 't2'      ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''    , ''
-   sn_rhm      =   'ecmwf'   ,        6          , 'rh'      ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''    , ''
-   sn_prec     =   'precip'  ,        6          , 'precip'  ,    .true.    , .false. , 'daily'  ,'bicubic'    , ''    , ''
-   cn_dir      =   './'      !  root directory for the location of the bulk files
+!              !  file name  ! frequency (hours) ! variable ! time interp. !  clim  ! 'yearly'/ ! weights     ! rotation ! land/sea mask !
+!              !             !  (if <0  months)  !   name   !   (logical)  !  (T/F) ! 'monthly' ! filename    ! pairing  ! filename      !
+   sn_wndi     =   'ecmwf'   ,        6          , 'u10'    ,    .true.    , .false., 'daily'   ,'bicubic.nc' ,   ''     ,   ''
+   sn_wndj     =   'ecmwf'   ,        6          , 'v10'    ,    .true.    , .false., 'daily'   ,'bicubic.nc' ,   ''     ,   ''
+   sn_clc      =   'ecmwf'   ,        6          , 'clc'    ,    .true.    , .false., 'daily'   ,'bilinear.nc',   ''     ,   ''
+   sn_msl      =   'ecmwf'   ,        6          , 'msl'    ,    .true.    , .false., 'daily'   ,'bicubic.nc' ,   ''     ,   ''
+   sn_tair     =   'ecmwf'   ,        6          , 't2'     ,    .true.    , .false., 'daily'   ,'bicubic.nc' ,   ''     ,   ''
+   sn_rhm      =   'ecmwf'   ,        6          , 'rh'     ,    .true.    , .false., 'daily'   ,'bilinear.nc',   ''     ,   ''
+   sn_prec     =   'ecmwf'   ,        6          , 'precip' ,    .true.    , .true. , 'daily'   ,'bicubic.nc' ,   ''     ,   ''
+
+   cn_dir      = './ECMWF/'      !  root directory for the location of the bulk files
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_rnf
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_rnf	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_rnf	(revision 6040)
@@ -1,15 +1,14 @@
 !-----------------------------------------------------------------------
-&namsbc_rnf    !   runoffs namelist surface boundary condition
+&namsbc_rnf    !   runoffs namelist surface boundary condition          (ln_rnf=T)
 !-----------------------------------------------------------------------
-!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
-!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_rnf      = 'runoff_core_monthly',    -1    , 'sorunoff',   .true.     , .true. , 'yearly'  , ''       , ''       , ''
-   sn_cnf      = 'runoff_core_monthly',     0    , 'socoefr0',   .false.    , .true. , 'yearly'  , ''       , ''       , ''
-   sn_s_rnf    = 'runoffs'            ,    24    , 'rosaline',   .true.     , .true. , 'yearly'  , ''       , ''       , ''
-   sn_t_rnf    = 'runoffs'            ,    24    , 'rotemper',   .true.     , .true. , 'yearly'  , ''       , ''       , ''
-   sn_dep_rnf  = 'runoffs'            ,     0    , 'rodepth' ,   .false.    , .true. , 'yearly'  , ''       , ''       , ''
+!              !  file name           ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !                      !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_rnf      = 'runoff_core_monthly',        -1         , 'sorunoff',   .true.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_cnf      = 'runoff_core_monthly',         0         , 'socoefr0',   .false.    , .true. , 'yearly'  , ''       , ''       , ''
+   sn_s_rnf    = 'runoffs'            ,        24         , 'rosaline',   .true.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_t_rnf    = 'runoffs'            ,        24         , 'rotemper',   .true.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_dep_rnf  = 'runoffs'            ,         0         , 'rodepth' ,   .false.    , .true. , 'yearly'  , ''       , ''       , ''
 
    cn_dir       = './'      !  root directory for the location of the runoff files
-   ln_rnf_emp   = .false.   !  runoffs included into precipitation field (T) or into a file (F)
    ln_rnf_mouth = .true.    !  specific treatment at rivers mouths
    rn_hrnf      =  15.e0    !  depth over which enhanced vertical mixing is used
@@ -19,3 +18,7 @@
    ln_rnf_tem   = .false.   !  read in temperature information for runoff
    ln_rnf_sal   = .false.   !  read in salinity information for runoff
+   ln_rnf_depth_ini = .false.  ! compute depth at initialisation from runoff file
+   rn_rnf_max   = 5.735e-4  !  max value of the runoff climatologie over global domain ( ln_rnf_depth_ini = .true )
+   rn_dep_max   = 150.      !  depth over which runoffs is spread ( ln_rnf_depth_ini = .true )
+   nn_rnf_depth_file = 0    !  create (=1) a runoff depth file or not (=0)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_sas
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_sas	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_sas	(revision 6040)
@@ -4,11 +4,14 @@
 !              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
 !              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_usp      = 'sas_grid_U' ,    120    , 'vozocrtx' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
-   sn_vsp      = 'sas_grid_V' ,    120    , 'vomecrty' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
-   sn_tem      = 'sas_grid_T' ,    120    , 'sosstsst' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
-   sn_sal      = 'sas_grid_T' ,    120    , 'sosaline' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
-   sn_ssh      = 'sas_grid_T' ,    120    , 'sossheig' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
+   sn_usp      = 'sas_grid_U' ,    120           , 'vozocrtx' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
+   sn_vsp      = 'sas_grid_V' ,    120           , 'vomecrty' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
+   sn_tem      = 'sas_grid_T' ,    120           , 'sosstsst' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
+   sn_sal      = 'sas_grid_T' ,    120           , 'sosaline' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
+   sn_ssh      = 'sas_grid_T' ,    120           , 'sossheig' ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
+   sn_e3t      = 'sas_grid_T' ,    120           , 'e3t_m'    ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
+   sn_frq      = 'sas_grid_T' ,    120           , 'frq_m'    ,  .true.    , .true. ,   'yearly'  , ''       , ''             , ''
 
-   ln_3d_uv    = .true.    !  specify whether we are supplying a 3D u,v field
+   ln_3d_uve   = .true.    !  specify whether we are supplying a 3D u,v and e3 field
+   ln_read_frq = .false.    !  specify whether we must read frq or not
    cn_dir      = './'      !  root directory for the location of the bulk files are
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_ssr
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_ssr	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_ssr	(revision 6040)
@@ -1,4 +1,4 @@
 !-----------------------------------------------------------------------
-&namsbc_ssr    !   surface boundary condition : sea surface restoring
+&namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr=T)
 !-----------------------------------------------------------------------
 !              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_wave
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_wave	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namsbc_wave	(revision 6040)
@@ -1,12 +1,14 @@
 !-----------------------------------------------------------------------
-&namsbc_wave   ! External fields from wave model
+&namsbc_wave   ! External fields from wave model                        (ln_wave=T)
 !-----------------------------------------------------------------------
-!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
-!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' , .true.    , .false. , 'daily'  , ''       , ''       , ''
-   sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     , .true.   , .false. , 'daily'  ,''         , ''       , ''
-   sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     , .true.   , .false. , 'daily'  ,''         , ''       , ''
-   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   , .true.   , .false. , 'daily'  ,''         , ''       , ''
+!              !  file name  ! frequency (hours) ! variable    ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !             !  (if <0  months)  !   name      !  (logical)   !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff',     .true.   , .false., 'daily'   ,  ''      , ''       , ''
+   sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'    ,     .true.   , .false., 'daily'   ,  ''      , ''       , ''
+   sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'    ,     .true.   , .false., 'daily'   ,  ''      , ''       , ''
+   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'  ,     .true.   , .false., 'daily'   ,  ''      , ''       , ''
 !
-   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+   cn_dir_cdg  = './'      !  root directory for the location of drag coefficient files
+   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model
+   ln_sdw  = .false.       !  Computation of 3D stokes drift               
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_adv
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_adv	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_adv	(revision 6040)
@@ -1,11 +1,17 @@
 !-----------------------------------------------------------------------
-&namtra_adv    !   advection scheme for tracer
+&namtra_adv    !   advection scheme for tracer                          (default: NO advection)
 !-----------------------------------------------------------------------
-   ln_traadv_cen2   =  .false.  !  2nd order centered scheme
-   ln_traadv_tvd    =  .true.   !  TVD scheme
-   ln_traadv_muscl  =  .false.  !  MUSCL scheme
-   ln_traadv_muscl2 =  .false.  !  MUSCL2 scheme + cen2 at boundaries
-   ln_traadv_ubs    =  .false.  !  UBS scheme
-   ln_traadv_qck    =  .false.  !  QUICKEST scheme
-   ln_traadv_msc_ups=  .false.  !  use upstream scheme within muscl
+   ln_traadv_cen =  .false.  !  2nd order centered scheme
+      nn_cen_h   =  4               !  =2/4, horizontal 2nd order CEN / 4th order CEN
+      nn_cen_v   =  4               !  =2/4, vertical   2nd order CEN / 4th order COMPACT
+   ln_traadv_fct =  .false.  !  FCT scheme
+      nn_fct_h   =  2               !  =2/4, horizontal 2nd / 4th order 
+      nn_fct_v   =  2               !  =2/4, vertical   2nd / COMPACT 4th order 
+      nn_fct_zts =  0               !  >=1,  2nd order FCT scheme with vertical sub-timestepping
+      !                             !        (number of sub-timestep = nn_fct_zts)
+   ln_traadv_mus =  .false.  !  MUSCL scheme
+      ln_mus_ups =  .false.         !  use upstream scheme near river mouths
+   ln_traadv_ubs =  .false.  !  UBS scheme
+      nn_ubs_v   =  2               !  =2  , vertical 2nd order FCT
+   ln_traadv_qck =  .false.  !  QUICKEST scheme
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_dmp
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_dmp	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_dmp	(revision 6040)
@@ -1,10 +1,9 @@
 !-----------------------------------------------------------------------
-&namtra_dmp    !   tracer: T & S newtonian damping
+&namtra_dmp    !   tracer: T & S newtonian damping                      (default: NO)
 !-----------------------------------------------------------------------
-   ln_tradmp   =  .true.     !  add a damping termn (T) or not (F)
-   nn_zdmp     =    0        !  vertical   shape =0    damping throughout the water column
-                             !                   =1 no damping in the mixing layer (kz  criteria)
-                             !                   =2 no damping in the mixed  layer (rho crieria)
-   cn_resto    = 'resto.nc'  ! Name of file containing restoration coefficient field (use dmp_tools to create this)
-
+   ln_tradmp   =  .true.   !  add a damping termn (T) or not (F)
+   nn_zdmp     =    0      !  vertical   shape =0    damping throughout the water column
+                           !                   =1 no damping in the mixing layer (kz  criteria)
+                           !                   =2 no damping in the mixed  layer (rho crieria)
+   cn_resto    ='resto.nc' !  Name of file containing restoration coeff. field (use dmp_tools to create this)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_ldf
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_ldf	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_ldf	(revision 6040)
@@ -1,28 +1,32 @@
-!----------------------------------------------------------------------------------
-&namtra_ldf    !   lateral diffusion scheme for tracers
-!----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------
+&namtra_ldf    !   lateral diffusion scheme for tracers                 (default: NO diffusion)
+!-----------------------------------------------------------------------
    !                       !  Operator type:
-   ln_traldf_lap    =  .true.   !  laplacian operator
-   ln_traldf_bilap  =  .false.  !  bilaplacian operator
+   !                           !  no diffusion: set ln_traldf_lap=..._blp=F 
+   ln_traldf_lap   =  .false.  !    laplacian operator
+   ln_traldf_blp   =  .false.  !  bilaplacian operator
    !                       !  Direction of action:
-   ln_traldf_level  =  .false.  !  iso-level
-   ln_traldf_hor    =  .false.  !  horizontal (geopotential)   (needs "key_ldfslp" when ln_sco=T)
-   ln_traldf_iso    =  .true.   !  iso-neutral                 (needs "key_ldfslp")
-   !		       	   !  Griffies parameters              (all need "key_ldfslp")
-   ln_traldf_grif   =  .false.  !  use griffies triads
-   ln_traldf_gdia   =  .false.  !  output griffies eddy velocities
-   ln_triad_iso     =  .false.  !  pure lateral mixing in ML
-   ln_botmix_grif   =  .false.  !  lateral mixing on bottom
-   !                       !  Coefficients
-   ! Eddy-induced (GM) advection always used with Griffies; otherwise needs "key_traldf_eiv"
-   ! Value rn_aeiv_0 is ignored unless = 0 with Held-Larichev spatially varying aeiv
-   !                                  (key_traldf_c2d & key_traldf_eiv & key_orca_r2, _r1 or _r05)
-   rn_aeiv_0        =  2000.    !  eddy induced velocity coefficient [m2/s]
-   rn_aht_0         =  2000.    !  horizontal eddy diffusivity for tracers [m2/s]
-   rn_ahtb_0        =     0.    !  background eddy diffusivity for ldf_iso [m2/s]
-   !                                           (normally=0; not used with Griffies)
-   rn_slpmax        =     0.01  !  slope limit
-   rn_chsmag        =     1.    !  multiplicative factor in Smagorinsky diffusivity
-   rn_smsh          =     1.    !  Smagorinsky diffusivity: = 0 - use only sheer
-   rn_aht_m         =  2000.    !  upper limit or stability criteria for lateral eddy diffusivity (m2/s)
+   ln_traldf_lev   =  .false.  !  iso-level
+   ln_traldf_hor   =  .false.  !  horizontal (geopotential)
+   ln_traldf_iso   =  .false.  !  iso-neutral (standard operator)
+   ln_traldf_triad =  .false.  !  iso-neutral (triad    operator)
+   !
+   !		       	         !  iso-neutral options:        
+   ln_traldf_msc   =  .false.  !  Method of Stabilizing Correction (both operators)
+   rn_slpmax       =   0.01    !  slope limit                      (both operators)
+   ln_triad_iso    =  .false.  !  pure horizontal mixing in ML              (triad only)
+   rn_sw_triad     =  1        !  =1 switching triad ; =0 all 4 triads used (triad only)
+   ln_botmix_triad =  .false.  !  lateral mixing on bottom                  (triad only)
+   !
+   !                       !  Coefficients:
+   nn_aht_ijk_t    = 0         !  space/time variation of eddy coef
+   !                                !   =-20 (=-30)    read in eddy_diffusivity_2D.nc (..._3D.nc) file
+   !                                !   =  0           constant 
+   !                                !   = 10 F(k)      =ldf_c1d 
+   !                                !   = 20 F(i,j)    =ldf_c2d 
+   !                                !   = 21 F(i,j,t)  =Treguier et al. JPO 1997 formulation
+   !                                !   = 30 F(i,j,k)  =ldf_c2d * ldf_c1d
+   !                                !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing)
+   rn_aht_0        = 2000.     !  lateral eddy diffusivity   (lap. operator) [m2/s]
+   rn_bht_0        = 1.e+12    !  lateral eddy diffusivity (bilap. operator) [m4/s]
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_qsr
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_qsr	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtra_qsr	(revision 6040)
@@ -1,4 +1,4 @@
 !-----------------------------------------------------------------------
-&namtra_qsr    !   penetrative solar radiation
+&namtra_qsr    !   penetrative solar radiation                          (ln_traqsr=T)
 !-----------------------------------------------------------------------
 !              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
@@ -7,5 +7,4 @@
 
    cn_dir      = './'      !  root directory for the location of the runoff files
-   ln_traqsr   = .true.    !  Light penetration (T) or not (F)
    ln_qsr_rgb  = .true.    !  RGB (Red-Green-Blue) light penetration
    ln_qsr_2bd  = .false.   !  2 bands              light penetration
@@ -15,3 +14,4 @@
    rn_si0      =   0.35    !  RGB & 2 bands: shortess depth of extinction
    rn_si1      =   23.0    !  2 bands: longest depth of extinction
+   ln_qsr_ice  = .true.    !  light penetration for ice-model LIM3
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtrd
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtrd	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtrd	(revision 6040)
@@ -1,12 +1,14 @@
 !-----------------------------------------------------------------------
-&namtrd        !   diagnostics on dynamics and/or tracer trends         ("key_trddyn" and/or "key_trdtra")
-!              !       or mixed-layer trends or barotropic vorticity    ("key_trdmld" or     "key_trdvor")
+&namtrd        !   diagnostics on dynamics and/or tracer trends         (default F)
+!              !   and/or mixed-layer trends and/or barotropic vorticity
 !-----------------------------------------------------------------------
-   nn_trd      = 365       !  time step frequency dynamics and tracers trends
-   nn_ctls     =   0       !  control surface type in mixed-layer trends (0,1 or n<jpk)
-   rn_ucf      =   1.      !  unit conversion factor (=1 -> /seconds ; =86400. -> /day)
-   cn_trdrst_in      = "restart_mld"   ! suffix of ocean restart name (input)
-   cn_trdrst_out     = "restart_mld"   ! suffix of ocean restart name (output)
-   ln_trdmld_restart = .false.         !  restart for ML diagnostics
-   ln_trdmld_instant = .false.         !  flag to diagnose trends of instantantaneous or mean ML T/S
+   ln_glo_trd  = .false.   ! (T) global domain averaged diag for T, T^2, KE, and PE
+   ln_dyn_trd  = .false.   ! (T) 3D momentum trend output
+   ln_dyn_mxl  = .false.   ! (T) 2D momentum trends averaged over the mixed layer (not coded yet)
+   ln_vor_trd  = .false.   ! (T) 2D barotropic vorticity trends (not coded yet)
+   ln_KE_trd   = .false.   ! (T) 3D Kinetic   Energy     trends
+   ln_PE_trd   = .false.   ! (T) 3D Potential Energy     trends
+   ln_tra_trd  = .false.   ! (T) 3D tracer trend output
+   ln_tra_mxl  = .false.   ! (T) 2D tracer trends averaged over the mixed layer (not coded yet)
+   nn_trd      = 365       !  print frequency (ln_glo_trd=T) (unit=time step)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtsd
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtsd	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namtsd	(revision 6040)
@@ -2,8 +2,8 @@
 &namtsd    !   data : Temperature  & Salinity
 !-----------------------------------------------------------------------
-!          ! file name ! frequency (hours)         ! variable     ! time interp. ! clim  !'yearly' or ! weights  ! rotation ! land/sea mask !
-!          !           !  (if <0  months)          !   name       !  (logical)   ! (T/F) ! 'monthly'  ! filename ! pairing  ! filename      !
-   sn_tem  = 'data_1m_potential_temperature_nomask', -1,'votemper',  .true.      , .true., 'yearly'   , ' '      , ' '      , ''
-   sn_sal  = 'data_1m_salinity_nomask'             , -1,'vosaline',  .true.      , .true., 'yearly'   , ''       , ' '      , ''
+!          !  file name                            ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!          !                                       !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_tem  = 'data_1m_potential_temperature_nomask',         -1        ,'votemper' ,    .true.    , .true. , 'yearly'   , ''       ,   ''    ,    ''
+   sn_sal  = 'data_1m_salinity_nomask'             ,         -1        ,'vosaline' ,    .true.    , .true. , 'yearly'   , ''       ,   ''    ,    ''
    !
    cn_dir        = './'     !  root directory for the location of the runoff files
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzdf_gls
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzdf_gls	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzdf_gls	(revision 6040)
@@ -2,17 +2,17 @@
 &namzdf_gls                !   GLS vertical diffusion                   ("key_zdfgls")
 !-----------------------------------------------------------------------
-   rn_emin       = 1.e-6   !  minimum value of e   [m2/s2]
+   rn_emin       = 1.e-7   !  minimum value of e   [m2/s2]
    rn_epsmin     = 1.e-12  !  minimum value of eps [m2/s3]
    ln_length_lim = .true.  !  limit on the dissipation rate under stable stratification (Galperin et al., 1988)
-   rn_clim_galp  = 0.53    !  galperin limit
-   ln_crban      = .true.  !  Use Craig & Banner (1994) surface wave mixing parametrisation
+   rn_clim_galp  = 0.267   !  galperin limit
    ln_sigpsi     = .true.  !  Activate or not Burchard 2001 mods on psi schmidt number in the wb case
    rn_crban      = 100.    !  Craig and Banner 1994 constant for wb tke flux
    rn_charn      = 70000.  !  Charnock constant for wb induced roughness length
-   nn_tkebc_surf =   1     !  surface tke condition (0/1/2=Dir/Neum/Dir Mellor-Blumberg)
-   nn_tkebc_bot  =   1     !  bottom tke condition (0/1=Dir/Neum)
-   nn_psibc_surf =   1     !  surface psi condition (0/1/2=Dir/Neum/Dir Mellor-Blumberg)
-   nn_psibc_bot  =   1     !  bottom psi condition (0/1=Dir/Neum)
-   nn_stab_func  =   2     !  stability function (0=Galp, 1= KC94, 2=CanutoA, 3=CanutoB)
-   nn_clos       =   1     !  predefined closure type (0=MY82, 1=k-eps, 2=k-w, 3=Gen)
+   rn_hsro       =  0.02   !  Minimum surface roughness
+   rn_frac_hs    =   1.3   !  Fraction of wave height as roughness (if nn_z0_met=2)
+   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2)
+   nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum)
+   nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum)
+   nn_stab_func  =     2   !  stability function (0=Galp, 1= KC94, 2=CanutoA, 3=CanutoB)
+   nn_clos       =     1   !  predefined closure type (0=MY82, 1=k-eps, 2=k-w, 3=Gen)
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzdf_tke
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzdf_tke	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzdf_tke	(revision 6040)
@@ -21,5 +21,5 @@
                            !        = 1 add a tke source below the ML
                            !        = 2 add a tke source just at the base of the ML
-                           !        = 3 as = 1 applied on HF part of the stress    ("key_coupled")
+                           !        = 3 as = 1 applied on HF part of the stress    ("key_oasis3")
    rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
    nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzgr
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzgr	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzgr	(revision 6040)
@@ -1,8 +1,9 @@
 !-----------------------------------------------------------------------
-&namzgr        !   vertical coordinate
+&namzgr        !   vertical coordinate                                  (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_zco      = .false.   !  z-coordinate - full    steps   (T/F)      ("key_zco" may also be defined)
-   ln_zps      = .true.    !  z-coordinate - partial steps   (T/F)
-   ln_sco      = .false.   !  s- or hybrid z-s-coordinate    (T/F)
-   ln_isfcav   = .false.   !  ice shelf cavity               (T/F)
+   ln_zco      = .false.   !  z-coordinate - full    steps
+   ln_zps      = .false.   !  z-coordinate - partial steps
+   ln_sco      = .false.   !  s- or hybrid z-s-coordinate
+   ln_isfcav   = .false.   !  ice shelf cavity
+   ln_linssh   = .false.   !  linear free surface
 /
Index: /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzgr_sco
===================================================================
--- /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzgr_sco	(revision 6039)
+++ /branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Namelist/namzgr_sco	(revision 6040)
@@ -2,5 +2,5 @@
 &namzgr_sco    !   s-coordinate or hybrid z-s-coordinate
 !-----------------------------------------------------------------------
-   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+   ln_s_sh94   = .false.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
    ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
    ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
@@ -13,5 +13,5 @@
                         !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
    rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
-   rn_bb       =    0.8    !  stretching with SH94 s-sigma   
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma
                         !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
    rn_alpha    =    4.4    !  stretching with SF12 s-sigma
@@ -22,4 +22,4 @@
    rn_zb_b     =   -0.2    !  offset for calculating Zb
                         !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
-   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
+   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)
 /
