Index: unk/DOC/TexFiles/Chapters/Abstracts_Foreword - save.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Abstracts_Foreword - save.tex	(revision 998)
+++ 	(revision )
@@ -1,52 +1,0 @@
-
-% ================================================================
-% Abstract (English)
-% ================================================================
-
-\chapter*{Abstract}
-
-OPA is the ocean engine of NEMO, a framework of ocean related components developed for both research and operational purposes. It is a primitive equation model adapted to regional and global ocean circulation. It is intended to be a flexible tool for studying the ocean and its interactions with the others components of the earth climate system (atmosphere, sea-ice, biogeochemical tracers, ...) over a wide range of space and time scales. Prognostic variables are the three-dimensional velocity field, the 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 z full or partial step 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 a 1.5 turbulent closure  or a KPP scheme for the vertical mixing. Within NEMO, OPA is interfaced with a sea-ice model (LIM), passive tracer and biogeochemical models (TOP) and, via the OASIS coupler, with several atmospheric general circulation models. In addition, it can be run on many different computers, including shared and distributed memory multiprocessor computers.
-
-% ================================================================
-% Abstract (French)
-% ================================================================
-\chapter*{R\'{e}sum\'{e}}
-
-OPA est le moteur oc\'{e}anique de NEMO, un syst\`{e}me de mod\'{e}lisation de l'oc\'{e}an d\'{e}velopp\'{e} \`{a} des fins de recherche et d'oc\'{e}anographie op\'{e}rationnelle. Ce moteur est un mod\`{e}le aux \'{e}quations primitives de la circulation oc\'{e}anique r\'{e}gionale et globale. Il se veut un outil flexible pour \'{e}tudier sur un vaste spectre spatiotemporel l'oc\'{e}an et ses interactions avec les autres composantes du syst\`{e}me climatique terrestre (atmosph\`{e}re, glace de mer, traceurs biog\'{e}ochimiques...). Les variables pronostiques sont le champ tridimensionnel de vitesse, la hauteur de la mer , la temperature et la salinit\'{e}. 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 une fermeture turbulente d'ordre 1.5 ou un sch\'{e}ma KPP pour le m\'{e}lange vertical. Via l'infrastructure NEMO, OPA est interfac\'{e} avec un mod\`{e}le de glace de mer, des mod\`{e}les biog\'{e}ochimiques et de traceur passif, et, via le coupleur OASIS, \`{a} plusieurs mod\`{e}les de circulation g\'{e}n\'{e}rale atmosph\'{e}rique. En outre, il peut \^etre ex\'{e}cut\'{e} sur de nombreux calculateurs, y compris des calculateurs multiprocesseurs \`{a} m\'{e}moire partag\'{e}e ou distribu\'{e}e.
-
-
-% ================================================================
-% Disclaimer
-% ================================================================
-\chapter*{Disclaimer}
-
-OPA (an acronym for Ocean PArall\'{e}lis\'{e}) is the ocean component of NEMO (Nucleus for European Modelling of the Ocean (www.locean-ipsl.upmc.fr/NEMO). Like all components of NEMO, it is developed under the CECILL license, which is a french adaptation of the GNU GPL (General Public license). Anyone may use OPA freely for research purposes, and is encouraged to communicate back to the NEMO team its own developments and improvements. The model and the present document have been made available as a service to the community. We cannot certify that the code and its manual are free of errors. Bugs are inevitable and some have undoubtedly survived the testing phase. Users are encouraged to bring them to our attention. The author assumes no responsibility for problems, errors, or incorrect usage of OPA.
-
- \vspace{1cm}
-The OPA OGCM reference in papers and other publications is as follows:
- \vspace{0.5cm}
-
-Madec, G., 2007: NEMO ocean engine. \textit{Note du P\^ole de mod\'{e}lisation}, Institut Pierre-Simon Laplace (IPSL), France, NXX, \colorbox{yellow}{YYpp}.\\
-
-% ================================================================
-% Foreword
-% ================================================================
-\chapter*{Foreword}
-
-OPA is the ocean engine of the NEMO. This model, like all research tools, is in perpetual evolution. The present document describes the OPA model include in the release 2.3 of NEMO. This is the release 9 of OPA. This version differs significantly from version 8, documented in \citet{Madec1998}. The major modifications are :\\
-(1) transition to full native \textsc{Fortran} 90, deep code restructuring and drastic reduction of CPP keys, \\
-(2) introduction of partial step representation of bottom topography \\
-(3) partial reactivation of a terrain-following vertical coordinate ($s$- and hybrid $s$-$z$) with the addition of several options for pressure gradient computation \footnote{Partial support of $s$-coordinate: there is presently no support for neutral physics in $s$-coordinate and for the new options for horizontal pressure gradient computation with true equation of state.}, \\ 
-(4) more choices for the treatment of the free surface: full explicit, split-explicit , filtered and rigid-lid \\
-(5) non linear free surface option (variable level thickness distributed on the whole water column) \\
-(6) additional schemes for vector and flux forms of the momentum  advection \\
-(7) addition of several advection schemes on tracers \\
- (8) implementation of the AGRIF package (Adaptative Grid Refinement in \textsc{Fortran} ) \\
-(9) online diagnostics : tracers trend in the mixed layer and vorticity balance \\
-(10) rewriting of the I/O management \\
-(11) OASIS 3 and 4 couplers interfacing with atmospheric global circulation models. 
-
-In addition, several minor modifications in the coding have been introduced with the constant concern of improving performance on both scalar and vector computers. 
-
-At the time of this writing, the current release is NEMO 2.3. The new surface module described in this document is not yet part of the current distribution.
-
Index: /trunk/DOC/TexFiles/Chapters/Annex_C.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Annex_C.tex	(revision 998)
+++ /trunk/DOC/TexFiles/Chapters/Annex_C.tex	(revision 999)
@@ -610,5 +610,5 @@
 % ================================================================
 \section{Conservation Properties on Lateral Momentum Physics}
-\label{Apdx_C.3}
+\label{Apdx_dynldf_properties}
 
 
Index: /trunk/DOC/TexFiles/Chapters/Chap_LDF.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Chap_LDF.tex	(revision 998)
+++ /trunk/DOC/TexFiles/Chapters/Chap_LDF.tex	(revision 999)
@@ -9,14 +9,15 @@
 $\ $\newline    % force a new ligne
 
-The lateral physics on momentum and tracer equations have been given 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 the choices that underlie each lateral physics option. 
-Choosing one lateral physics 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 nam\_traldf and nam\_dynldf 
-below).
+The lateral physics terms in the momentum and tracer equations have been 
+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 nam\_traldf and nam\_dynldf below).
+
 %-----------------------------------nam_traldf - nam_dynldf--------------------------------------------
 \namdisplay{nam_traldf} 
@@ -28,35 +29,51 @@
 % Lateral Mixing Coefficients
 % ================================================================
-\section{Lateral Mixing Coefficient (\textbf{key\_ldftra\_c.d} - \textbf{key\_ldfdyn\_c.d)} }
+\section {Lateral Mixing Coefficient (\mdl{ldftra}, \mdl{ldfdyn)} }
 \label{LDF_coef}
 
 
 Introducing a space variation in the lateral eddy mixing coefficients changes 
-the model core memory requirement, adding up to four three-dimensional 
-arrays for geopotential or isopycnal second order operator applied to 
-momentum. Six cpp keys control the space variation of eddy 
-coefficients: three for momentum and three for tracer. They 
-allow to specify a space variation in the three space directions, in the 
-horizontal plane, or in the vertical only. The default option is a constant 
-value over the whole ocean on momentum and tracers. 
+the model core memory requirement, adding up to four extra three-dimensional 
+arrays for the geopotential or isopycnal second order operator applied to 
+momentum. Six CPP keys control the space variation of eddy coefficients: 
+three for momentum and three for tracer. The three choices allow: 
+a space variation in the three space directions, in the horizontal plane, 
+or in the vertical only. The default option is a constant value over the whole 
+ocean on both momentum and tracers. 
 
 The number of additional arrays that have to be defined and the gridpoint 
 position at which they are defined depend on both the space variation chosen 
-and the type of operator used. The resulting eddy viscosity and 
-diffusivity coefficients can be either single or multiple valued functions. 
-Changes in the computer code when switching from one option to another have 
-been minimized by introducing the eddy coefficients as statement function 
-(include file \hf{ldftra\_substitute} and \hf{ldfdyn\_substitute}). The functions are replaced by their actual meaning during the preprocessing step (cpp capability). The specification of the space variation of the coefficient is settled in \mdl{ldftra} and \mdl{ldfdyn}, or more precisely in include files \textit{ldftra\_cNd.h90} and \textit{ldfdyn\_cNd.h90}, with N=1, 2 or 3. The user can change these include files following his desiderata. The way the mixing coefficient are set in the reference version can be briefly described as follows:
+and the type of operator used. The resulting eddy viscosity and diffusivity 
+coefficients can be a function of more than one variable. Changes in the 
+computer code when switching from one option to another have been 
+minimized by introducing the eddy coefficients as statement functions
+(include file \hf{ldftra\_substitute} and \hf{ldfdyn\_substitute}). The functions 
+are replaced by their actual meaning during the preprocessing step (CPP). 
+The specification of the space variation of the coefficient is made in 
+\mdl{ldftra} and \mdl{ldfdyn}, or more precisely in include files 
+\textit{ldftra\_cNd.h90} and \textit{ldfdyn\_cNd.h90}, with N=1, 2 or 3. 
+The user can modify these include files as he/she wishes. The way the 
+mixing coefficient are set in the reference version can be briefly described 
+as follows:
 
 \subsubsection{Constant Mixing Coefficients (default option)}
-When none of the \textbf{key\_ldfdyn\_...} and \textbf{key\_ldftra\_...} keys are defined, a constant value over the whole ocean on momentum and tracers that is specified through 
-\np{ahm0} and \np{aht0} namelist parameters.
+When none of the \textbf{key\_ldfdyn\_...} and \textbf{key\_ldftra\_...} keys are 
+defined, a constant value is used over the whole ocean for momentum and 
+tracers, which is specified through the \np{ahm0} and \np{aht0} namelist 
+parameters.
 
 \subsubsection{Vertically varying Mixing Coefficients (\key{ldftra\_c1d} and \key{ldfdyn\_c1d})} 
-The 1D option is only available in $z$-coordinate with full step. Indeed in all the other type of vertical coordinate, the depth is a 3D function of (\textbf{i},\textbf{j},\textbf{j}) and therefore, introducing depth-dependant mixing coefficients will requires 3D arrays, $i.e.$ \key{ldftra\_c3d} and \key{ldftra\_c3d}.  In the 1D option, a hyperbolic variation of the lateral mixing coefficient is introduced in which the surface value is \np{aht0} (\np{ahm0}), the bottom value is 1/4 of the surface value, and the transition is round z=300~m with a width of 300~m ($i.e.$ both the depth and the width of the inflection point are set to 300~m). This profile is hard coded in \hf{ldftra\_c1d} file, but can be easily modified by users.
+The 1D option is only available when using the $z$-coordinate with full step. 
+Indeed in all the other types of vertical coordinate, the depth is a 3D function 
+of (\textbf{i},\textbf{j},\textbf{k}) and therefore, introducing depth-dependent 
+mixing coefficients will require 3D arrays, $i.e.$ \key{ldftra\_c3d} and \key{ldftra\_c3d}.  
+In the 1D option, a hyperbolic variation of the lateral mixing coefficient is introduced 
+in which the surface value is \np{aht0} (\np{ahm0}), the bottom value is 1/4 of 
+the surface value, and the transition takes place around z=300~m with a width 
+of 300~m ($i.e.$ both the depth and the width of the inflection point are set to 300~m). 
+This profile is hard coded in file \hf{ldftra\_c1d}, but can be easily modified by users.
 
 \subsubsection{Horizontally Varying Mixing Coefficients (\key{ldftra\_c2d} and \key{ldfdyn\_c2d})}
-
-By default the horizontal variation of the eddy coefficient depend on the local mesh size and the type of operator used:
+By default the horizontal variation of the eddy coefficient depends on the local mesh size and the type of operator used:
 \begin{equation} \label{Eq_title}
   A_l = \left\{     
@@ -67,37 +84,60 @@
 \quad \text{comments}
 \end{equation}
-where $e_{max}$ is the max of $e_1$ and $e_2$ taken over the whole masked ocean domain, and $A_o^l$ is \np{ahm0} (momentum) or \np{aht0} (tracer) namelist parameters. This variation is intended to reflect the lesser need for subgrid scale eddy mixing where the grid size is smaller in the domain. It was introduced in the context of the DYNAMO modelling project \citep{Willebrand2001}. 
+where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked ocean domain, and $A_o^l$ is the \np{ahm0} (momentum) or \np{aht0} (tracer) namelist parameter. This variation is intended to reflect the lesser need for subgrid scale eddy mixing where the grid size is smaller in the domain. It was introduced in the context of the DYNAMO modelling project \citep{Willebrand2001}. 
 %%%
 \gmcomment { not only that! stability reasons: with non uniform grid size, it is common to face a blow up of the model due to to large diffusive coefficient compare to the smallest grid size... this is especially true for bilaplacian (to be added in the text!)  }
 
-Other formulations can be introduced by the user for a given configuration. For example, in the ORCA2 global ocean model (\key{orca\_r2}), the laplacian viscous operator uses \np{ahm0}~=~$4.10^4 m^2.s^{-1}$ poleward of 20$^{\circ}$ north and south and decreases to \np{aht0}~=~$2.10^3 m^2.s^{-1}$ at the equator \citep{Madec1996, Delecluse_Madec_Bk00}. This specification can be found in \rou{ldf\_dyn\_c2d\_orca} routine defined in \mdl{ldfdyn\_c2d}. Similar specific horizontal variation can be found for Antarctic or Arctic sub-domain of ORCA2 and ORCA05 (\key{antarctic} or \key{arctic} defined, see \hf{ldfdyn\_antarctic} and \hf{ldfdyn\_arctic}).
+Other formulations can be introduced by the user for a given configuration. 
+For example, in the ORCA2 global ocean model (\key{orca\_r2}), the laplacian 
+viscosity operator uses \np{ahm0}~=~$4.10^4 m^2/s$ poleward of 20$^{\circ}$ 
+north and south and decreases linearly to \np{aht0}~=~$2.10^3 m^2/s$ 
+at the equator \citep{Madec1996, Delecluse_Madec_Bk00}. This modification 
+can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}. 
+Similar modified horizontal variations can be found with the Antarctic or Arctic 
+sub-domain options of ORCA2 and ORCA05 (\key{antarctic} or \key{arctic} 
+defined, see \hf{ldfdyn\_antarctic} and \hf{ldfdyn\_arctic}).
 
 \subsubsection{Space Varying Mixing Coefficients (\key{ldftra\_c3d} and \key{ldfdyn\_c3d})}
 
-The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases, $i.e.$ a hyperbolic tangent variation with depth associated with a grid size dependence of the magnitude of the coefficient. 
-
-\subsubsection{Space and time Varying Mixing Coefficients}
-
-There is no default specification of space and time varying mixing coefficient. The only case available is specific to ORCA2 and ORCA05 global ocean configurations (\key{orca\_r2} or \key{orca\_r05}). It provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and eddy induced velocity (ORCA05) that depends on the local growth rate of baroclinic instability. This specification is actually used when a ORCA key plus \key{traldf\_eiv} plus \key{traldf\_c2d} are defined.
+The 3D space variation of the mixing coefficient is simply the combination of the 
+1D and 2D cases, $i.e.$ a hyperbolic tangent variation with depth associated with 
+a grid size dependence of the magnitude of the coefficient. 
+
+\subsubsection{Space and Time Varying Mixing Coefficients}
+
+There is no default specification of space and time varying mixing coefficient. 
+The only case available is specific to the ORCA2 and ORCA05 global ocean 
+configurations (\key{orca\_r2} or \key{orca\_r05}). It provides only a tracer 
+mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and 
+eddy induced velocity (ORCA05) that depends on the local growth rate of 
+baroclinic instability. This specification is actually used when an ORCA key 
+and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined.
 
 A space variation in the eddy coefficient appeals several remarks:
 
-(1) the momentum diffusive operator acting along model level surfaces is 
+(1) the momentum diffusion operator acting along model level surfaces is 
 written in terms of curl and divergent components of the horizontal current 
 (see \S\ref{PE_ldf}). Although the eddy coefficient can be set to different values 
 in these two terms, this option is not available. 
 
-(2) with a horizontal varying viscosity, the quadratic integral constraints 
+(2) with an horizontally varying viscosity, the quadratic integral constraints 
 on enstrophy and on the square of the horizontal divergence for operators 
-acting along model-surfaces are no more satisfied (\colorbox{yellow}{Appendix C}).
+acting along model-surfaces are no longer satisfied 
+(Appendix~\ref{Apdx_dynldf_properties}).
 
 (3) for isopycnal diffusion on momentum or tracers, an additional purely 
 horizontal background diffusion with uniform coefficient can be added by 
-setting a non zero value of \np{ahmb0} or \np{ahtb0}, a background horizontal eddy 
-viscosity or diffusivity coefficient (\textbf{namelist parameters} which default value are $0$). Nevertheless, the technique used to compute the isopycnal slopes allows to get rid of such a background diffusion which introduces spurious diapycnal diffusion (see {\S\ref{LDF_slp}).
-
-(4) when an eddy induced advection is used (\key{trahdf\_eiv}), $A^{eiv}$ , the eddy induced coefficient has to be defined. Its space variations are controlled by the same CPP variable as for the eddy diffusivity coefficient (i.e. \textbf{key\_traldf\_cNd}). 
-
-(5) the eddy coefficient associated to a biharmonic operator must be set to a \emph{negative} value.
+setting a non zero value of \np{ahmb0} or \np{ahtb0}, a background horizontal 
+eddy viscosity or diffusivity coefficient (namelist parameters whose default 
+values are $0$). However, the technique used to compute the isopycnal 
+slopes is intended to get rid of such a background diffusion, since it introduces 
+spurious diapycnal diffusion (see {\S\ref{LDF_slp}).
+
+(4) when an eddy induced advection term is used (\key{trahdf\_eiv}), $A^{eiv}$, 
+the eddy induced coefficient has to be defined. Its space variations are controlled 
+by the same CPP variable as for the eddy diffusivity coefficient ($i.e.$ 
+\textbf{key\_traldf\_cNd}). 
+
+(5) the eddy coefficient associated with a biharmonic operator must be set to a \emph{negative} value.
 
 
@@ -112,12 +152,26 @@
 \gmcomment{  we should emphasize here that the implementation is a rather old one. Better work can be achieved by using \citet{Griffies1998, Griffies2004} 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$-coordinate, 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 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$. 
+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 $s$-coordinate}
-
-In $s$-coordinates, geopotential mixing ($i.e.$ horizontal one) $r_1$ and $r_2$ are the slopes between the geopotential and computational surfaces. Their discrete formulation is found by locally 
-vanishing the diffusive fluxes when $T$ is horizontally uniform, i.e. by replacing in \eqref{Eq_tra_ldf_iso} $T$ by $z_T$, the depth of $T$-point, and setting to zero the diffusive fluxes in the three directions. This leads to the following expression for the slopes:
+\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}
@@ -141,9 +195,16 @@
 
 %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}=T and \np{ln\_traldf\_hor}=T or \np{ln\_dynldf\_hor}=T. 
+
+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}
-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 to zero diffusive fluxes in the three directions leads to the following definition for the neutral slopes:
+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}
@@ -171,15 +232,32 @@
 %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 $z$-coordinate, the derivative of the  \eqref{Eq_ldfslp_iso} numerator is evaluated at a same depth ($T$-level which is also $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} have 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:
+%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. it 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 : ] the technique 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, there is no specific treatment for iso-neutral mixing in $s$-coordinate. In other word, iso-neutral mixing will only be accurately represented with a linear equation of state (\np{neos}=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 a wrong evaluation of the neutral slopes.
+\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, 
+there is no specific treatment for iso-neutral mixing in the $s$-coordinate. 
+In other words, iso-neutral mixing will only be accurately represented with a 
+linear equation of state (\np{neos}=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% 
@@ -188,4 +266,5 @@
 \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:
@@ -215,16 +294,32 @@
 \end{split}
 \end{equation}
-where $\alpha$ and $\beta$, the thermal expansion and saline contracion coefficients introduced in \S\ref{TRA_bn2}, have to be evaluated at the three velocity point. Inorder 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 $z$ and $zps$ cases.
-
+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 OGCM are modified \citep{Weaver1997, Griffies1998}. 
-Here, another strategy has been chosen \citep{Lazar1997}: 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 diffusive 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 diffusive operator that acts along large-scale (2~$\Delta$x) 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 a horizontal background mixing. 
+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{Weaver1997, 
+Griffies1998}. Here, another strategy has been chosen \citep{Lazar1997}: 
+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{Griffies1998} operator which has that property. 
@@ -237,5 +332,5 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-%There is three additional questions about the slope calculation. First the expression of the rotation tensor used have been obtain assuming the "small slope" approximation, so a bound has to be specified on slopes. Second, numerical stability issues also require a bound on slopes. Third, the question of boundary condition spefified on slopes...
+%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....
@@ -247,9 +342,24 @@
 flattening of isopycnals near the surface).
 
+For numerical stability reasons \citep{Cox1987, Griffies2004}, 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] \label{Fig_eiv_slp}  \begin{center}
 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_eiv_slp.pdf}
-\caption {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 and their have to adjust to the surface boundary (i.e. tend to zero at the surface as there is no mixing across the air-sea interface: wall boundary condition). Nevertheless, the profile between surface zero value and interior iso-neutral one is unknown, and especially the value at the based 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 actuelly used in \NEMO: 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 based of the mixed layer between  \textit{(b)}  and \textit{(c)}. 
-.}
+\caption {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}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -259,5 +369,12 @@
 \subsection{slopes for momentum iso-neutral mixing}
 
-The diffusive iso-neutral operator on momentum is the same as the on used on tracer but applied to each component of the velocity (see \eqref{Eq_dyn_ldf_iso} in section~\ref{DYN_ldf_iso}). The slopes between the surface along which the diffusive 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 $f-T$-, \textit{vw}-points for the $v$-component. They are computed as follows from the slopes used for tracer diffusion, i.e. \eqref{Eq_ldfslp_geo} and \eqref{Eq_ldfslp_iso} :
+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}
@@ -270,6 +387,6 @@
 \end{equation}
 
-The major issue remains in the specification of the boundary conditions. The 
-choice made consists in keeping the same boundary conditions as for lateral 
+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 
@@ -284,5 +401,9 @@
 \label{LDF_eiv}
 
-When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined), an eddy induced tracer advection term is added, the formulation of which depends on the slopes of iso-neutral surfaces. Contrary to iso-neutral mixing, the slopes use here are referenced to the geopotential surfaces, i.e. \eqref{Eq_ldfslp_geo} is used in $z$-coordinates, and the sum \eqref{Eq_ldfslp_geo}  + \eqref{Eq_ldfslp_iso} in $s$-coordinates. The eddy induced velocity is given by:
+When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined), 
+an eddy induced tracer advection term is added, the formulation of which 
+depends on the slopes of iso-neutral surfaces. Contrary to the case of iso-neutral 
+mixing, the slopes used here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo} is used in $z$-coordinates, and the sum \eqref{Eq_ldfslp_geo}  
++ \eqref{Eq_ldfslp_iso} in $s$-coordinates. The eddy induced velocity is given by: 
 \begin{equation} \label{Eq_ldfeiv}
 \begin{split}
@@ -292,12 +413,20 @@
 \end{split}
 \end{equation}
-where $A^{eiv}$ is the eddy induced velocity coefficient set through \np{aeiv}, a \textit{nam\_traldf} namelist parameter. 
-The three components of the eddy induced velocity are computed and add to the eulerian velocity in the mdl{traadv\_eiv}. This has been preferred to a separate computation of the advective trends associated to the eiv velocity as it allows to take advantage of all the advection schemes offered for the tracers (see \S\ref{TRA_adv}) and not only the $2^{nd}$ order advection scheme as in previous release of OPA \citep{Madec1998}. This is particularly useful for passive tracers where \emph{positivity}of the advection scheme is of paramount importance. 
-
-At surface, lateral and bottom boundaries, the eddy induced velocity and thus the 
-advective eddy fluxes of heat and salt are set to zero. 
-
-
-
-
-
+where $A^{eiv}$ is the eddy induced velocity coefficient whose value is set 
+through \np{aeiv}, a \textit{nam\_traldf} namelist parameter. 
+The three components of the eddy induced velocity are computed and add 
+to the eulerian velocity in \mdl{traadv\_eiv}. This has been preferred to a 
+separate computation of the advective trends associated with the eiv velocity, 
+since it allows us to take advantage of all the advection schemes offered for 
+the tracers (see \S\ref{TRA_adv}) and not just the $2^{nd}$ order advection 
+scheme as in previous releases of OPA \citep{Madec1998}. This is particularly 
+useful for passive tracers where \emph{positivity} of the advection scheme is 
+of paramount importance. 
+
+At the surface, lateral and bottom boundaries, the eddy induced velocity, 
+and thus the advective eddy fluxes of heat and salt, are set to zero. 
+
+
+
+
+
