Changeset 10354

Timestamp:

2018-11-21T17:59:55+01:00 (6 years ago)

Author:

nicolasmartin

Message:

Vast edition of LaTeX subfiles to improve the readability by cutting sentences in a more suitable way
Every sentence begins in a new line and if necessary is splitted around 110 characters lenght for side-by-side visualisation,
this setting may not be adequate for everyone but something has to be set.
The punctuation was the primer trigger for the cutting process, otherwise subordinators and coordinators, in order to mostly keep a meaning for each line

Location:

NEMO/trunk/doc/latex/NEMO/subfiles

Files:

• abstract_foreword.tex (modified) (2 diffs)
• annex_A.tex (modified) (25 diffs)
• annex_B.tex (modified) (15 diffs)
• annex_C.tex (modified) (45 diffs)
• annex_D.tex (modified) (7 diffs)
• annex_E.tex (modified) (36 diffs)
• annex_iso.tex (modified) (47 diffs)
• chap_ASM.tex (modified) (8 diffs)
• chap_CONFIG.tex (modified) (9 diffs)
• chap_DIA.tex (modified) (70 diffs)
• chap_DIU.tex (modified) (7 diffs)
• chap_DOM.tex (modified) (38 diffs)
• chap_DYN.tex (modified) (62 diffs)
• chap_LBC.tex (modified) (28 diffs)
• chap_LDF.tex (modified) (15 diffs)
• chap_OBS.tex (modified) (51 diffs)
• chap_SBC.tex (modified) (46 diffs)
• chap_STO.tex (modified) (8 diffs)
• chap_TRA.tex (modified) (50 diffs)
• chap_ZDF.tex (modified) (50 diffs)
• chap_conservation.tex (modified) (15 diffs)
• chap_misc.tex (modified) (10 diffs)
• chap_model_basics.tex (modified) (39 diffs)
• chap_model_basics_zstar.tex (modified) (15 diffs)
• chap_time_domain.tex (modified) (14 diffs)
• introduction.tex (modified) (3 diffs)

NEMO/trunk/doc/latex/NEMO/subfiles/abstract_foreword.tex

@@ -11 +11 @@
 
 {\small
-The ocean engine of NEMO (Nucleus for European Modelling of the Ocean) is a primitive 
-equation model adapted to regional and global ocean circulation problems. It is intended to 
-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 non-linear sea surface height, 
-the \textit{Conservative} Temperature and the \textit{Absolute} 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.
+  The ocean engine of NEMO (Nucleus for European Modelling of the Ocean) is a primitive equation model adapted to
+  regional and global ocean circulation problems.
+  It is intended to 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 non-linear sea surface height,
+  the \textit{Conservative} Temperature and the \textit{Absolute} 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 +32 @@
 \chapter*{Disclaimer}
 
-Like all components of NEMO, the ocean component is developed under the \href{http://www.cecill.info/}{CECILL license}, 
-which is a French adaptation of the GNU GPL (General Public License). Anyone may use it 
-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 NEMO.
+Like all components of NEMO,
+the ocean component is developed under the \href{http://www.cecill.info/}{CECILL license},
+which is a French adaptation of the GNU GPL (General Public License).
+Anyone may use it 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 NEMO.
 
 \vspace{1cm}

NEMO/trunk/doc/latex/NEMO/subfiles/annex_A.tex

@@ -19 +19 @@
 
 In order to establish the set of Primitive Equation in curvilinear $s$-coordinates
-($i.e.$ an orthogonal curvilinear coordinate in the horizontal and an Arbitrary Lagrangian 
-Eulerian (ALE) coordinate in the vertical), we start from the set of equations established 
-in \autoref{subsec:PE_zco_Eq} for the special case $k = z$ and thus $e_3 = 1$, and we introduce 
-an arbitrary vertical coordinate $a = a(i,j,z,t)$. Let us define a new vertical scale factor by 
-$e_3 = \partial z / \partial s$ (which now depends on $(i,j,z,t)$) and the horizontal 
-slope of $s-$surfaces by :
+($i.e.$ an orthogonal curvilinear coordinate in the horizontal and
+an Arbitrary Lagrangian Eulerian (ALE) coordinate in the vertical),
+we start from the set of equations established in \autoref{subsec:PE_zco_Eq} for
+the special case $k = z$ and thus $e_3 = 1$,
+and we introduce an arbitrary vertical coordinate $a = a(i,j,z,t)$.
+Let us define a new vertical scale factor by $e_3 = \partial z / \partial s$ (which now depends on $(i,j,z,t)$) and
+the horizontal slope of $s-$surfaces by:
 \begin{equation} \label{apdx:A_s_slope}
 \sigma _1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s 
@@ -31 +32 @@
 \end{equation}
 
-The chain rule to establish the model equations in the curvilinear $s-$coordinate 
-system is:
+The chain rule to establish the model equations in the curvilinear $s-$coordinate system is:
 \begin{equation} \label{apdx:A_s_chain_rule}
 \begin{aligned}
@@ -52 +52 @@
 \end{equation}
 
-In particular applying the time derivative chain rule to $z$ provides the expression 
-for $w_s$,  the vertical velocity of the $s-$surfaces referenced to a fix z-coordinate:
+In particular applying the time derivative chain rule to $z$ provides the expression for $w_s$,
+the vertical velocity of the $s-$surfaces referenced to a fix z-coordinate:
 \begin{equation} \label{apdx:A_w_in_s}
 w_s   =  \left.   \frac{\partial z }{\partial t}   \right|_s 
@@ -67 +67 @@
 \label{sec:A_continuity}
 
-Using (\autoref{apdx:A_s_chain_rule}) and the fact that the horizontal scale factors 
-$e_1$ and $e_2$ do not depend on the vertical coordinate, the divergence of 
-the velocity relative to the ($i$,$j$,$z$) coordinate system is transformed as follows
-in order to obtain its expression in the curvilinear $s-$coordinate system:
+Using (\autoref{apdx:A_s_chain_rule}) and
+the fact that the horizontal scale factors $e_1$ and $e_2$ do not depend on the vertical coordinate,
+the divergence of the velocity relative to the ($i$,$j$,$z$) coordinate system is transformed as follows in order to
+obtain its expression in the curvilinear $s-$coordinate system:
 
 \begin{subequations} 
@@ -128 +128 @@
 \end{subequations}
 
-Here, $w$ is the vertical velocity relative to the $z-$coordinate system. 
-Introducing the dia-surface velocity component, $\omega $, defined as 
-the volume flux across the moving $s$-surfaces per unit horizontal area:
+Here, $w$ is the vertical velocity relative to the $z-$coordinate system.
+Introducing the dia-surface velocity component,
+$\omega $, defined as the volume flux across the moving $s$-surfaces per unit horizontal area:
 \begin{equation} \label{apdx:A_w_s}
 \omega  = w - w_s - \sigma _1 \,u - \sigma _2 \,v    \\
 \end{equation}
-with $w_s$ given by \autoref{apdx:A_w_in_s}, we obtain the expression for 
-the divergence of the velocity in the curvilinear $s-$coordinate system:
+with $w_s$ given by \autoref{apdx:A_w_in_s},
+we obtain the expression for the divergence of the velocity in the curvilinear $s-$coordinate system:
 \begin{subequations} 
 \begin{align*} {\begin{array}{*{20}l} 
@@ -167 +167 @@
 \end{subequations}
 
-As a result, the continuity equation \autoref{eq:PE_continuity} in the 
-$s-$coordinates is:
+As a result, the continuity equation \autoref{eq:PE_continuity} in the $s-$coordinates is:
 \begin{equation} \label{apdx:A_sco_Continuity}
 \frac{1}{e_3 } \frac{\partial e_3}{\partial t} 
@@ -176 +175 @@
  +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0   
 \end{equation}
-A additional term has appeared that take into account the contribution of the time variation 
-of the vertical coordinate to the volume budget.
+A additional term has appeared that take into account
+the contribution of the time variation of the vertical coordinate to the volume budget.
 
 
@@ -186 +185 @@
 \label{sec:A_momentum}
 
-Here we only consider the first component of the momentum equation, 
+Here we only consider the first component of the momentum equation,
 the generalization to the second one being straightforward.
 
@@ -193 +192 @@
 $\bullet$ \textbf{Total derivative in vector invariant form}
 
-Let us consider \autoref{eq:PE_dyn_vect}, the first component of the momentum 
-equation in the vector invariant form. Its total $z-$coordinate time derivative, 
-$\left. \frac{D u}{D t} \right|_z$ can be transformed as follows in order to obtain 
+Let us consider \autoref{eq:PE_dyn_vect}, the first component of the momentum equation in the vector invariant form.
+Its total $z-$coordinate time derivative,
+$\left. \frac{D u}{D t} \right|_z$ can be transformed as follows in order to obtain
 its expression in the curvilinear $s-$coordinate system:
 
@@ -258 +257 @@
 \end{subequations}
 %
-Applying the time derivative chain rule (first equation of (\autoref{apdx:A_s_chain_rule}))
-to $u$ and using (\autoref{apdx:A_w_in_s}) provides the expression of the last term 
-of the right hand side,
+Applying the time derivative chain rule (first equation of (\autoref{apdx:A_s_chain_rule})) to $u$ and
+using (\autoref{apdx:A_w_in_s}) provides the expression of the last term of the right hand side,
 \begin{equation*} {\begin{array}{*{20}l} 
 w_s  \;\frac{\partial u}{\partial s} 
@@ -267 +265 @@
 \end{array} }     
 \end{equation*}
-leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative, 
+leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative,
 $i.e.$ the total $s-$coordinate time derivative :
 \begin{align} \label{apdx:A_sco_Dt_vect}
@@ -276 +274 @@
   + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s}   
 \end{align}
-Therefore, the vector invariant form of the total time derivative has exactly the same 
-mathematical form in $z-$ and $s-$coordinates. This is not the case for the flux form
-as shown in next paragraph.
+Therefore, the vector invariant form of the total time derivative has exactly the same mathematical form in
+$z-$ and $s-$coordinates.
+This is not the case for the flux form as shown in next paragraph.
 
 $\ $\newline    % force a new ligne
@@ -284 +282 @@
 $\bullet$ \textbf{Total derivative in flux form}
 
-Let us start from the total time derivative in the curvilinear $s-$coordinate system 
-we have just establish. Following the procedure used to establish (\autoref{eq:PE_flux_form}), 
-it can be transformed into :
+Let us start from the total time derivative in the curvilinear $s-$coordinate system we have just establish.
+Following the procedure used to establish (\autoref{eq:PE_flux_form}), it can be transformed into :
 %\begin{subequations} 
 \begin{align*} {\begin{array}{*{20}l} 
@@ -355 +352 @@
 \end{subequations}
 which leads to the $s-$coordinate flux formulation of the total $s-$coordinate time derivative, 
-$i.e.$ the total $s-$coordinate time derivative in flux form :
+$i.e.$ the total $s-$coordinate time derivative in flux form:
 \begin{flalign}\label{apdx:A_sco_Dt_flux}
 \left. \frac{D u}{D t} \right|_s   = \frac{1}{e_3}  \left. \frac{\partial ( e_3\,u)}{\partial t} \right|_s  
@@ -363 +360 @@
 \end{flalign}
 which is the total time derivative expressed in the curvilinear $s-$coordinate system.
-It has the same form as in the $z-$coordinate but for the vertical scale factor 
-that has appeared inside the time derivative which comes from the modification 
-of (\autoref{apdx:A_sco_Continuity}), the continuity equation.
+It has the same form as in the $z-$coordinate but for
+the vertical scale factor that has appeared inside the time derivative which
+comes from the modification of (\autoref{apdx:A_sco_Continuity}),
+the continuity equation.
 
 $\ $\newline    % force a new ligne
@@ -380 +378 @@
 \end{split}
 \end{equation*}
-Applying similar manipulation to the second component and replacing 
-$\sigma _1$ and $\sigma _2$ by their expression \autoref{apdx:A_s_slope}, it comes:
+Applying similar manipulation to the second component and
+replacing $\sigma _1$ and $\sigma _2$ by their expression \autoref{apdx:A_s_slope}, it comes:
 \begin{equation} \label{apdx:A_grad_p_1}
 \begin{split}
@@ -394 +392 @@
 \end{equation}
 
-An additional term appears in (\autoref{apdx:A_grad_p_1}) which accounts for the 
-tilt of $s-$surfaces with respect to geopotential $z-$surfaces.
-
-As in $z$-coordinate, the horizontal pressure gradient can be split in two parts
-following \citet{Marsaleix_al_OM08}. Let defined a density anomaly, $d$, by $d=(\rho - \rho_o)/ \rho_o$,
-and a hydrostatic pressure anomaly, $p_h'$, by $p_h'= g \; \int_z^\eta d \; e_3 \; dk$. 
+An additional term appears in (\autoref{apdx:A_grad_p_1}) which accounts for
+the tilt of $s-$surfaces with respect to geopotential $z-$surfaces.
+
+As in $z$-coordinate,
+the horizontal pressure gradient can be split in two parts following \citet{Marsaleix_al_OM08}.
+Let defined a density anomaly, $d$, by $d=(\rho - \rho_o)/ \rho_o$,
+and a hydrostatic pressure anomaly, $p_h'$, by $p_h'= g \; \int_z^\eta d \; e_3 \; dk$.
 The pressure is then given by:
 \begin{equation*} 
@@ -416 +415 @@
 \end{equation*}
 
-Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and using the definition of 
-the density anomaly it comes the expression in two parts:
+Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and
+using the definition of the density anomaly it comes the expression in two parts:
 \begin{equation} \label{apdx:A_grad_p_2}
 \begin{split}
@@ -429 +428 @@
 \end{split}
 \end{equation}
-This formulation of the pressure gradient is characterised by the appearance of a term depending on the 
-the sea surface height only (last term on the right hand side of expression \autoref{apdx:A_grad_p_2}).
-This term will be loosely termed \textit{surface pressure gradient}
-whereas the first term will be termed the 
-\textit{hydrostatic pressure gradient} by analogy to the $z$-coordinate formulation. 
-In fact, the the true surface pressure gradient is $1/\rho_o \nabla (\rho \eta)$, and 
-$\eta$ is implicitly included in the computation of $p_h'$ through the upper bound of 
-the vertical integration.
- 
+This formulation of the pressure gradient is characterised by the appearance of
+a term depending on the sea surface height only
+(last term on the right hand side of expression \autoref{apdx:A_grad_p_2}).
+This term will be loosely termed \textit{surface pressure gradient} whereas
+the first term will be termed the \textit{hydrostatic pressure gradient} by analogy to
+the $z$-coordinate formulation.
+In fact, the true surface pressure gradient is $1/\rho_o \nabla (\rho \eta)$,
+and $\eta$ is implicitly included in the computation of $p_h'$ through the upper bound of the vertical integration.
+
 
 $\ $\newline    % force a new ligne
@@ -443 +442 @@
 $\bullet$ \textbf{The other terms of the momentum equation}
 
-The coriolis and forcing terms as well as the the vertical physics remain unchanged 
-as they involve neither time nor space derivatives. The form of the lateral physics is 
-discussed in \autoref{apdx:B}.
+The coriolis and forcing terms as well as the the vertical physics remain unchanged as
+they involve neither time nor space derivatives.
+The form of the lateral physics is discussed in \autoref{apdx:B}.
 
 
@@ -452 +451 @@
 $\bullet$ \textbf{Full momentum equation}
 
-To sum up, in a curvilinear $s$-coordinate system, the vector invariant momentum equation 
-solved by the model has the same mathematical expression as the one in a curvilinear 
-$z-$coordinate, except for the pressure gradient term :
+To sum up, in a curvilinear $s$-coordinate system,
+the vector invariant momentum equation solved by the model has the same mathematical expression as
+the one in a curvilinear $z-$coordinate, except for the pressure gradient term:
 \begin{subequations} \label{apdx:A_dyn_vect}
 \begin{multline} \label{apdx:A_PE_dyn_vect_u}
@@ -475 +474 @@
 \end{multline}
 \end{subequations}
-whereas the flux form momentum equation differ from it by the formulation of both
-the time derivative and the pressure gradient term  :
+whereas the flux form momentum equation differs from it by
+the formulation of both the time derivative and the pressure gradient term:
 \begin{subequations} \label{apdx:A_dyn_flux}
 \begin{multline} \label{apdx:A_PE_dyn_flux_u}
@@ -503 +502 @@
 \end{equation}
 
-It is important to realize that the change in coordinate system has only concerned
-the position on the vertical. It has not affected (\textbf{i},\textbf{j},\textbf{k}), the 
-orthogonal curvilinear set of unit vectors. ($u$,$v$) are always horizontal velocities
-so that their evolution is driven by \emph{horizontal} forces, in particular 
-the pressure gradient. By contrast, $\omega$ is not $w$, the third component of the velocity,
-but the dia-surface velocity component, $i.e.$ the volume flux across the moving 
-$s$-surfaces per unit horizontal area. 
+It is important to realize that the change in coordinate system has only concerned the position on the vertical.
+It has not affected (\textbf{i},\textbf{j},\textbf{k}), the orthogonal curvilinear set of unit vectors.
+($u$,$v$) are always horizontal velocities so that their evolution is driven by \emph{horizontal} forces,
+in particular the pressure gradient.
+By contrast, $\omega$ is not $w$, the third component of the velocity, but the dia-surface velocity component,
+$i.e.$ the volume flux across the moving $s$-surfaces per unit horizontal area. 
 
 
@@ -518 +516 @@
 \label{sec:A_tracer}
 
-The tracer equation is obtained using the same calculation as for the continuity 
-equation and then regrouping the time derivative terms in the left hand side :
+The tracer equation is obtained using the same calculation as for the continuity equation and then
+regrouping the time derivative terms in the left hand side :
 
 \begin{multline} \label{apdx:A_tracer}
@@ -531 +529 @@
 
 
-The expression for the advection term is a straight consequence of (A.4), the 
-expression of the 3D divergence in the $s-$coordinates established above. 
+The expression for the advection term is a straight consequence of (A.4),
+the expression of the 3D divergence in the $s-$coordinates established above. 
 
 \end{document}

NEMO/trunk/doc/latex/NEMO/subfiles/annex_B.tex

@@ -19 +19 @@
 
 \subsubsection*{In z-coordinates}
-In $z$-coordinates, the horizontal/vertical second order tracer diffusion operator
-is given by:
+In $z$-coordinates, the horizontal/vertical second order tracer diffusion operator is given by:
 \begin{eqnarray} \label{apdx:B1}
  &D^T = \frac{1}{e_1 \, e_2}      \left[
@@ -30 +29 @@
 
 \subsubsection*{In generalized vertical coordinates}
-In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and
-$\sigma_2$ by \autoref{apdx:A_s_slope} and the vertical/horizontal ratio of diffusion
-coefficient by $\epsilon = A^{vT} / A^{lT}$. The diffusion operator is given by:
+In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and $\sigma_2$ by \autoref{apdx:A_s_slope} and
+the vertical/horizontal ratio of diffusion coefficient by $\epsilon = A^{vT} / A^{lT}$.
+The diffusion operator is given by:
 
 \begin{equation} \label{apdx:B2}
@@ -56 +55 @@
 \end{subequations}
 
-Equation \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any
-additional assumption. Indeed, for the special case $k=z$ and thus $e_3 =1$,
-we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A}
-and use \autoref{apdx:A_s_slope} and \autoref{apdx:A_s_chain_rule}.
-Since no cross horizontal derivative $\partial _i \partial _j $ appears in
-\autoref{apdx:B1}, the ($i$,$z$) and ($j$,$z$) planes are independent.
-The derivation can then be demonstrated for the ($i$,$z$)~$\to$~($j$,$s$)
-transformation without any loss of generality:
+Equation \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption.
+Indeed, for the special case $k=z$ and thus $e_3 =1$,
+we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A} and
+use \autoref{apdx:A_s_slope} and \autoref{apdx:A_s_chain_rule}.
+Since no cross horizontal derivative $\partial _i \partial _j $ appears in \autoref{apdx:B1},
+the ($i$,$z$) and ($j$,$z$) planes are independent.
+The derivation can then be demonstrated for the ($i$,$z$)~$\to$~($j$,$s$) transformation without
+any loss of generality:
 
 \begin{subequations}
@@ -143 +142 @@
 \subsubsection*{In z-coordinates}
 
-The iso/diapycnal diffusive tensor $\textbf {A}_{\textbf I}$ expressed in the ($i$,$j$,$k$)
-curvilinear coordinate system in which the equations of the ocean circulation model are
-formulated, takes the following form \citep{Redi_JPO82}:
+The iso/diapycnal diffusive tensor $\textbf {A}_{\textbf I}$ expressed in
+the ($i$,$j$,$k$) curvilinear coordinate system in which
+the equations of the ocean circulation model are formulated,
+takes the following form \citep{Redi_JPO82}:
 
 \begin{equation} \label{apdx:B3}
@@ -155 +155 @@
 \end{array} }} \right]
 \end{equation}
-where ($a_1$, $a_2$) are the isopycnal slopes in ($\textbf{i}$,
-$\textbf{j}$) directions, relative to geopotentials:
+where ($a_1$, $a_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, relative to geopotentials:
 \begin{equation*}
 a_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1}
@@ -164 +163 @@
 \end{equation*}
 
-In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean, so
-$\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}:
+In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean,
+so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}:
 \begin{subequations} \label{apdx:B4}
 \begin{equation} \label{apdx:B4a}
@@ -183 +182 @@
 
 
-Physically, the full tensor \autoref{apdx:B3}
-represents strong isoneutral diffusion on a plane parallel to the isoneutral
-surface and weak dianeutral diffusion perpendicular to this plane.
-However, the approximate `weak-slope' tensor \autoref{apdx:B4a} represents strong
-diffusion along the isoneutral surface, with weak
-\emph{vertical}  diffusion -- the principal axes of the tensor are no
-longer orthogonal. This simplification also decouples
-the ($i$,$z$) and ($j$,$z$) planes of the tensor. The weak-slope operator therefore takes the same
-form, \autoref{apdx:B4}, as \autoref{apdx:B2}, the diffusion operator for geopotential
-diffusion written in non-orthogonal $i,j,s$-coordinates. Written out
-explicitly,
+Physically, the full tensor \autoref{apdx:B3} represents strong isoneutral diffusion on a plane parallel to
+the isoneutral surface and weak dianeutral diffusion perpendicular to this plane.
+However,
+the approximate `weak-slope' tensor \autoref{apdx:B4a} represents strong diffusion along the isoneutral surface,
+with weak \emph{vertical} diffusion -- the principal axes of the tensor are no longer orthogonal.
+This simplification also decouples the ($i$,$z$) and ($j$,$z$) planes of the tensor.
+The weak-slope operator therefore takes the same form, \autoref{apdx:B4}, as \autoref{apdx:B2},
+the diffusion operator for geopotential diffusion written in non-orthogonal $i,j,s$-coordinates.
+Written out explicitly,
 
 \begin{multline} \label{apdx:B_ldfiso}
@@ -204 +201 @@
 
 The isopycnal diffusion operator \autoref{apdx:B4},
-\autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its
-square. The demonstration of the first property is trivial as \autoref{apdx:B4} is the divergence
-of fluxes. Let us demonstrate the second one:
+\autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its square.
+The demonstration of the first property is trivial as \autoref{apdx:B4} is the divergence of fluxes.
+Let us demonstrate the second one:
 \begin{equation*}
 \iiint\limits_D T\;\nabla .\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv
@@ -229 +226 @@
 \end{subequations}
 \addtocounter{equation}{-1}
- the property becomes obvious.
+the property becomes obvious.
 
 \subsubsection*{In generalized vertical coordinates}
 
-Because the weak-slope operator \autoref{apdx:B4}, \autoref{apdx:B_ldfiso} is decoupled
-in the ($i$,$z$) and ($j$,$z$) planes, it may be transformed into
-generalized $s$-coordinates in the same way as \autoref{sec:B_1} was transformed into
-\autoref{sec:B_2}. The resulting operator then takes the simple form
+Because the weak-slope operator \autoref{apdx:B4},
+\autoref{apdx:B_ldfiso} is decoupled in the ($i$,$z$) and ($j$,$z$) planes,
+it may be transformed into generalized $s$-coordinates in the same way as
+\autoref{sec:B_1} was transformed into \autoref{sec:B_2}.
+The resulting operator then takes the simple form
 
 \begin{equation} \label{apdx:B_ldfiso_s}
@@ -249 +247 @@
 \end{equation}
 
-where ($r_1$, $r_2$) are the isopycnal slopes in ($\textbf{i}$,
-$\textbf{j}$) directions, relative to $s$-coordinate surfaces:
+where ($r_1$, $r_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions,
+relative to $s$-coordinate surfaces:
 \begin{equation*}
 r_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial s}} \right)^{-1}
@@ -258 +256 @@
 \end{equation*}
 
-To prove  \autoref{apdx:B5}  by direct re-expression of \autoref{apdx:B_ldfiso} is
-straightforward, but laborious. An easier way is first to note (by reversing the
-derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that the
-weak-slope operator may be \emph{exactly} reexpressed in 
-non-orthogonal $i,j,\rho$-coordinates as
+To prove \autoref{apdx:B5} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious.
+An easier way is first to note (by reversing the derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that
+the weak-slope operator may be \emph{exactly} reexpressed in non-orthogonal $i,j,\rho$-coordinates as
 
 \begin{equation} \label{apdx:B5}
@@ -273 +269 @@
 \end{array} }} \right).
 \end{equation}
-Then direct transformation from $i,j,\rho$-coordinates to
-$i,j,s$-coordinates gives \autoref{apdx:B_ldfiso_s} immediately.
-
-Note that the weak-slope approximation is only made in
-transforming from the (rotated,orthogonal) isoneutral axes to the
-non-orthogonal $i,j,\rho$-coordinates. The further transformation
-into $i,j,s$-coordinates is exact, whatever the steepness of
-the  $s$-surfaces, in the same way as the transformation of
-horizontal/vertical Laplacian diffusion in $z$-coordinates,
+Then direct transformation from $i,j,\rho$-coordinates to $i,j,s$-coordinates gives
+\autoref{apdx:B_ldfiso_s} immediately.
+
+Note that the weak-slope approximation is only made in transforming from
+the (rotated,orthogonal) isoneutral axes to the non-orthogonal $i,j,\rho$-coordinates.
+The further transformation into $i,j,s$-coordinates is exact, whatever the steepness of the $s$-surfaces,
+in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates,
 \autoref{sec:B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces.
 
@@ -291 +285 @@
 \label{sec:B_3}
 
-The second order momentum diffusion operator (Laplacian) in the $z$-coordinate
-is found by applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian
-of a vector,  to the horizontal velocity vector :
+The second order momentum diffusion operator (Laplacian) in the $z$-coordinate is found by
+applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian of a vector,
+to the horizontal velocity vector:
 \begin{align*}
 \Delta {\textbf{U}}_h
@@ -329 +323 @@
 \end{array} }} \right)
 \end{align*}
-Using \autoref{eq:PE_div}, the definition of the horizontal divergence, the third
-componant of the second vector is obviously zero and thus :
+Using \autoref{eq:PE_div}, the definition of the horizontal divergence,
+the third componant of the second vector is obviously zero and thus :
 \begin{equation*}
 \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right)
 \end{equation*}
 
-Note that this operator ensures a full separation between the vorticity and horizontal
-divergence fields (see \autoref{apdx:C}). It is only equal to a Laplacian
-applied to each component in Cartesian coordinates, not on the sphere.
-
-The horizontal/vertical second order (Laplacian type) operator used to diffuse
-horizontal momentum in the $z$-coordinate therefore takes the following form :
+Note that this operator ensures a full separation between
+the vorticity and horizontal divergence fields (see \autoref{apdx:C}).
+It is only equal to a Laplacian applied to each component in Cartesian coordinates, not on the sphere.
+
+The horizontal/vertical second order (Laplacian type) operator used to diffuse horizontal momentum in
+the $z$-coordinate therefore takes the following form:
 \begin{equation} \label{apdx:B_Lap_U}
  {\textbf{D}}^{\textbf{U}} =
@@ -360 +354 @@
 \end{align*}
 
-Note Bene: introducing a rotation in \autoref{apdx:B_Lap_U} does not lead to a
-useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate.
-Similarly, we did not found an expression of practical use for the geopotential
-horizontal/vertical Laplacian operator in the $s$-coordinate. Generally,
-\autoref{apdx:B_Lap_U} is used in both $z$- and $s$-coordinate systems, that is
-a Laplacian diffusion is applied on momentum along the coordinate directions.
+Note Bene: introducing a rotation in \autoref{apdx:B_Lap_U} does not lead to
+a useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate.
+Similarly, we did not found an expression of practical use for
+the geopotential horizontal/vertical Laplacian operator in the $s$-coordinate.
+Generally, \autoref{apdx:B_Lap_U} is used in both $z$- and $s$-coordinate systems,
+that is a Laplacian diffusion is applied on momentum along the coordinate directions.
 \end{document}

NEMO/trunk/doc/latex/NEMO/subfiles/annex_C.tex

@@ -22 +22 @@
 \label{sec:C.0}
 
-Notation used in this appendix in the demonstations :
+Notation used in this appendix in the demonstations:
 
 fluxes at the faces of a $T$-box:
@@ -37 +37 @@
 
 $dv=e_1\,e_2\,e_3 \,di\,dj\,dk$  is the volume element, with only $e_3$ that depends on time.
-$D$ and $S$ are the ocean domain volume and surface, respectively. 
-No wetting/drying is allow ($i.e.$ $\frac{\partial S}{\partial t} = 0$) 
-Let $k_s$ and $k_b$ be the ocean surface and bottom, resp. 
+$D$ and $S$ are the ocean domain volume and surface, respectively.
+No wetting/drying is allow ($i.e.$ $\frac{\partial S}{\partial t} = 0$).
+Let $k_s$ and $k_b$ be the ocean surface and bottom, resp.
 ($i.e.$ $s(k_s) = \eta$ and $s(k_b)=-H$, where $H$ is the bottom depth).
 \begin{flalign*}
@@ -60 +60 @@
                                                        =  \int_D { \frac{1}{e_3} \partial_t \left( e_3 \, Q \right) dv } =0
 \end{equation*}
-equation of evolution of $Q$ written as the time evolution of the vertical content of $Q$
-like for tracers, or momentum in flux form, the quadratic quantity $\frac{1}{2}Q^2$ is conserved when :
+equation of evolution of $Q$ written as
+the time evolution of the vertical content of $Q$ like for tracers, or momentum in flux form,
+the quadratic quantity $\frac{1}{2}Q^2$ is conserved when:
 \begin{flalign*}
 \partial_t \left(   \int_D{ \frac{1}{2} \,Q^2\;dv }   \right)
@@ -74 +75 @@
   -   \frac{1}{2} \int_D {  \frac{Q^2}{e_3} \partial_t (e_3) \;dv }
 \end{flalign}
-equation of evolution of $Q$ written as the time evolution of $Q$
-like for momentum in vector invariant form, the quadratic quantity $\frac{1}{2}Q^2$ is conserved when :
+equation of evolution of $Q$ written as the time evolution of $Q$ like for momentum in vector invariant form,
+the quadratic quantity $\frac{1}{2}Q^2$ is conserved when:
 \begin{flalign*}
 \partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right)
@@ -82 +83 @@
 +  \int_D { \frac{1}{2} Q^2 \, \partial_t e_3  \;e_1e_2\;di\,dj\,dk } \\
 \end{flalign*}
-that is in a more compact form :
+that is in a more compact form:
 \begin{flalign} \label{eq:Q2_vect}
 \partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right)
@@ -97 +98 @@
 
 
-The discretization of pimitive equation in $s$-coordinate ($i.e.$ time and space varying 
-vertical coordinate) must be chosen so that the discrete equation of the model satisfy 
-integral constrains on energy and enstrophy. 
+The discretization of pimitive equation in $s$-coordinate ($i.e.$ time and space varying vertical coordinate)
+must be chosen so that the discrete equation of the model satisfy integral constrains on energy and enstrophy. 
 
 
 Let us first establish those constraint in the continuous world.
-The total energy ($i.e.$ kinetic plus potential energies) is conserved :
+The total energy ($i.e.$ kinetic plus potential energies) is conserved:
 \begin{flalign} \label{eq:Tot_Energy}
   \partial_t \left( \int_D \left( \frac{1}{2} {\textbf{U}_h}^2 +  \rho \, g \, z \right) \;dv \right)  = & 0
 \end{flalign}
-under the following assumptions: no dissipation, no forcing 
-(wind, buoyancy flux, atmospheric pressure variations), mass 
-conservation, and closed domain. 
-
-This equation can be transformed to obtain several sub-equalities. 
-The transformation for the advection term depends on whether 
-the vector invariant form or the flux form is used for the momentum equation.
-Using \autoref{eq:Q2_vect} and introducing \autoref{apdx:A_dyn_vect} in \autoref{eq:Tot_Energy} 
-for the former form and
-Using \autoref{eq:Q2_flux} and introducing \autoref{apdx:A_dyn_flux} in \autoref{eq:Tot_Energy} 
-for the latter form  leads to:
+under the following assumptions: no dissipation, no forcing (wind, buoyancy flux, atmospheric pressure variations),
+mass conservation, and closed domain. 
+
+This equation can be transformed to obtain several sub-equalities.
+The transformation for the advection term depends on whether the vector invariant form or
+the flux form is used for the momentum equation.
+Using \autoref{eq:Q2_vect} and introducing \autoref{apdx:A_dyn_vect} in
+\autoref{eq:Tot_Energy} for the former form and
+using \autoref{eq:Q2_flux} and introducing \autoref{apdx:A_dyn_flux} in
+\autoref{eq:Tot_Energy} for the latter form leads to:
 
 \begin{subequations} \label{eq:E_tot}
@@ -348 +347 @@
 
 Substituting the discrete expression of the time derivative of the velocity either in vector invariant,
-leads to the discrete equivalent of the four equations \autoref{eq:E_tot_flux}. 
+leads to the discrete equivalent of the four equations \autoref{eq:E_tot_flux}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -356 +355 @@
 \label{subsec:C_vor}
 
-Let $q$, located at $f$-points, be either the relative ($q=\zeta / e_{3f}$), or  
-the planetary ($q=f/e_{3f}$), or the total potential vorticity ($q=(\zeta +f) /e_{3f}$). 
-Two discretisation of the vorticity term (ENE and EEN) allows the conservation of 
-the kinetic energy.
+Let $q$, located at $f$-points, be either the relative ($q=\zeta / e_{3f}$),
+or the planetary ($q=f/e_{3f}$), or the total potential vorticity ($q=(\zeta +f) /e_{3f}$).
+Two discretisation of the vorticity term (ENE and EEN) allows the conservation of the kinetic energy.
 % -------------------------------------------------------------------------------------------------------------
 %       Vorticity Term with ENE scheme
@@ -366 +364 @@
 \label{subsec:C_vorENE} 
 
-For the ENE scheme, the two components of the vorticity term are given by :
+For the ENE scheme, the two components of the vorticity term are given by:
 \begin{equation*}
 - e_3 \, q \;{\textbf{k}}\times {\textbf {U}}_h    \equiv 
@@ -377 +375 @@
 \end{equation*}
 
-This formulation does not conserve the enstrophy but it does conserve the 
-total kinetic energy. Indeed, the kinetic energy tendency associated to the 
-vorticity term and averaged over the ocean domain can be transformed as
-follows:
+This formulation does not conserve the enstrophy but it does conserve the total kinetic energy.
+Indeed, the kinetic energy tendency associated to the vorticity term and
+averaged over the ocean domain can be transformed as follows:
 \begin{flalign*}
 &\int\limits_D -  \left(  e_3 \, q \;\textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv &&  \\
@@ -412 +409 @@
 \end{aligned}   } \right.
 \end{equation} 
-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$,
+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: 
 \begin{equation} \tag{\ref{eq:Q_triads}}
@@ -420 +416 @@
 \end{equation}
 
-This formulation does conserve the total kinetic energy. Indeed,
+This formulation does conserve the total kinetic energy.
+Indeed,
 \begin{flalign*}
 &\int\limits_D - \textbf{U}_h \cdot   \left(  \zeta \;\textbf{k} \times \textbf{U}_h  \right)  \;  dv &&  \\
@@ -473 +470 @@
 \label{subsec:C_zad} 
 
-The change of Kinetic Energy (KE) due to the vertical advection is exactly 
-balanced by the change of KE due to the horizontal gradient of KE~:
+The change of Kinetic Energy (KE) due to the vertical advection is exactly balanced by the change of KE due to the horizontal gradient of KE~:
 \begin{equation*}
     \int_D \textbf{U}_h \cdot \frac{1}{e_3 } \omega \partial_k \textbf{U}_h \;dv
@@ -480 +476 @@
    +   \frac{1}{2} \int_D {  \frac{{\textbf{U}_h}^2}{e_3} \partial_t ( e_3) \;dv }  \\
 \end{equation*}
-Indeed, using successively \autoref{eq:DOM_di_adj} ($i.e.$ the skew symmetry 
-property of the $\delta$ operator) and the continuity equation, then 
-\autoref{eq:DOM_di_adj} again, then the commutativity of operators 
-$\overline {\,\cdot \,}$ and $\delta$, and finally \autoref{eq:DOM_mi_adj} 
-($i.e.$ the symmetry property of the $\overline {\,\cdot \,}$ operator) 
+Indeed, using successively \autoref{eq:DOM_di_adj} ($i.e.$ the skew symmetry property of the $\delta$ operator)
+and the continuity equation, then \autoref{eq:DOM_di_adj} again,
+then the commutativity of operators $\overline {\,\cdot \,}$ and $\delta$, and finally \autoref{eq:DOM_mi_adj}
+($i.e.$ the symmetry property of the $\overline {\,\cdot \,}$ operator)
 applied in the horizontal and vertical directions, it becomes:
 \begin{flalign*}
@@ -543 +538 @@
 \end{flalign*}
 
-There is two main points here. First, the satisfaction of this property links the choice of 
-the discrete formulation of the vertical advection and of the horizontal gradient 
-of KE. Choosing one imposes the other. For example KE can also be discretized 
-as $1/2\,({\overline u^{\,i}}^2 + {\overline v^{\,j}}^2)$. This leads to the following 
-expression for the vertical advection:
+There is two main points here.
+First, the satisfaction of this property links the choice of the discrete formulation of the vertical advection and
+of the horizontal gradient of KE.
+Choosing one imposes the other.
+For example KE can also be discretized as $1/2\,({\overline u^{\,i}}^2 + {\overline v^{\,j}}^2)$.
+This leads to the following expression for the vertical advection:
 \begin{equation*}
 \frac{1} {e_3 }\; \omega\; \partial_k \textbf{U}_h
@@ -557 +553 @@
 \end{array}} } \right)
 \end{equation*}
-a formulation that requires an additional horizontal mean in contrast with 
-the one used in NEMO. Nine velocity points have to be used instead of 3. 
+a formulation that requires an additional horizontal mean in contrast with the one used in NEMO.
+Nine velocity points have to be used instead of 3.
 This is the reason why it has not been chosen.
 
-Second, as soon as the chosen $s$-coordinate depends on time, an extra constraint
-arises on the time derivative of the volume at $u$- and $v$-points:
+Second, as soon as the chosen $s$-coordinate depends on time,
+an extra constraint arises on the time derivative of the volume at $u$- and $v$-points:
 \begin{flalign*}
 e_{1u}\,e_{2u}\,\partial_t (e_{3u}) =\overline{ e_{1t}\,e_{2t}\;\partial_t (e_{3t}) }^{\,i+1/2}    \\
@@ -583 +579 @@
 
 \gmcomment{
-A pressure gradient has no contribution to the evolution of the vorticity as the 
-curl of a gradient is zero. In the $z$-coordinate, this property is satisfied locally 
-on a C-grid with 2nd order finite differences (property \autoref{eq:DOM_curl_grad}). 
+  A pressure gradient has no contribution to the evolution of the vorticity as the curl of a gradient is zero.
+  In the $z$-coordinate, this property is satisfied locally on a C-grid with 2nd order finite differences
+  (property \autoref{eq:DOM_curl_grad}). 
 }
 
-When the equation of state is linear ($i.e.$ when an advection-diffusion equation 
-for density can be derived from those of temperature and salinity) the change of 
-KE due to the work of pressure forces is balanced by the change of potential 
-energy due to buoyancy forces: 
+When the equation of state is linear
+($i.e.$ when an advection-diffusion equation for density can be derived from those of temperature and salinity)
+the change of KE due to the work of pressure forces is balanced by
+the change of potential energy due to buoyancy forces: 
 \begin{equation*}
 - \int_D  \left. \nabla p \right|_z \cdot \textbf{U}_h \;dv 
@@ -598 +594 @@
 \end{equation*}
 
-This property can be satisfied in a discrete sense for both $z$- and $s$-coordinates. 
-Indeed, defining the depth of a $T$-point, $z_t$, as the sum of the vertical scale 
-factors at $w$-points starting from the surface, the work of pressure forces can be 
-written as:
+This property can be satisfied in a discrete sense for both $z$- and $s$-coordinates.
+Indeed, defining the depth of a $T$-point, $z_t$,
+as the sum of the vertical scale factors at $w$-points starting from the surface,
+the work of pressure forces can be written as:
 \begin{flalign*}
 &- \int_D  \left. \nabla p \right|_z \cdot \textbf{U}_h \;dv   
@@ -658 +654 @@
 \end{flalign*}
 The first term is exactly the first term of the right-hand-side of \autoref{eq:KE+PE_vect_discrete}.
-It remains to demonstrate that the last term, which is obviously a discrete analogue of 
-$\int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to the last term of \autoref{eq:KE+PE_vect_discrete}.
+It remains to demonstrate that the last term,
+which is obviously a discrete analogue of $\int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to
+the last term of \autoref{eq:KE+PE_vect_discrete}.
 In other words, the following property must be satisfied:
 \begin{flalign*}
@@ -666 +663 @@
 \end{flalign*}
 
-Let introduce $p_w$ the pressure at $w$-point such that $\delta_k [p_w] = - \rho \,g\,e_{3t}$. 
+Let introduce $p_w$ the pressure at $w$-point such that $\delta_k [p_w] = - \rho \,g\,e_{3t}$.
 The right-hand-side of the above equation can be transformed as follows:
 
@@ -718 +715 @@
 
 
-Note that this property strongly constrains the discrete expression of both 
-the depth of $T-$points and of the term added to the pressure gradient in the
-$s$-coordinate. Nevertheless, it is almost never satisfied since a linear equation 
-of state is rarely used.
+Note that this property strongly constrains the discrete expression of both the depth of $T-$points and
+of the term added to the pressure gradient in the $s$-coordinate.
+Nevertheless, it is almost never satisfied since a linear equation of state is rarely used.
 
 
@@ -755 +751 @@
 \end{flalign*}
 
-Substituting the discrete expression of the time derivative of the velocity either in vector invariant or in flux form,
-leads to the discrete equivalent of the 
+Substituting the discrete expression of the time derivative of the velocity either in
+vector invariant or in flux form, leads to the discrete equivalent of the ????
 
 
@@ -771 +767 @@
 \label{subsec:C.3.3} 
 
-In flux from the vorticity term reduces to a Coriolis term in which the Coriolis 
-parameter has been modified to account for the ``metric'' term. This altered 
-Coriolis parameter is discretised at an f-point. It is given by: 
+In flux from the vorticity term reduces to a Coriolis term in which
+the Coriolis parameter has been modified to account for the ``metric'' term.
+This altered Coriolis parameter is discretised at an f-point.
+It is given by:
 \begin{equation*}
 f+\frac{1} {e_1 e_2 } \left( v \frac{\partial e_2 } {\partial i} - u \frac{\partial e_1 } {\partial j}\right)\;
@@ -781 +778 @@
 \end{equation*}
 
-Either the ENE or EEN scheme is then applied to obtain the vorticity term in flux form. 
-It therefore conserves the total KE. The derivation is the same as for the 
-vorticity term in the vector invariant form (\autoref{subsec:C_vor}).
+Either the ENE or EEN scheme is then applied to obtain the vorticity term in flux form.
+It therefore conserves the total KE.
+The derivation is the same as for the vorticity term in the vector invariant form (\autoref{subsec:C_vor}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -791 +788 @@
 \label{subsec:C.3.4} 
 
-The flux form operator of the momentum advection is evaluated using a 
-centered second order finite difference scheme. Because of the flux form, 
-the discrete operator does not contribute to the global budget of linear 
-momentum. Because of the centered second order scheme, it conserves 
-the horizontal kinetic energy, that is :
+The flux form operator of the momentum advection is evaluated using
+a centered second order finite difference scheme.
+Because of the flux form, the discrete operator does not contribute to the global budget of linear momentum.
+Because of the centered second order scheme, it conserves the horizontal kinetic energy, that is:
 
 \begin{equation} \label{eq:C_ADV_KE_flux}
@@ -804 +800 @@
 \end{equation}
 
-Let us first consider the first term of the scalar product ($i.e.$ just the the terms 
-associated with the i-component of the advection) :
+Let us first consider the first term of the scalar product
+($i.e.$ just the the terms associated with the i-component of the advection):
 \begin{flalign*}
 &  - \int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv   \\
@@ -845 +841 @@
    \biggl\{     \left(   \frac{1}{e_{3t}} \frac{\partial e_{3t}}{\partial t}   \right) \; b_t     \biggr\}    &&& \\
 \end{flalign*}
-Applying similar manipulation applied to the second term of the scalar product 
-leads to :
+Applying similar manipulation applied to the second term of the scalar product leads to:
 \begin{equation*}
  -  \int_D \textbf{U}_h \cdot     \left(                 {{\begin{array} {*{20}c}
@@ -854 +849 @@
    \biggl\{     \left(   \frac{1}{e_{3t}} \frac{\partial e_{3t}}{\partial t}   \right) \; b_t     \biggr\}    
 \end{equation*}
-which is the discrete form of 
-$ \frac{1}{2} \int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv $. 
+which is the discrete form of $ \frac{1}{2} \int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv $.
 \autoref{eq:C_ADV_KE_flux} is thus satisfied.
 
 
-When the UBS scheme is used to evaluate the flux form momentum advection, 
-the discrete operator does not contribute to the global budget of linear momentum 
-(flux form). The horizontal kinetic energy is not conserved, but forced to decay 
-($i.e.$ the scheme is diffusive). 
+When the UBS scheme is used to evaluate the flux form momentum advection,
+the discrete operator does not contribute to the global budget of linear momentum (flux form).
+The horizontal kinetic energy is not conserved, but forced to decay ($i.e.$ the scheme is diffusive). 
 
 
@@ -894 +887 @@
 \end{equation} 
 
-The scheme does not allow but the conservation of the total kinetic energy but the conservation 
-of $q^2$, the potential enstrophy for a horizontally non-divergent flow ($i.e.$ when $\chi$=$0$). 
-Indeed, using the symmetry or skew symmetry properties of the operators ( \autoref{eq:DOM_mi_adj} 
-and \autoref{eq:DOM_di_adj}), it can be shown that:
+The scheme does not allow but the conservation of the total kinetic energy but the conservation of $q^2$,
+the potential enstrophy for a horizontally non-divergent flow ($i.e.$ when $\chi$=$0$).
+Indeed, using the symmetry or skew symmetry properties of the operators
+( \autoref{eq:DOM_mi_adj} and \autoref{eq:DOM_di_adj}),
+it can be shown that:
 \begin{equation} \label{eq:C_1.1}
 \int_D {q\,\;{\textbf{k}}\cdot \frac{1} {e_3} \nabla \times \left( {e_3 \, q \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} \equiv 0
 \end{equation}
-where $dv=e_1\,e_2\,e_3 \; di\,dj\,dk$ is the volume element. Indeed, using 
-\autoref{eq:dynvor_ens}, the discrete form of the right hand side of \autoref{eq:C_1.1} 
-can be transformed as follow:
+where $dv=e_1\,e_2\,e_3 \; di\,dj\,dk$ is the volume element.
+Indeed, using \autoref{eq:dynvor_ens},
+the discrete form of the right hand side of \autoref{eq:C_1.1} can be transformed as follow:
 \begin{flalign*} 
 &\int_D q \,\; \textbf{k} \cdot \frac{1} {e_3 } \nabla \times 
@@ -955 +949 @@
 \end{aligned}   } \right.
 \end{equation} 
-where the indices $i_p$ and $k_p$ take the following value: 
+where the indices $i_p$ and $k_p$ take the following values: 
 $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: 
@@ -966 +960 @@
 This formulation does conserve the potential enstrophy for a horizontally non-divergent flow ($i.e.$ $\chi=0$). 
 
-Let consider one of the vorticity triad, for example ${^{i}_j}\mathbb{Q}^{+1/2}_{+1/2} $, 
-similar manipulation can be done for the 3 others. The discrete form of the right hand 
-side of \autoref{eq:C_1.1} applied to this triad only can be transformed as follow:
+Let consider one of the vorticity triad, for example ${^{i}_j}\mathbb{Q}^{+1/2}_{+1/2} $,
+similar manipulation can be done for the 3 others.
+The discrete form of the right hand side of \autoref{eq:C_1.1} applied to
+this triad only can be transformed as follow:
 
 \begin{flalign*} 
@@ -1020 +1015 @@
 
 
-All the numerical schemes used in NEMO are written such that the tracer content 
-is conserved by the internal dynamics and physics (equations in flux form). 
-For advection, only the CEN2 scheme ($i.e.$ $2^{nd}$ order finite different scheme) 
-conserves the global variance of tracer. Nevertheless the other schemes ensure 
-that the global variance decreases ($i.e.$ they are at least slightly diffusive). 
-For diffusion, all the schemes ensure the decrease of the total tracer variance, 
-except the iso-neutral operator. There is generally no strict conservation of mass, 
-as the equation of state is non linear with respect to $T$ and $S$. In practice, 
-the mass is conserved to a very high accuracy. 
+All the numerical schemes used in NEMO are written such that the tracer content is conserved by
+the internal dynamics and physics (equations in flux form).
+For advection,
+only the CEN2 scheme ($i.e.$ $2^{nd}$ order finite different scheme) conserves the global variance of tracer.
+Nevertheless the other schemes ensure that the global variance decreases
+($i.e.$ they are at least slightly diffusive).
+For diffusion, all the schemes ensure the decrease of the total tracer variance, except the iso-neutral operator.
+There is generally no strict conservation of mass,
+as the equation of state is non linear with respect to $T$ and $S$.
+In practice, the mass is conserved to a very high accuracy. 
 % -------------------------------------------------------------------------------------------------------------
 %       Advection Term
@@ -1049 +1045 @@
 
 
-Whatever the advection scheme considered it conserves of the tracer content as all 
-the scheme are written in flux form. Indeed,  let $T$ be the tracer and $\tau_u$, $\tau_v$, 
-and $\tau_w$ its interpolated values at velocity point (whatever the interpolation is), 
+Whatever the advection scheme considered it conserves of the tracer content as
+all the scheme are written in flux form.
+Indeed, let $T$ be the tracer and its $\tau_u$, $\tau_v$, and $\tau_w$ interpolated values at velocity point
+(whatever the interpolation is),
 the conservation of the tracer content due to the advection tendency is obtained as follows: 
 \begin{flalign*}
@@ -1067 +1064 @@
 \end{flalign*}
 
-The conservation of the variance of tracer due to the advection tendency 
-can be achieved only with the CEN2 scheme, $i.e.$ when 
-$\tau_u= \overline T^{\,i+1/2}$, $\tau_v= \overline T^{\,j+1/2}$, and $\tau_w= \overline T^{\,k+1/2}$. 
+The conservation of the variance of tracer due to the advection tendency can be achieved only with the CEN2 scheme,
+$i.e.$ when $\tau_u= \overline T^{\,i+1/2}$, $\tau_v= \overline T^{\,j+1/2}$, and $\tau_w= \overline T^{\,k+1/2}$. 
 It can be demonstarted as follows:
 \begin{flalign*}
@@ -1103 +1099 @@
 
 
-The discrete formulation of the horizontal diffusion of momentum ensures the 
-conservation of potential vorticity and the horizontal divergence, and the 
-dissipation of the square of these quantities ($i.e.$ enstrophy and the 
-variance of the horizontal divergence) as well as the dissipation of the 
-horizontal kinetic energy. In particular, when the eddy coefficients are 
-horizontally uniform, it ensures a complete separation of vorticity and 
-horizontal divergence fields, so that diffusion (dissipation) of vorticity 
-(enstrophy) does not generate horizontal divergence (variance of the 
-horizontal divergence) and \textit{vice versa}. 
-
-These properties of the horizontal diffusion operator are a direct consequence 
-of properties \autoref{eq:DOM_curl_grad} and \autoref{eq:DOM_div_curl}. 
-When the vertical curl of the horizontal diffusion of momentum (discrete sense) 
-is taken, the term associated with the horizontal gradient of the divergence is 
-locally zero. 
+The discrete formulation of the horizontal diffusion of momentum ensures
+the conservation of potential vorticity and the horizontal divergence,
+and the dissipation of the square of these quantities
+($i.e.$ enstrophy and the variance of the horizontal divergence) as well as
+the dissipation of the horizontal kinetic energy.
+In particular, when the eddy coefficients are horizontally uniform,
+it ensures a complete separation of vorticity and horizontal divergence fields,
+so that diffusion (dissipation) of vorticity (enstrophy) does not generate horizontal divergence
+(variance of the horizontal divergence) and \textit{vice versa}. 
+
+These properties of the horizontal diffusion operator are a direct consequence of
+properties \autoref{eq:DOM_curl_grad} and \autoref{eq:DOM_div_curl}.
+When the vertical curl of the horizontal diffusion of momentum (discrete sense) is taken,
+the term associated with the horizontal gradient of the divergence is locally zero. 
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1125 +1120 @@
 \label{subsec:C.6.1}
 
-The lateral momentum diffusion term conserves the potential vorticity :
+The lateral momentum diffusion term conserves the potential vorticity:
 \begin{flalign*}
 &\int \limits_D \frac{1} {e_3 } \textbf{k} \cdot \nabla \times 
@@ -1211 +1206 @@
 \label{subsec:C.6.3}
 
-The lateral momentum diffusion term dissipates the enstrophy when the eddy 
-coefficients are horizontally uniform:
+The lateral momentum diffusion term dissipates the enstrophy when the eddy coefficients are horizontally uniform:
 \begin{flalign*}
 &\int\limits_D  \zeta \; \textbf{k} \cdot \nabla \times 
@@ -1236 +1230 @@
 \label{subsec:C.6.4}
 
-When the horizontal divergence of the horizontal diffusion of momentum 
-(discrete sense) is taken, the term associated with the vertical curl of the 
-vorticity is zero locally, due to \autoref{eq:DOM_div_curl}. 
-The resulting term conserves the $\chi$ and dissipates $\chi^2$ 
-when the eddy coefficients are horizontally uniform.
+When the horizontal divergence of the horizontal diffusion of momentum (discrete sense) is taken,
+the term associated with the vertical curl of the vorticity is zero locally, due to \autoref{eq:DOM_div_curl}. 
+The resulting term conserves the $\chi$ and dissipates $\chi^2$ when the eddy coefficients are horizontally uniform.
 \begin{flalign*}
 & \int\limits_D  \nabla_h \cdot 
@@ -1291 +1283 @@
 \label{sec:C.7}
 
-As for the lateral momentum physics, the continuous form of the vertical diffusion 
-of momentum satisfies several integral constraints. The first two are associated 
-with the conservation of momentum and the dissipation of horizontal kinetic energy:
+As for the lateral momentum physics,
+the continuous form of the vertical diffusion of momentum satisfies several integral constraints.
+The first two are associated with the conservation of momentum and the dissipation of horizontal kinetic energy:
 \begin{align*}
  \int\limits_D   \frac{1} {e_3 }\; \frac{\partial } {\partial k} 
@@ -1306 +1298 @@
 \end{align*}
 
-The first property is obvious. The second results from:
+The first property is obvious.
+The second results from:
 \begin{flalign*}
 \int\limits_D 
@@ -1326 +1319 @@
 \end{flalign*}
 
-The vorticity is also conserved. Indeed:
+The vorticity is also conserved.
+Indeed:
 \begin{flalign*}
 \int \limits_D 
@@ -1346 +1340 @@
 \end{flalign*}
 
-If the vertical diffusion coefficient is uniform over the whole domain, the 
-enstrophy is dissipated, $i.e.$
+If the vertical diffusion coefficient is uniform over the whole domain, the enstrophy is dissipated, $i.e.$
 \begin{flalign*}
 \int\limits_D \zeta \, \textbf{k} \cdot \nabla \times 
@@ -1378 +1371 @@
       &\left[  \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2} \left[ \delta_{j+1/2} \left[ e_{1u}\,u \right] \right]  \right]  \biggr\}  &&\\ 
 \end{flalign*}
-Using the fact that the vertical diffusion coefficients are uniform, and that in 
-$z$-coordinate, the vertical scale factors do not depend on $i$ and $j$ so 
-that: $e_{3f} =e_{3u} =e_{3v} =e_{3t} $ and $e_{3w} =e_{3uw} =e_{3vw} $, 
-it follows:
+Using the fact that the vertical diffusion coefficients are uniform,
+and that in $z$-coordinate, the vertical scale factors do not depend on $i$ and $j$ so that:
+$e_{3f} =e_{3u} =e_{3v} =e_{3t} $ and $e_{3w} =e_{3uw} =e_{3vw} $, it follows:
 \begin{flalign*}
 \equiv A^{\,vm} \sum\limits_{i,j,k} \zeta \;\delta_k 
@@ -1398 +1390 @@
       \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} \right) \right)\; dv = 0    &&&\\
 \end{flalign*}
-and the square of the horizontal divergence decreases ($i.e.$ the horizontal 
-divergence is dissipated) if the vertical diffusion coefficient is uniform over the 
-whole domain:
+and the square of the horizontal divergence decreases ($i.e.$ the horizontal divergence is dissipated) if
+the vertical diffusion coefficient is uniform over the whole domain:
 
 \begin{flalign*}
@@ -1463 +1454 @@
 \label{sec:C.8}
 
-The numerical schemes used for tracer subgridscale physics are written such 
-that the heat and salt contents are conserved (equations in flux form). 
-Since a flux form is used to compute the temperature and salinity, 
-the quadratic form of these quantities ($i.e.$ their variance) globally tends to diminish. 
+The numerical schemes used for tracer subgridscale physics are written such that
+the heat and salt contents are conserved (equations in flux form).
+Since a flux form is used to compute the temperature and salinity,
+the quadratic form of these quantities ($i.e.$ their variance) globally tends to diminish.
 As for the advection term, there is conservation of mass only if the Equation Of Seawater is linear. 
 

NEMO/trunk/doc/latex/NEMO/subfiles/annex_D.tex

@@ -13 +13 @@
 
 
-A "model life" is more than ten years. Its software, composed of a few 
-hundred modules, is used by many people who are scientists or students 
-and do not necessarily know every aspect of computing very well. 
-Moreover, a well thought-out program is easier to read and understand, 
-less difficult to modify, produces fewer bugs and is easier to maintain. 
-Therefore, it is essential that the model development follows some rules :
+A "model life" is more than ten years.
+Its software, composed of a few hundred modules, is used by many people who are scientists or students and
+do not necessarily know every aspect of computing very well.
+Moreover, a well thought-out program is easier to read and understand, less difficult to modify,
+produces fewer bugs and is easier to maintain.
+Therefore, it is essential that the model development follows some rules:
 
 - well planned and designed
@@ -32 +32 @@
 - flexible.
 
-To satisfy part of these aims, \NEMO is written with a coding standard which 
-is close to the ECMWF rules, named DOCTOR \citep{Gibson_TR86}. 
-These rules present some advantages like :
+To satisfy part of these aims, \NEMO is written with a coding standard which is close to the ECMWF rules,
+named DOCTOR \citep{Gibson_TR86}. 
+These rules present some advantages like:
 
 - to provide a well presented program
 
-- to use rules for variable names which allow recognition of their type 
+- to use rules for variable names which allow recognition of their type
 (integer, real, parameter, local or shared variables, etc. ). 
 
@@ -49 +49 @@
 \label{sec:D_structure}
 
-Each program begins with a set of headline comments containing :
+Each program begins with a set of headline comments containing:
 
 - the program title
@@ -65 +65 @@
 - the author name(s), the date of creation and any updates.
 
-- Each program is split into several well separated sections and 
+- Each program is split into several well separated sections and
 sub-sections with an underlined title and specific labelled statements.
 
 - A program has not more than 200 to 300 lines.
 
-A template of a module style can be found on the NEMO depository 
-in the following file : NEMO/OPA\_SRC/module\_example.
+A template of a module style can be found on the NEMO depository in the following file:
+NEMO/OPA\_SRC/module\_example.
 % ================================================================
 % Coding conventions
@@ -78 +78 @@
 \label{sec:D_coding}
 
-- Use of the universal language \textsc{Fortran} 90, and try to avoid obsolescent
-features like statement functions, do not use GO TO and EQUIVALENCE statements.
-
-- A continuation line begins with the character {\&} indented by three spaces 
-compared to the previous line, while the previous line ended with the character {\&}.
-
-- All the variables must be declared. The code is usually compiled with implicit none.
+- Use of the universal language \textsc{Fortran} 90, and try to avoid obsolescent features like statement functions,
+do not use GO TO and EQUIVALENCE statements.
+
+- A continuation line begins with the character {\&} indented by three spaces compared to the previous line,
+while the previous line ended with the character {\&}.
+
+- All the variables must be declared.
+The code is usually compiled with implicit none.
  
-- Never use continuation lines in the declaration of a variable. When searching a 
-variable in the code through a \textit{grep} command, the declaration line will be found.
-
-- In the declaration of a PUBLIC variable, the comment part at the end of the line 
-should start with the two characters "\verb?!:?". the following UNIX command, \\
+- Never use continuation lines in the declaration of a variable.
+When searching a variable in the code through a \textit{grep} command, the declaration line will be found.
+
+- In the declaration of a PUBLIC variable, the comment part at the end of the line should start with
+the two characters "\verb?!:?".
+The following UNIX command, \\
 \verb?grep var_name *90 \ grep \!: ? \\
 will display the module name and the line where the var\_name declaration is.
 
-- Always use a three spaces indentation in DO loop, CASE, or IF-ELSEIF-ELSE-ENDIF 
-statements.
+- Always use a three spaces indentation in DO loop, CASE, or IF-ELSEIF-ELSE-ENDIF statements.
 
 - use a space after a comma, except when it appears to separate the indices of an array.
@@ -109 +110 @@
 \label{sec:D_naming}
 
-The purpose of the naming conventions is to use prefix letters to classify 
-model variables. These conventions allow the variable type to be easily 
-known and rapidly identified. The naming conventions are summarised 
-in the Table below:
+The purpose of the naming conventions is to use prefix letters to classify model variables.
+These conventions allow the variable type to be easily known and rapidly identified.
+The naming conventions are summarised in the Table below:
 
 
@@ -192 +192 @@
 %--------------------------------------------------------------------------------------------------------------
 
-N.B.   Parameter here, in not only parameter in the \textsc{Fortran} acceptation, it is also used for code variables 
-that are read in namelist and should never been modified during a simulation. 
+N.B. Parameter here, in not only parameter in the \textsc{Fortran} acceptation,
+it is also used for code variables that are read in namelist and should never been modified during a simulation. 
 It is the case, for example, for the size of a domain (jpi,jpj,jpk).
 

NEMO/trunk/doc/latex/NEMO/subfiles/annex_E.tex

@@ -10 +10 @@
 \newpage
 $\ $\newline    % force a new ligne
- 
- This appendix some on going consideration on algorithms used or planned to be used
-in \NEMO. 
+
+This appendix some on going consideration on algorithms used or planned to be used in \NEMO. 
 
 $\ $\newline    % force a new ligne
@@ -22 +21 @@
 \label{sec: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 Cell Averaged 
-QUICK scheme (Quadratic Upstream Interpolation for Convective 
-Kinematics). For example, in the $i$-direction :
+The UBS advection scheme is an upstream biased third order scheme based on
+an upstream-biased parabolic interpolation.
+It is also known as Cell Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective Kinematics).
+For example, in the $i$-direction:
 \begin{equation} \label{eq:tra_adv_ubs2}
 \tau _u^{ubs} = \left\{  \begin{aligned}
@@ -38 +37 @@
 - \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
 \end{equation}
-where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and 
-$\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. 
-By choosing this expression for $\tau "$ we consider a fourth order approximation 
-of $\partial_i^2$ with a constant i-grid spacing ($\Delta i=1$). 
+where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and
+$\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$.
+By choosing this expression for $\tau "$ we consider a fourth order approximation of $\partial_i^2$ with
+a constant i-grid spacing ($\Delta i=1$).
 
 Alternative choice: introduce the scale factors:  
@@ -47 +46 @@
 
 
-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}. 
-It is a relatively good compromise between accuracy and smoothness. It is 
-not a \emph{positive} scheme meaning false extrema are permitted but the 
-amplitude of such are significantly reduced over the centred second order 
-method. Nevertheless it is not recommended to apply it 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. 
-It has therefore been preferred to evaluate the vertical flux using the TVD 
-scheme when \np{ln\_traadv\_ubs}\forcode{ = .true.}.
-
-For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds 
-to a second order centred scheme is evaluated using the \textit{now} velocity 
-(centred in time) while the second term which is the diffusive part of the scheme, 
-is evaluated using the \textit{before} velocity (forward in time. This is discussed 
-by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. UBS and QUICK 
-schemes only differ by one coefficient. Substituting 1/6 with 1/8 in 
-(\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 
-This option is not available through a namelist parameter, since the 1/6 
-coefficient is hard coded. Nevertheless it is quite easy to make the 
-substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme
-
-NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can 
-be controlled by vertical advection (not vertical diffusion which is usually 
-solved using an implicit scheme). Computer time can be saved by using a 
-time-splitting technique on vertical advection. This possibility have been 
-implemented and validated in ORCA05-L301. It is not currently offered in the 
-current reference version. 
-
-NB 2 : In a forthcoming release four options will be proposed for the 
-vertical component used in the UBS scheme. $\tau _w^{ubs}$ will 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 
-parabolic splines following \citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, 
-or  \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an 
-eight-order accurate conventional scheme.
-
-NB 3 : It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows:
+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}.
+It is a relatively good compromise between accuracy and smoothness.
+It is not a \emph{positive} scheme meaning false extrema are permitted but
+the amplitude of such are significantly reduced over the centred second order method.
+Nevertheless it is not recommended to apply it 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.
+It has therefore been preferred to evaluate the vertical flux using the TVD scheme when
+\np{ln\_traadv\_ubs}\forcode{ = .true.}.
+
+For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds to
+a second order centred scheme is evaluated using the \textit{now} velocity (centred in time) while
+the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity
+(forward in time).
+This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme.
+UBS and QUICK schemes only differ by one coefficient.
+Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
+This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded.
+Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme.
+
+NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can be controlled by vertical advection
+(not vertical diffusion which is usually solved using an implicit scheme).
+Computer time can be saved by using a time-splitting technique on vertical advection.
+This possibility have been implemented and validated in ORCA05-L301.
+It is not currently offered in the current reference version. 
+
+NB 2: In a forthcoming release four options will be proposed for the vertical component used in the UBS scheme.
+$\tau _w^{ubs}$ will be evaluated using either \textit{(a)} a centered $2^{nd}$ order scheme,
+or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative parabolic splines following
+\citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS.
+The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme.
+
+NB 3: It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows:
 \begin{equation} \label{eq:tra_adv_ubs2}
 \tau _u^{ubs} = \left\{  \begin{aligned}
@@ -102 +96 @@
 \end{split}
 \end{equation}
-\autoref{eq:tra_adv_ubs2} has several advantages. First it clearly evidence that 
-the UBS scheme is based on the fourth order scheme to which is added an 
-upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order 
-part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order 
-part as stated above using \autoref{eq:tra_adv_ubs}. Third, the diffusive term is 
-in fact a biharmonic operator with a eddy coefficient with is simply proportional 
-to the velocity.
+\autoref{eq:tra_adv_ubs2} has several advantages.
+First it clearly evidences that the UBS scheme is based on the fourth order scheme to which
+is added an upstream biased diffusive term.
+Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step,
+not only the $2^{th}$ order part as stated above using \autoref{eq:tra_adv_ubs}.
+Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient which
+is simply proportional to the velocity.
 
 laplacian diffusion:
@@ -135 +129 @@
 with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$, 
 $i.e.$ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$
-it comes :
+it comes:
 \begin{equation} \label{eq:tra_ldf_lap}
 \begin{split}
@@ -163 +157 @@
 \end{split}
 \end{equation}
-if the velocity is uniform ($i.e.$ $|u|=cst$) and choosing $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right]$
+if the velocity is uniform ($i.e.$ $|u|=cst$) and
+choosing $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right]$
 
 sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$):
@@ -191 +186 @@
 \label{sec:LF}
 
-We adopt the following semi-discrete notation for time derivative. Given the values of a variable $q$ at successive time step, the time derivation and averaging operators at the mid time step are:
+We adopt the following semi-discrete notation for time derivative.
+Given the values of a variable $q$ at successive time step,
+the time derivation and averaging operators at the mid time step are:
 \begin{subequations} \label{eq:dt_mt}
 \begin{align}
@@ -198 +195 @@
 \end{align}
 \end{subequations}
-As for space operator, the adjoint of the derivation and averaging time operators are 
-$\delta_t^*=\delta_{t+\rdt/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$
-, respectively. 
+As for space operator,
+the adjoint of the derivation and averaging time operators are $\delta_t^*=\delta_{t+\rdt/2}$ and
+$\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$, respectively.
 
 The Leap-frog time stepping given by \autoref{eq:DOM_nxt} can be defined as:
@@ -208 +205 @@
       =         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt}
 \end{equation} 
-Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$, not $2\rdt$ 
-as it can be found sometime in literature. 
-The leap-Frog time stepping is a second order centered scheme. As such it respects 
-the quadratic invariant in integral forms, $i.e.$ the following continuous property,
+Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$,
+not $2\rdt$ as it can be found sometimes in literature.
+The leap-Frog time stepping is a second order centered scheme.
+As such it respects the quadratic invariant in integral forms, $i.e.$ the following continuous property,
 \begin{equation} \label{eq:Energy}
 \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} 
@@ -217 +214 @@
    =  \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) ,
 \end{equation}
-is satisfied in discrete form. Indeed, 
+is satisfied in discrete form.
+Indeed, 
 \begin{equation} \begin{split}
 \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} 
@@ -228 +226 @@
       \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right)
 \end{split} \end{equation}
-NB here pb of boundary condition when applying the adjoin! In space, setting to 0 
-the quantity in land area is sufficient to get rid of the boundary condition 
-(equivalently of the boundary value of the integration by part). In time this boundary 
-condition is not physical and \textbf{add something here!!!}
+NB here pb of boundary condition when applying the adjoint!
+In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition 
+(equivalently of the boundary value of the integration by part).
+In time this boundary condition is not physical and \textbf{add something here!!!}
 
 
@@ -249 +247 @@
 \subsection{Griffies iso-neutral diffusion operator}
 
-Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98}
-scheme, but is formulated within the \NEMO framework ($i.e.$ using scale 
-factors rather than grid-size and having a position of $T$-points that is not 
-necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}).
-
-In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO, 
-the off-diagonal terms of the small angle diffusion tensor contain several double 
-spatial averages of a gradient, for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$. 
-It is apparent that the combination of a $k$ average and a $k$ derivative of the tracer 
-allows for the presence of grid point oscillation structures that will be invisible 
-to the operator. These structures are \textit{computational modes}. They
-will not be damped by the iso-neutral operator, and even possibly amplified by it. 
-In other word, the operator applied to a tracer does not warranties the decrease of 
-its global average variance. To circumvent this, we have introduced a smoothing of 
-the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). Nevertheless, this technique 
-works fine for $T$ and $S$ as they are active tracers ($i.e.$ they enter the computation 
-of density), but it does not work for a passive tracer.   \citep{Griffies_al_JPO98} introduce 
-a different way to discretise the off-diagonal terms that nicely solve the problem. 
-The idea is to get rid of combinations of an averaged in one direction combined 
-with a derivative in the same direction by considering triads. For example in the 
-(\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows:
+Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98} scheme,
+but is formulated within the \NEMO framework
+($i.e.$ using scale factors rather than grid-size and having a position of $T$-points that
+is not necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}).
+
+In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO,
+the off-diagonal terms of the small angle diffusion tensor contain several double spatial averages of a gradient,
+for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$.
+It is apparent that the combination of a $k$ average and a $k$ derivative of the tracer allows for
+the presence of grid point oscillation structures that will be invisible to the operator.
+These structures are \textit{computational modes}.
+They will not be damped by the iso-neutral operator, and even possibly amplified by it.
+In other word, the operator applied to a tracer does not warranties the decrease of its global average variance.
+To circumvent this, we have introduced a smoothing of the slopes of the iso-neutral surfaces
+(see \autoref{chap:LDF}).
+Nevertheless, this technique works fine for $T$ and $S$ as they are active tracers
+($i.e.$ they enter the computation of density), but it does not work for a passive tracer.
+\citep{Griffies_al_JPO98} introduce a different way to discretise the off-diagonal terms that
+nicely solve the problem.
+The idea is to get rid of combinations of an averaged in one direction combined with
+a derivative in the same direction by considering triads.
+For example in the (\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows:
 \begin{equation} \label{eq:Gf_triads}
 _i^k \mathbb{T}_{i_p}^{k_p} (T)
@@ -277 +276 @@
                              \right)
 \end{equation}
-where the indices $i_p$ and $k_p$ define the four triads and take the following value: 
-$i_p = -1/2$ or $1/2$ and $k_p = -1/2$ or $1/2$, 
-$b_u= e_{1u}\,e_{2u}\,e_{3u}$ is the volume of $u$-cells, 
+where the indices $i_p$ and $k_p$ define the four triads and take the following value:
+$i_p = -1/2$ or $1/2$ and $k_p = -1/2$ or $1/2$,
+$b_u= e_{1u}\,e_{2u}\,e_{3u}$ is the volume of $u$-cells,
 $A_i^k$ is the lateral eddy diffusivity coefficient defined at $T$-point,
-and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad :
+and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad:
 \begin{equation} \label{eq:Gf_slopes}
 _i^k \mathbb{R}_{i_p}^{k_p} 
@@ -288 +287 @@
 {\left(\alpha / \beta \right)_i^k  \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }
 \end{equation}
-Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of 
-multiplying the temperature derivative by $\alpha$ and the salinity derivative 
-by $\beta$. This is more efficient as the ratio $\alpha / \beta$ can to be 
-evaluated directly.
-
-Note that in \autoref{eq:Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of 
-${b_{uw}}_{\,i+i_p}^{\,k+k_p}$. This choice has been motivated by the decrease 
-of tracer variance and the presence of partial cell at the ocean bottom 
-(see \autoref{apdx:Gf_operator}).
+Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of
+multiplying the temperature derivative by $\alpha$ and the salinity derivative by $\beta$.
+This is more efficient as the ratio $\alpha / \beta$ can to be evaluated directly.
+
+Note that in \autoref{eq:Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of ${b_{uw}}_{\,i+i_p}^{\,k+k_p}$.
+This choice has been motivated by the decrease of tracer variance and
+the presence of partial cell at the ocean bottom (see \autoref{apdx:Gf_operator}).
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht] \begin{center}
 \includegraphics[width=0.70\textwidth]{Fig_ISO_triad}
-\caption{  \protect\label{fig:ISO_triad}   
-Triads used in the Griffies's like iso-neutral diffision scheme for 
-$u$-component (upper panel) and $w$-component (lower panel).}
+\caption{  \protect\label{fig:ISO_triad}
+  Triads used in the Griffies's like iso-neutral diffision scheme for
+  $u$-component (upper panel) and $w$-component (lower panel).}
 \end{center}
 \end{figure}
@@ -309 +306 @@
 
 The four iso-neutral fluxes associated with the triads are defined at $T$-point. 
-They take the following expression :
+They take the following expression:
 \begin{flalign} \label{eq:Gf_fluxes}
 \begin{split}
@@ -320 +317 @@
 \end{flalign}
 
-The resulting iso-neutral fluxes at $u$- and $w$-points are then given by the 
-sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}):
+The resulting iso-neutral fluxes at $u$- and $w$-points are then given by
+the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}):
 \begin{flalign} \label{eq:iso_flux} 
 \textbf{F}_{iso}(T) 
@@ -350 +347 @@
 %    \end{pmatrix}      
 \end{flalign}
-resulting in a iso-neutral diffusion tendency on temperature given by the divergence 
-of the sum of all the four triad fluxes :
+resulting in a iso-neutral diffusion tendency on temperature given by
+the divergence of the sum of all the four triad fluxes:
 \begin{equation} \label{eq:Gf_operator}
 D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{  
@@ -359 +356 @@
 where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. 
 
-This expression of the iso-neutral diffusion has been chosen in order to satisfy 
-the following six properties:
+This expression of the iso-neutral diffusion has been chosen in order to satisfy the following six properties:
 \begin{description}
-\item[$\bullet$ horizontal diffusion] The discretization of the diffusion operator 
-recovers the traditional five-point Laplacian in the limit of flat iso-neutral direction :
+\item[$\bullet$ horizontal diffusion]
+  The discretization of the diffusion operator recovers the traditional five-point Laplacian in
+  the limit of flat iso-neutral direction:
 \begin{equation} \label{eq:Gf_property1a}
 D_l^T = \frac{1}{b_T}  \ \delta_{i} 
@@ -371 +368 @@
 \end{equation}
 
-\item[$\bullet$ implicit treatment in the vertical]  In the diagonal term associated 
-with the vertical divergence of the iso-neutral fluxes (i.e. the term associated 
-with a second order vertical derivative) appears only tracer values associated 
-with a single water column. This is of paramount importance since it means
-that the implicit in time algorithm for solving the vertical diffusion equation can 
-be used to evaluate this term. It is a necessity since the vertical eddy diffusivity 
-associated with this term,  
+\item[$\bullet$ implicit treatment in the vertical]
+  In the diagonal term associated with the vertical divergence of the iso-neutral fluxes
+  (i.e. the term associated with a second order vertical derivative)
+  appears only tracer values associated with a single water column.
+  This is of paramount importance since it means that
+  the implicit in time algorithm for solving the vertical diffusion equation can be used to evaluate this term.
+  It is a necessity since the vertical eddy diffusivity associated with this term,  
 \begin{equation}
     \sum_{\substack{i_p, \,k_p}} \left\{  
@@ -385 +382 @@
 can be quite large.
 
-\item[$\bullet$ pure iso-neutral operator]  The iso-neutral flux of locally referenced 
-potential density is zero, $i.e.$
+\item[$\bullet$ pure iso-neutral operator]
+  The iso-neutral flux of locally referenced potential density is zero, $i.e.$
 \begin{align} \label{eq:Gf_property2}
 \begin{matrix}
@@ -397 +394 @@
 \end{matrix}
 \end{align}
-This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$ 
-and the definition of the triads' slopes \autoref{eq:Gf_slopes}.
-
-\item[$\bullet$ conservation of tracer] The iso-neutral diffusion term conserve the 
-total tracer content, $i.e.$
+This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$ and
+the definition of the triads' slopes \autoref{eq:Gf_slopes}.
+
+\item[$\bullet$ conservation of tracer]
+  The iso-neutral diffusion term conserve the total tracer content, $i.e.$
 \begin{equation} \label{eq:Gf_property1}
 \sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0
 \end{equation}
-This property is trivially satisfied since the iso-neutral diffusive operator 
-is written in flux form.
-
-\item[$\bullet$ decrease of tracer variance] The iso-neutral diffusion term does 
-not increase the total tracer variance, $i.e.$
+This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form.
+
+\item[$\bullet$ decrease of tracer variance]
+  The iso-neutral diffusion term does not increase the total tracer variance, $i.e.$
 \begin{equation} \label{eq:Gf_property1}
 \sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0
 \end{equation}
-The property is demonstrated in the \autoref{apdx:Gf_operator}. It is a 
-key property for a diffusion term. It means that the operator is also a dissipation 
-term, $i.e.$ it is a sink term for the square of the quantity on which it is applied. 
-It therfore ensure that, when the diffusivity coefficient is large enough, the field 
-on which it is applied become free of grid-point noise.
-
-\item[$\bullet$ self-adjoint operator] The iso-neutral diffusion operator is self-adjoint, 
-$i.e.$
+The property is demonstrated in the \autoref{apdx:Gf_operator}.
+It is a key property for a diffusion term.
+It means that the operator is also a dissipation term,
+$i.e.$ it is a sink term for the square of the quantity on which it is applied.
+It therfore ensures that, when the diffusivity coefficient is large enough,
+the field on which it is applied become free of grid-point noise.
+
+\item[$\bullet$ self-adjoint operator]
+  The iso-neutral diffusion operator is self-adjoint, $i.e.$
 \begin{equation} \label{eq:Gf_property1}
 \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 
 \end{equation}
-In other word, there is no needs to develop a specific routine from the adjoint of this 
-operator. We just have to apply the same routine. This properties can be demonstrated 
-quite easily in a similar way the "non increase of tracer variance" property has been 
-proved (see \autoref{apdx:Gf_operator}).
+In other word, there is no needs to develop a specific routine from the adjoint of this operator.
+We just have to apply the same routine.
+This properties can be demonstrated quite easily in a similar way the "non increase of tracer variance" property
+has been proved (see \autoref{apdx:Gf_operator}).
 \end{description}
 
@@ -437 +434 @@
 \subsection{Eddy induced velocity and skew flux formulation}
 
-When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined), 
-an additional advection term is added. The associated velocity is the so called 
-eddy induced velocity, 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.$ \autoref{eq:ldfslp_geo} 
-is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo}
-+ \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. 
+When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined),
+an additional advection term is added.
+The associated velocity is the so called eddy induced velocity,
+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.$ \autoref{eq:ldfslp_geo} is used in $z$-coordinate,
+and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. 
 
 The eddy induced velocity is given by: 
@@ -456 +453 @@
 \end{split}
 \end{equation}
-where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the 
-slopes between the iso-neutral and the geopotential surfaces.
+where $A_{e}$ is the eddy induced velocity coefficient,
+and $r_i$ and $r_j$ the slopes between the iso-neutral and the geopotential surfaces.
 %%gm wrong: to be modified with 2 2D streamfunctions
- In other words,
-the eddy induced velocity can be derived from a vector streamfuntion, $\phi$, which
-is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^*  = \textbf{k} \times \nabla \phi$
+In other words, the eddy induced velocity can be derived from a vector streamfuntion, $\phi$,
+which is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^*  = \textbf{k} \times \nabla \phi$.
 %%end gm
 
-A traditional way to implement this additional advection is to add it to the eulerian 
-velocity prior to compute the tracer advection. This allows us to take advantage of 
-all the advection schemes offered for the tracers (see \autoref{sec:TRA_adv}) and not just 
-a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers 
-where \emph{positivity} of the advection scheme is of paramount importance. 
+A traditional way to implement this additional advection is to add it to the eulerian velocity prior to
+compute the tracer advection.
+This allows us to take advantage of all the advection schemes offered for the tracers
+(see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme.
+This is particularly useful for passive tracers where
+\emph{positivity} of the advection scheme is of paramount importance. 
 % give here the expression using the triads. It is different from the one given in \autoref{eq:ldfeiv}
 % see just below a copy of this equation:
@@ -490 +487 @@
 \end{equation}
 
-\citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, 
-the so-called skew form. It is based on a transformation of the advective fluxes 
-using the non-divergent nature of the eddy induced velocity. 
-For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be 
-transformed as follows:
+\citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form.
+It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity.
+For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows:
 \begin{flalign*}
 \begin{split}
@@ -519 +514 @@
 \end{split}
 \end{flalign*}
-and since the eddy induces velocity field is no-divergent, we end up with the skew 
-form of the eddy induced advective fluxes:
+and since the eddy induces velocity field is no-divergent,
+we end up with the skew form of the eddy induced advective fluxes:
 \begin{equation} \label{eq:eiv_skew_continuous}
 \textbf{F}_{eiv}^T = \begin{pmatrix} 
@@ -527 +522 @@
                                  \end{pmatrix}
 \end{equation}
-The tendency associated with eddy induced velocity is then simply the divergence 
-of the \autoref{eq:eiv_skew_continuous} fluxes. It naturally conserves the tracer 
-content, as it is expressed in flux form and, as the advective form, it preserve the 
-tracer variance. Another interesting property of \autoref{eq:eiv_skew_continuous} 
-form is that when $A=A_e$, a simplification occurs in the sum of the iso-neutral 
-diffusion and eddy induced velocity terms:
+The tendency associated with eddy induced velocity is then simply the divergence of
+the \autoref{eq:eiv_skew_continuous} fluxes.
+It naturally conserves the tracer content, as it is expressed in flux form and,
+as the advective form, it preserves the tracer variance.
+Another interesting property of \autoref{eq:eiv_skew_continuous} form is that when $A=A_e$,
+a simplification occurs in the sum of the iso-neutral diffusion and eddy induced velocity terms:
 \begin{flalign} \label{eq:eiv_skew+eiv_continuous}
 \textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &= 
@@ -549 +544 @@
 \end{pmatrix}
 \end{flalign}
-The horizontal component reduces to the one use for an horizontal laplacian 
-operator and the vertical one keep the same complexity, but not more. This property
-has been used to reduce the computational time \citep{Griffies_JPO98}, but it is 
-not of practical use as usually $A \neq A_e$. Nevertheless this property can be used to 
-choose a discret form of  \autoref{eq:eiv_skew_continuous} which is consistent with the 
-iso-neutral operator \autoref{eq:Gf_operator}. Using the slopes \autoref{eq:Gf_slopes} 
-and defining $A_e$ at $T$-point($i.e.$ as $A$, the eddy diffusivity coefficient),
-the resulting discret form is given by:
+The horizontal component reduces to the one use for an horizontal laplacian operator and
+the vertical one keeps the same complexity, but not more.
+This property has been used to reduce the computational time \citep{Griffies_JPO98},
+but it is not of practical use as usually $A \neq A_e$.
+Nevertheless this property can be used to choose a discret form of \autoref{eq:eiv_skew_continuous} which
+is consistent with the iso-neutral operator \autoref{eq:Gf_operator}.
+Using the slopes \autoref{eq:Gf_slopes} and defining $A_e$ at $T$-point($i.e.$ as $A$,
+the eddy diffusivity coefficient), the resulting discret form is given by:
 \begin{equation} \label{eq:eiv_skew}  
 \textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left( \begin{aligned}                                
@@ -569 +564 @@
 \end{equation}
 Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells. 
-In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces 
-must be added to $\mathbb{R}$ for the discret form to be exact. 
-
-Such a choice of discretisation is consistent with the iso-neutral operator as it uses the 
-same definition for the slopes. It also ensures the conservation of the tracer variance 
-(see Appendix \autoref{apdx:eiv_skew}), $i.e.$ it does not include a diffusive component 
-but is a "pure" advection term.
+In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces must be added to
+$\mathbb{R}$ for the discret form to be exact. 
+
+Such a choice of discretisation is consistent with the iso-neutral operator as
+it uses the same definition for the slopes.
+It also ensures the conservation of the tracer variance (see Appendix \autoref{apdx:eiv_skew}),
+$i.e.$ it does not include a diffusive component but is a "pure" advection term.
 
 
@@ -591 +586 @@
 This part will be moved in an Appendix.
 
-The continuous property to be demonstrated is :
+The continuous property to be demonstrated is:
 \begin{align*}
 \int_D  D_l^T \; T \;dv   \leq 0
@@ -642 +637 @@
 %
 \allowdisplaybreaks
-\intertext{The summation is done over all $i$ and $k$ indices, it is therefore possible to introduce a shift of $-1$ either in $i$ or $k$ direction in order to regroup all the terms of the summation by triad at a ($i$,$k$) point. In other words, we regroup all the terms in the neighbourhood  that contain a triad at the same ($i$,$k$) indices. It becomes: }
+  \intertext{The summation is done over all $i$ and $k$ indices,
+  it is therefore possible to introduce a shift of $-1$ either in $i$ or $k$ direction in order to
+  regroup all the terms of the summation by triad at a ($i$,$k$) point.
+  In other words, we regroup all the terms in the neighbourhood that contain a triad at the same ($i$,$k$) indices.
+  It becomes: }
 %
 &\equiv -\sum_{i,k}
@@ -672 +671 @@
 %
 \allowdisplaybreaks
-\intertext{Then outing in factor the triad in each of the four terms of the summation and substituting the triads by their expression given in \autoref{eq:Gf_triads}. It becomes: }
+  \intertext{Then outing in factor the triad in each of the four terms of the summation and
+  substituting the triads by their expression given in \autoref{eq:Gf_triads}.
+  It becomes: }
 %
 &\equiv -\sum_{i,k}
@@ -710 +711 @@
 The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities.
 
-Note that, if instead of multiplying $D_l^T$ by $T$, we were using another tracer field, let say $S$, then the previous demonstration would have let to:
+Note that, if instead of multiplying $D_l^T$ by $T$, we were using another tracer field, let say $S$,
+then the previous demonstration would have let to:
 \begin{align*}
 \int_D  S \; D_l^T  \;dv &\equiv  \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\}    \\
@@ -729 +731 @@
 &\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\}   
 \end{align*} 
-This means that the iso-neutral operator is self-adjoint. There is no need to develop a specific to obtain it.
+This means that the iso-neutral operator is self-adjoint.
+There is no need to develop a specific to obtain it.
 
 
@@ -745 +748 @@
 This have to be moved in an Appendix.
 
-The continuous property to be demonstrated is :
+The continuous property to be demonstrated is:
 \begin{align*}
 \int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0
@@ -797 +800 @@
 \end{matrix}   
 \end{align*}
-The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the 
-same but of opposite signs, they cancel out. 
-Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$. 
-The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the 
-same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$ 
-they cancel out with the neighbouring grid points. 
-Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the 
-$i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the 
-tracer is preserved by the discretisation of the skew fluxes.
+The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs,
+they cancel out. 
+Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$.
+The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the same but both of opposite signs and
+shifted by 1 in $k$ direction.
+When summing over $k$ they cancel out with the neighbouring grid points.
+Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the $i$ direction.
+Therefore the sum over the domain is zero,
+$i.e.$ the variance of the tracer is preserved by the discretisation of the skew fluxes.
 
 \end{document}

NEMO/trunk/doc/latex/NEMO/subfiles/annex_iso.tex

@@ -15 +15 @@
 %---------------------------------------------------------------------------------------------------------
 
-Two scheme are available to perform the iso-neutral diffusion. 
-If the namelist logical \np{ln\_traldf\_triad} is set true, 
-\NEMO updates both active and passive tracers using the Griffies triad representation 
-of iso-neutral diffusion and the eddy-induced advective skew (GM) fluxes. 
-If the namelist logical \np{ln\_traldf\_iso} is set true, 
-the filtered version of Cox's original scheme (the Standard scheme) is employed (\autoref{sec:LDF_slp}). 
-In the present implementation of the Griffies scheme, 
+Two scheme are available to perform the iso-neutral diffusion.
+If the namelist logical \np{ln\_traldf\_triad} is set true,
+\NEMO updates both active and passive tracers using the Griffies triad representation of iso-neutral diffusion and
+the eddy-induced advective skew (GM) fluxes.
+If the namelist logical \np{ln\_traldf\_iso} is set true,
+the filtered version of Cox's original scheme (the Standard scheme) is employed (\autoref{sec:LDF_slp}).
+In the present implementation of the Griffies scheme,
 the advective skew fluxes are implemented even if \np{ln\_traldf\_eiv} is false.
 
-Values of iso-neutral diffusivity and GM coefficient are set as
-described in \autoref{sec:LDF_coef}. Note that when GM fluxes are used, 
-the eddy-advective (GM) velocities are output for diagnostic purposes using xIOS, 
+Values of iso-neutral diffusivity and GM coefficient are set as described in \autoref{sec:LDF_coef}.
+Note that when GM fluxes are used, the eddy-advective (GM) velocities are output for diagnostic purposes using xIOS,
 even though the eddy advection is accomplished by means of the skew fluxes.
 
@@ -32 +31 @@
 The options specific to the Griffies scheme include:
 \begin{description}[font=\normalfont]
-\item[\np{ln\_triad\_iso}] See \autoref{sec:taper}. If this is set false (the default), then
-  `iso-neutral' mixing is accomplished within the surface mixed-layer
-  along slopes linearly decreasing with depth from the value immediately below
-  the mixed-layer to zero (flat) at the surface (\autoref{sec:lintaper}). 
-  This is the same treatment as used in the default implementation \autoref{subsec:LDF_slp_iso}; \autoref{fig:eiv_slp}.  
-  Where \np{ln\_triad\_iso} is set true, the vertical skew flux is further reduced 
-  to ensure no vertical buoyancy flux, giving an almost pure
-  horizontal diffusive tracer flux within the mixed layer. This is similar to
-  the tapering suggested by \citet{Gerdes1991}. See \autoref{subsec:Gerdes-taper}
-\item[\np{ln\_botmix\_triad}] See \autoref{sec:iso_bdry}. 
+\item[\np{ln\_triad\_iso}]
+  See \autoref{sec:taper}.
+  If this is set false (the default),
+  then `iso-neutral' mixing is accomplished within the surface mixed-layer along slopes linearly decreasing with
+  depth from the value immediately below the mixed-layer to zero (flat) at the surface (\autoref{sec:lintaper}).
+  This is the same treatment as used in the default implementation
+  \autoref{subsec:LDF_slp_iso}; \autoref{fig:eiv_slp}.
+  Where \np{ln\_triad\_iso} is set true,
+  the vertical skew flux is further reduced to ensure no vertical buoyancy flux,
+  giving an almost pure horizontal diffusive tracer flux within the mixed layer.
+  This is similar to the tapering suggested by \citet{Gerdes1991}. See \autoref{subsec:Gerdes-taper}
+\item[\np{ln\_botmix\_triad}]
+  See \autoref{sec:iso_bdry}. 
   If this is set false (the default) then the lateral diffusive fluxes
   associated with triads partly masked by topography are neglected. 
   If it is set true, however, then these lateral diffusive fluxes are applied, 
   giving smoother bottom tracer fields at the cost of introducing diapycnal mixing.
-\item[\np{rn\_sw\_triad}]  blah blah to be added....
+\item[\np{rn\_sw\_triad}]
+  blah blah to be added....
 \end{description}
 The options shared with the Standard scheme include:
@@ -56 +59 @@
 \section{Triad formulation of iso-neutral diffusion}
 \label{sec:iso}
-We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, 
+We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98},
 but formulated within the \NEMO framework, using scale factors rather than grid-sizes.
 
 \subsection{Iso-neutral diffusion operator}
-The iso-neutral second order tracer diffusive operator for small
-angles between iso-neutral surfaces and geopotentials is given by
-\autoref{eq:iso_tensor_1}:
+The iso-neutral second order tracer diffusive operator for small angles between
+iso-neutral surfaces and geopotentials is given by \autoref{eq:iso_tensor_1}:
 \begin{subequations} \label{eq:iso_tensor_1}
   \begin{equation}
@@ -94 +96 @@
 %  {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\
 % \end{array} }} \right)
- Here \autoref{eq:PE_iso_slopes} 
+Here \autoref{eq:PE_iso_slopes} 
 \begin{align*}
   r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i}
@@ -104 +106 @@
     }{\partial k} \right)^{-1}
 \end{align*}
-is the $i$-component of the slope of the iso-neutral surface relative to the computational
-surface, and $r_2$ is the $j$-component.
-
-We will find it useful to consider the fluxes per unit area in $i,j,k$
-space; we write
+is the $i$-component of the slope of the iso-neutral surface relative to the computational surface,
+and $r_2$ is the $j$-component.
+
+We will find it useful to consider the fluxes per unit area in $i,j,k$ space; we write
 \begin{equation}
   \label{eq:Fijk}
   \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right).
 \end{equation}
-Additionally, we will sometimes write the contributions towards the
-fluxes $\vect{f}$ and $\vect{F}_{\mathrm{iso}}$ from the component
-$R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$, with
-$f_{ij}=R_{ij}e_i^{-1}\partial T/\partial x_i$ (no summation) etc.
+Additionally, we will sometimes write the contributions towards the fluxes $\vect{f}$ and
+$\vect{F}_{\mathrm{iso}}$ from the component $R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$,
+with $f_{ij}=R_{ij}e_i^{-1}\partial T/\partial x_i$ (no summation) etc.
 
 The off-diagonal terms of the small angle diffusion tensor
-\autoref{eq:iso_tensor_1}, \autoref{eq:iso_tensor_2} produce skew-fluxes along the
-$i$- and $j$-directions resulting from the vertical tracer gradient:
+\autoref{eq:iso_tensor_1}, \autoref{eq:iso_tensor_2} produce skew-fluxes along
+the $i$- and $j$-directions resulting from the vertical tracer gradient:
 \begin{align}
   \label{eq:i13c}
@@ -129 +129 @@
 \end{align}
 
-The vertical diffusive flux associated with the $_{33}$
-component of the small angle diffusion tensor is
+The vertical diffusive flux associated with the $_{33}$ component of the small angle diffusion tensor is
 \begin{equation}
   \label{eq:i33c}
@@ -136 +135 @@
 \end{equation}
 
-Since there are no cross terms involving $r_1$ and $r_2$ in the above, we can
-consider the iso-neutral diffusive fluxes separately in the $i$-$k$ and $j$-$k$
-planes, just adding together the vertical components from each
-plane. The following description will describe the fluxes on the $i$-$k$
-plane.
-
-There is no natural discretization for the $i$-component of the
-skew-flux, \autoref{eq:i13c}, as
-although it must be evaluated at $u$-points, it involves vertical
-gradients (both for the tracer and the slope $r_1$), defined at
-$w$-points. Similarly, the vertical skew flux, \autoref{eq:i31c}, is evaluated at
-$w$-points but involves horizontal gradients defined at $u$-points.
+Since there are no cross terms involving $r_1$ and $r_2$ in the above,
+we can consider the iso-neutral diffusive fluxes separately in the $i$-$k$ and $j$-$k$ planes,
+just adding together the vertical components from each plane.
+The following description will describe the fluxes on the $i$-$k$ plane.
+
+There is no natural discretization for the $i$-component of the skew-flux, \autoref{eq:i13c},
+as although it must be evaluated at $u$-points,
+it involves vertical gradients (both for the tracer and the slope $r_1$), defined at $w$-points.
+Similarly, the vertical skew flux, \autoref{eq:i31c},
+is evaluated at $w$-points but involves horizontal gradients defined at $u$-points.
 
 \subsection{Standard discretization}
 The straightforward approach to discretize the lateral skew flux
-\autoref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995
-into OPA, \autoref{eq:tra_ldf_iso}, is to calculate a mean vertical
-gradient at the $u$-point from the average of the four surrounding
-vertical tracer gradients, and multiply this by a mean slope at the
-$u$-point, calculated from the averaged surrounding vertical density
-gradients. The total area-integrated skew-flux (flux per unit area in
-$ijk$ space) from tracer cell $i,k$
-to $i+1,k$, noting that the $e_{{3}_{i+1/2}^k}$ in the area
-$e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with
-the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer
-gradient, is then \autoref{eq:tra_ldf_iso}
+\autoref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA,
+\autoref{eq:tra_ldf_iso}, is to calculate a mean vertical gradient at the $u$-point from
+the average of the four surrounding vertical tracer gradients, and multiply this by a mean slope at the $u$-point,
+calculated from the averaged surrounding vertical density gradients.
+The total area-integrated skew-flux (flux per unit area in $ijk$ space) from tracer cell $i,k$ to $i+1,k$,
+noting that the $e_{{3}_{i+1/2}^k}$ in the area $e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with
+the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer gradient, is then \autoref{eq:tra_ldf_iso}
 \begin{equation*}
   \left(F_u^{13} \right)_{i+\hhalf}^k = \Alts_{i+\hhalf}^k
@@ -173 +166 @@
   \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}},
 \end{equation*}
-and here and in the following we drop the $^{lT}$ superscript from
-$\Alt$ for simplicity.
-Unfortunately the resulting combination $\overline{\overline{\delta_k
-    \bullet}}^{\,i,k}$ of a $k$ average and a $k$ difference %of the tracer
-reduces to $\bullet_{k+1}-\bullet_{k-1}$, so two-grid-point oscillations are
-invisible to this discretization of the iso-neutral operator. These
-\emph{computational modes} will not be damped by this operator, and
-may even possibly be amplified by it.  Consequently, applying this
-operator to a tracer does not guarantee the decrease of its
-global-average variance. To correct this, we introduced a smoothing of
-the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). This
-technique works for $T$ and $S$ in so far as they are active tracers
-($i.e.$ they enter the computation of density), but it does not work
-for a passive tracer.
+and here and in the following we drop the $^{lT}$ superscript from $\Alt$ for simplicity.
+Unfortunately the resulting combination $\overline{\overline{\delta_k\bullet}}^{\,i,k}$ of a $k$ average and
+a $k$ difference of the tracer reduces to $\bullet_{k+1}-\bullet_{k-1}$,
+so two-grid-point oscillations are invisible to this discretization of the iso-neutral operator.
+These \emph{computational modes} will not be damped by this operator, and may even possibly be amplified by it.
+Consequently, applying this operator to a tracer does not guarantee the decrease of its global-average variance.
+To correct this, we introduced a smoothing of the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}).
+This technique works for $T$ and $S$ in so far as they are active tracers
+($i.e.$ they enter the computation of density), but it does not work for a passive tracer.
 
 \subsection{Expression of the skew-flux in terms of triad slopes}
-\citep{Griffies_al_JPO98} introduce a different discretization of the
-off-diagonal terms that nicely solves the problem.
+\citep{Griffies_al_JPO98} introduce a different discretization of the off-diagonal terms that
+nicely solves the problem.
 % Instead of multiplying the mean slope calculated at the $u$-point by
 % the mean vertical gradient at the $u$-point,
@@ -203 +191 @@
  \end{center} \end{figure}
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
-They get the skew flux from the products of the vertical gradients at
-each $w$-point surrounding the $u$-point with the corresponding `triad'
-slope calculated from the lateral density gradient across the $u$-point
-divided by the vertical density gradient at the same $w$-point as the
-tracer gradient. See \autoref{fig:ISO_triad}a, where the thick lines
-denote the tracer gradients, and the thin lines the corresponding
-triads, with slopes $s_1, \dotsc s_4$. The total area-integrated
-skew-flux from tracer cell $i,k$ to $i+1,k$
+They get the skew flux from the products of the vertical gradients at each $w$-point surrounding the $u$-point with
+the corresponding `triad' slope calculated from the lateral density gradient across the $u$-point divided by
+the vertical density gradient at the same $w$-point as the tracer gradient.
+See \autoref{fig:ISO_triad}a, where the thick lines denote the tracer gradients,
+and the thin lines the corresponding triads, with slopes $s_1, \dotsc s_4$.
+The total area-integrated skew-flux from tracer cell $i,k$ to $i+1,k$
 \begin{multline}
   \label{eq:i13}
@@ -222 +208 @@
   _{k-\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}},
 \end{multline}
-where the contributions of the triad fluxes are weighted by areas
-$a_1, \dotsc a_4$, and $\Alts$ is now defined at the tracer points
-rather than the $u$-points. This discretization gives a much closer
-stencil, and disallows the two-point computational modes.
-
- The vertical skew flux \autoref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at the
-$w$-point $i,k+\hhalf$ is constructed similarly (\autoref{fig:ISO_triad}b)
-by multiplying lateral tracer gradients from each of the four
-surrounding $u$-points by the appropriate triad slope:
+where the contributions of the triad fluxes are weighted by areas $a_1, \dotsc a_4$,
+and $\Alts$ is now defined at the tracer points rather than the $u$-points.
+This discretization gives a much closer stencil, and disallows the two-point computational modes.
+
+The vertical skew flux \autoref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at
+the $w$-point $i,k+\hhalf$ is constructed similarly (\autoref{fig:ISO_triad}b) by
+multiplying lateral tracer gradients from each of the four surrounding $u$-points by the appropriate triad slope:
 \begin{multline}
   \label{eq:i31}
@@ -241 +225 @@
 
 We notate the triad slopes $s_i$ and $s'_i$ in terms of the `anchor point' $i,k$
-(appearing in both the vertical and lateral gradient), and the $u$- and
-$w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the
-triad as follows (see also \autoref{fig:ISO_triad}):
+(appearing in both the vertical and lateral gradient),
+and the $u$- and $w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the triad as follows
+(see also \autoref{fig:ISO_triad}):
 \begin{equation}
   \label{eq:R}
@@ -253 +237 @@
   { \alpha_i^k  \ \delta_{k+k_p}[T^i] - \beta_i^k \ \delta_{k+k_p}[S^i] }.
 \end{equation}
-In calculating the slopes of the local neutral surfaces, 
-the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad, 
+In calculating the slopes of the local neutral surfaces,
+the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad,
 while the metrics are calculated at the $u$- and $w$-points on the arms.
 
@@ -261 +245 @@
     \includegraphics[width=0.80\textwidth]{Fig_GRIFF_qcells}
     \caption{   \protect\label{fig:qcells}
-    Triad notation for quarter cells. $T$-cells are inside
-      boxes, while the  $i+\half,k$ $u$-cell is shaded in green and the
-      $i,k+\half$ $w$-cell is shaded in pink.}
+      Triad notation for quarter cells. $T$-cells are inside boxes,
+      while the  $i+\half,k$ $u$-cell is shaded in green and
+      the $i,k+\half$ $w$-cell is shaded in pink.}
   \end{center} \end{figure}
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:qcells}) with the quarter
-cell that is the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ $u$-cell and the $i,k+k_p$ $w$-cell. 
-Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:i13} and \autoref{eq:i31} in this notation, 
-we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. 
-Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) 
-to calculate the lateral flux along its $u$-arm, at $(i+i_p,k)$, 
-and then again as an $s'$ to calculate the vertical flux along its $w$-arm at $(i,k+k_p)$. 
-Each vertical area $a_i$ used to calculate the lateral flux and horizontal area $a'_i$ used 
-to calculate the vertical flux can also be identified as the area across the $u$- and $w$-arms 
-of a unique triad, and we notate these areas, similarly to the triad slopes, 
-as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, 
-where $e.g.$ in \autoref{eq:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, 
+Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:qcells}) with the quarter cell that is
+the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ $u$-cell and the $i,k+k_p$ $w$-cell.
+Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:i13} and \autoref{eq:i31} in this notation,
+we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$.
+Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to
+calculate the lateral flux along its $u$-arm, at $(i+i_p,k)$,
+and then again as an $s'$ to calculate the vertical flux along its $w$-arm at $(i,k+k_p)$.
+Each vertical area $a_i$ used to calculate the lateral flux and horizontal area $a'_i$ used to
+calculate the vertical flux can also be identified as the area across the $u$- and $w$-arms of a unique triad,
+and we notate these areas, similarly to the triad slopes,
+as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$,
+where $e.g.$ in \autoref{eq:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,
 and in \autoref{eq:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.
 
 \subsection{Full triad fluxes}
-A key property of iso-neutral diffusion is that it should not affect
-the (locally referenced) density. In particular there should be no
-lateral or vertical density flux. The lateral density flux disappears so long as the
-area-integrated lateral diffusive flux from tracer cell $i,k$ to
-$i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the
-form
+A key property of iso-neutral diffusion is that it should not affect the (locally referenced) density.
+In particular there should be no lateral or vertical density flux.
+The lateral density flux disappears so long as the area-integrated lateral diffusive flux from
+tracer cell $i,k$ to $i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the form
 \begin{equation}
   \label{eq:i11}
@@ -295 +277 @@
   \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}},
 \end{equation}
-where the areas $a_i$ are as in \autoref{eq:i13}. In this case,
-separating the total lateral flux, the sum of \autoref{eq:i13} and
-\autoref{eq:i11}, into triad components, a lateral tracer
-flux
+where the areas $a_i$ are as in \autoref{eq:i13}.
+In this case, separating the total lateral flux, the sum of \autoref{eq:i13} and \autoref{eq:i11},
+into triad components, a lateral tracer flux
 \begin{equation}
   \label{eq:latflux-triad}
@@ -308 +289 @@
   \right)
 \end{equation}
-can be identified with each triad. Then, because the
-same metric factors ${e_{3w}}_{\,i}^{\,k+k_p}$ and
-${e_{1u}}_{\,i+i_p}^{\,k}$ are employed for both the density gradients
-in $ _i^k \mathbb{R}_{i_p}^{k_p}$ and the tracer gradients, the lateral
-density flux associated with each triad separately disappears.
+can be identified with each triad.
+Then, because the same metric factors ${e_{3w}}_{\,i}^{\,k+k_p}$ and ${e_{1u}}_{\,i+i_p}^{\,k}$ are employed for both
+the density gradients in $ _i^k \mathbb{R}_{i_p}^{k_p}$ and the tracer gradients,
+the lateral density flux associated with each triad separately disappears.
 \begin{equation}
   \label{eq:latflux-rho}
   {\mathbb{F}_u}_{i_p}^{k_p} (\rho)=-\alpha _i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (S)=0
 \end{equation}
-Thus the total flux $\left( F_u^{31} \right) ^i _{i,k+\frac{1}{2}} +
-\left( F_u^{11} \right) ^i _{i,k+\frac{1}{2}}$ from tracer cell $i,k$
-to $i+1,k$ must also vanish since it is a sum of four such triad fluxes.
-
-The squared slope $r_1^2$ in the expression \autoref{eq:i33c} for the
-$_{33}$ component is also expressed in terms of area-weighted
-squared triad slopes, so the area-integrated vertical flux from tracer
-cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is
+Thus the total flux $\left( F_u^{31} \right) ^i _{i,k+\frac{1}{2}} + \left( F_u^{11} \right) ^i _{i,k+\frac{1}{2}}$ from
+tracer cell $i,k$ to $i+1,k$ must also vanish since it is a sum of four such triad fluxes.
+
+The squared slope $r_1^2$ in the expression \autoref{eq:i33c} for the $_{33}$ component is also expressed in
+terms of area-weighted squared triad slopes,
+so the area-integrated vertical flux from tracer cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is
 \begin{equation}
   \label{eq:i33}
@@ -333 +311 @@
     + \Alts_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right],
 \end{equation}
-where the areas $a'$ and slopes $s'$ are the same as in
-\autoref{eq:i31}.
-Then, separating the total vertical flux, the sum of \autoref{eq:i31} and
-\autoref{eq:i33}, into triad components,  a vertical flux
+where the areas $a'$ and slopes $s'$ are the same as in \autoref{eq:i31}.
+Then, separating the total vertical flux, the sum of \autoref{eq:i31} and \autoref{eq:i33},
+into triad components, a vertical flux
 \begin{align}
   \label{eq:vertflux-triad}
@@ -349 +326 @@
    {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2}
 \end{align}
-may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$
-associated with a triad then separately disappears (because the
-lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (\rho)$
-disappears). Consequently the total vertical density flux $\left( F_w^{31} \right)_i ^{k+\frac{1}{2}} +
-\left( F_w^{33} \right)_i^{k+\frac{1}{2}}$ from tracer cell $i,k$
-to $i,k+1$ must also vanish since it is a sum of four such triad
-fluxes.
-
-We can explicitly identify (\autoref{fig:qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of
-the $u$-fluxes and $w$-fluxes in
-\autoref{eq:i31}, \autoref{eq:i13}, \autoref{eq:i11} \autoref{eq:i33} and
-\autoref{fig:ISO_triad} to  write out the iso-neutral fluxes at $u$- and
-$w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces:
+may be associated with each triad.
+Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$ associated with a triad then
+separately disappears (because the lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (\rho)$ disappears).
+Consequently the total vertical density flux
+$\left( F_w^{31} \right)_i ^{k+\frac{1}{2}} + \left( F_w^{33} \right)_i^{k+\frac{1}{2}}$ from
+tracer cell $i,k$ to $i,k+1$ must also vanish since it is a sum of four such triad fluxes.
+
+We can explicitly identify (\autoref{fig:qcells}) the triads associated with the $s_i$, $a_i$,
+and $s'_i$, $a'_i$ used in the definition of the $u$-fluxes and $w$-fluxes in \autoref{eq:i31},
+\autoref{eq:i13}, \autoref{eq:i11} \autoref{eq:i33} and \autoref{fig:ISO_triad} to write out
+the iso-neutral fluxes at $u$- and $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces:
 %(\autoref{fig:ISO_triad}):
 \begin{flalign} \label{eq:iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv
@@ -375 +350 @@
 \label{subsec:variance}
 
-We now require that this operator should not increase the
-globally-integrated tracer variance.
+We now require that this operator should not increase the globally-integrated tracer variance.
 %This changes according to
 % \begin{align*}
@@ -387 +361 @@
 %              + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\
 % \end{align*}
-Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux
-$_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the $u$-point $i+i_p,k$ and
-a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the
-$w$-point $i,k+k_p$.  The lateral flux drives a net rate of change of
-variance, summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
+Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across
+the $u$-point $i+i_p,k$ and a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the $w$-point $i,k+k_p$.
+The lateral flux drives a net rate of change of variance,
+summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
 \begin{multline}
   {b_T}_{i+i_p-1/2}^k\left(\frac{\partial T}{\partial t}T\right)_{i+i_p-1/2}^k+
@@ -402 +375 @@
  \end{aligned}
 \end{multline}
-while the vertical flux similarly drives a net rate of change of
-variance summed over the $T$-points $i,k+k_p-\half$ (above) and
-$i,k+k_p+\half$ (below) of
+while the vertical flux similarly drives a net rate of change of variance summed over
+the $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of
 \begin{equation}
 \label{eq:dvar_iso_k}
   _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
 \end{equation}
-The total variance tendency driven by the triad is the sum of these
-two. Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and
-$_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \autoref{eq:latflux-triad} and
-\autoref{eq:vertflux-triad}, it is
+The total variance tendency driven by the triad is the sum of these two.
+Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with
+\autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad}, it is
 \begin{multline*}
 -\Alts_i^k\left \{
@@ -428 +399 @@
 \right \}.
 \end{multline*}
-The key point is then that if we require
-$_i^k{\mathbb{A}_u}_{i_p}^{k_p}$ and $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$
-to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by
+The key point is then that if we require $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$ and $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$ to
+be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by
 \begin{equation}
   \label{eq:V-A}
@@ -447 +417 @@
   \right)^2\leq 0.
 \end{equation}
-Thus, the constraint \autoref{eq:V-A} ensures that the fluxes (\autoref{eq:latflux-triad}, \autoref{eq:vertflux-triad}) associated
-with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase
-the net variance. Since the total fluxes are sums of such fluxes from
-the various triads, this constraint, applied to all triads, is
-sufficient to ensure that the globally integrated variance does not
-increase.
-
-The expression \autoref{eq:V-A} can be interpreted as a discretization
-of the global integral
+Thus, the constraint \autoref{eq:V-A} ensures that the fluxes
+(\autoref{eq:latflux-triad}, \autoref{eq:vertflux-triad}) associated with
+a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase the net variance.
+Since the total fluxes are sums of such fluxes from the various triads, this constraint, applied to all triads,
+is sufficient to ensure that the globally integrated variance does not increase.
+
+The expression \autoref{eq:V-A} can be interpreted as a discretization of the global integral
 \begin{equation}
   \label{eq:cts-var}
@@ -461 +429 @@
   \int\!\mathbf{F}\cdot\nabla T\, dV,
 \end{equation}
-where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the
-lateral and vertical fluxes/unit area
+where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the lateral and vertical fluxes/unit area
 \[
 \mathbf{F}=\left(
@@ -477 +444 @@
 
 \subsection{Triad volumes in Griffes's scheme and in \NEMO}
-To complete the discretization we now need only specify the triad
-volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. \citet{Griffies_al_JPO98} identify
-these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter
-cells, defined in terms of the distances between $T$, $u$,$f$ and
-$w$-points. This is the natural discretization of
-\autoref{eq:cts-var}. The \NEMO model, however, operates with scale
-factors instead of grid sizes, and scale factors for the quarter
-cells are not defined. Instead, therefore we simply choose
+To complete the discretization we now need only specify the triad volumes $_i^k\mathbb{V}_{i_p}^{k_p}$.
+\citet{Griffies_al_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells,
+defined in terms of the distances between $T$, $u$,$f$ and $w$-points.
+This is the natural discretization of \autoref{eq:cts-var}.
+The \NEMO model, however, operates with scale factors instead of grid sizes,
+and scale factors for the quarter cells are not defined.
+Instead, therefore we simply choose
 \begin{equation}
   \label{eq:V-NEMO}
   _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k,
 \end{equation}
-as a quarter of the volume of the $u$-cell inside which the triad
-quarter-cell lies. This has the nice property that when the slopes
-$\mathbb{R}$ vanish, the lateral flux from tracer cell $i,k$ to
-$i+1,k$ reduces to the classical form
+as a quarter of the volume of the $u$-cell inside which the triad quarter-cell lies.
+This has the nice property that when the slopes $\mathbb{R}$ vanish,
+the lateral flux from tracer cell $i,k$ to $i+1,k$ reduces to the classical form
 \begin{equation}
   \label{eq:lat-normal}
@@ -500 +465 @@
  = -\overline\Alts_{\,i+1/2}^k\;\frac{{e_{1w}}_{\,i + 1/2}^{\,k}\:{e_{1v}}_{\,i + 1/2}^{\,k}\;\delta_{i+ 1/2}[T^k]}{{e_{1u}}_{\,i + 1/2}^{\,k}}.
 \end{equation}
-In fact if the diffusive coefficient is defined at $u$-points, so that
-we employ $\Alts_{i+i_p}^k$ instead of  $\Alts_i^k$ in the definitions of the
-triad fluxes \autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad},
+In fact if the diffusive coefficient is defined at $u$-points,
+so that we employ $\Alts_{i+i_p}^k$ instead of  $\Alts_i^k$ in the definitions of the triad fluxes
+\autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad},
 we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above.
 
 \subsection{Summary of the scheme}
-The iso-neutral fluxes at $u$- and
-$w$-points are the sums of the triad fluxes that cross the $u$- and
-$w$-faces \autoref{eq:iso_flux}:
+The iso-neutral fluxes at $u$- and $w$-points are the sums of the triad fluxes that
+cross the $u$- and $w$-faces \autoref{eq:iso_flux}:
 \begin{subequations}\label{eq:alltriadflux}
   \begin{flalign}\label{eq:vect_isoflux}
@@ -545 +509 @@
 \end{subequations}
 
- The divergence of the expression \autoref{eq:iso_flux} for the fluxes gives the iso-neutral diffusion tendency at
+The divergence of the expression \autoref{eq:iso_flux} for the fluxes gives the iso-neutral diffusion tendency at
 each tracer point:
 \begin{equation} \label{eq:iso_operator} D_l^T = \frac{1}{b_T}
@@ -555 +519 @@
 The diffusion scheme satisfies the following six properties:
 \begin{description}
-\item[$\bullet$ horizontal diffusion] The discretization of the
-  diffusion operator recovers \autoref{eq:lat-normal} the traditional five-point Laplacian in
-  the limit of flat iso-neutral direction :
+\item[$\bullet$ horizontal diffusion]
+  The discretization of the diffusion operator recovers the traditional five-point Laplacian
+  \autoref{eq:lat-normal} in the limit of flat iso-neutral direction:
   \begin{equation} \label{eq:iso_property0} D_l^T = \frac{1}{b_T} \
     \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \;
@@ -564 +528 @@
   \end{equation}
 
-\item[$\bullet$ implicit treatment in the vertical] Only tracer values
-  associated with a single water column appear in the expression
-  \autoref{eq:i33} for the $_{33}$ fluxes, vertical fluxes driven by
-  vertical gradients. This is of paramount importance since it means
-  that a time-implicit algorithm can be used to solve the vertical
-  diffusion equation. This is necessary
- since the vertical eddy
-  diffusivity associated with this term,
+\item[$\bullet$ implicit treatment in the vertical]
+  Only tracer values associated with a single water column appear in the expression \autoref{eq:i33} for
+  the $_{33}$ fluxes, vertical fluxes driven by vertical gradients.
+  This is of paramount importance since it means that a time-implicit algorithm can be used to
+  solve the vertical diffusion equation.
+  This is necessary since the vertical eddy diffusivity associated with this term,
   \begin{equation}
     \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{
@@ -579 +541 @@
       {b_u}_{i+i_p}^k\: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
     \right\},
- \end{equation}
+  \end{equation}
   (where $b_w= e_{1w}\,e_{2w}\,e_{3w}$ is the volume of $w$-cells) can be quite large.
 
-\item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of
-  locally referenced potential density is zero. See
-  \autoref{eq:latflux-rho} and \autoref{eq:vertflux-triad2}.
-
-\item[$\bullet$ conservation of tracer] The iso-neutral diffusion
-  conserves tracer content, $i.e.$
+\item[$\bullet$ pure iso-neutral operator]
+  The iso-neutral flux of locally referenced potential density is zero.
+  See \autoref{eq:latflux-rho} and \autoref{eq:vertflux-triad2}.
+
+\item[$\bullet$ conservation of tracer]
+  The iso-neutral diffusion conserves tracer content, $i.e.$
   \begin{equation} \label{eq:iso_property1} \sum_{i,j,k} \left\{ D_l^T \
       b_T \right\} = 0
   \end{equation}
-  This property is trivially satisfied since the iso-neutral diffusive
-  operator is written in flux form.
-
-\item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion
-  does not increase the tracer variance, $i.e.$
+  This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form.
+
+\item[$\bullet$ no increase of tracer variance]
+  The iso-neutral diffusion does not increase the tracer variance, $i.e.$
   \begin{equation} \label{eq:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T
       \ b_T \right\} \leq 0
   \end{equation}
-  The property is demonstrated in
-  \autoref{subsec:variance} above. It is a key property for a diffusion
-  term. It means that it is also a dissipation term, $i.e.$ it
-  dissipates the square of the quantity on which it is applied.  It
-  therefore ensures that, when the diffusivity coefficient is large
-  enough, the field on which it is applied becomes free of grid-point
-  noise.
-
-\item[$\bullet$ self-adjoint operator] The iso-neutral diffusion
-  operator is self-adjoint, $i.e.$
+  The property is demonstrated in \autoref{subsec:variance} above.
+  It is a key property for a diffusion term.
+  It means that it is also a dissipation term,
+  $i.e.$ it dissipates the square of the quantity on which it is applied.
+  It therefore ensures that, when the diffusivity coefficient is large enough,
+  the field on which it is applied becomes free of grid-point noise.
+
+\item[$\bullet$ self-adjoint operator]
+  The iso-neutral diffusion operator is self-adjoint, $i.e.$
   \begin{equation} \label{eq:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T
       \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}
   \end{equation}
-  In other word, there is no need to develop a specific routine from
-  the adjoint of this operator. We just have to apply the same
-  routine. This property can be demonstrated similarly to the proof of
-  the `no increase of tracer variance' property. The contribution by a
-  single triad towards the left hand side of \autoref{eq:iso_property3}, can
-  be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:dvar_iso_i}
-  and \autoref{eq:dvar_iso_k}. This results in a term similar to
-  \autoref{eq:perfect-square},
+  In other word, there is no need to develop a specific routine from the adjoint of this operator.
+  We just have to apply the same routine.
+  This property can be demonstrated similarly to the proof of the `no increase of tracer variance' property.
+  The contribution by a single triad towards the left hand side of \autoref{eq:iso_property3},
+  can be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:dvar_iso_i} and \autoref{eq:dvar_iso_k}.
+  This results in a term similar to \autoref{eq:perfect-square},
 \begin{equation}
   \label{eq:TScovar}
@@ -634 +592 @@
   \right).
 \end{equation}
-This is symmetrical in $T $ and $S$, so exactly the same term arises
-from the discretization of this triad's contribution towards the
-RHS of \autoref{eq:iso_property3}.
+This is symmetrical in $T $ and $S$, so exactly the same term arises from
+the discretization of this triad's contribution towards the RHS of \autoref{eq:iso_property3}.
 \end{description}
 
 \subsection{Treatment of the triads at the boundaries}\label{sec:iso_bdry}
-The triad slope can only be defined where both the grid boxes centred at
-the end of the arms exist. Triads that would poke up
-through the upper ocean surface into the atmosphere, or down into the
-ocean floor, must be masked out. See \autoref{fig:bdry_triads}. Surface layer triads
-$\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and
-$\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be
-specified above the ocean surface are masked (\autoref{fig:bdry_triads}a): this ensures that lateral
-tracer gradients produce no flux through the ocean surface. However,
-to prevent surface noise, it is customary to retain the $_{11}$ contributions towards
-the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and
-$\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer
-fluxes. Similar comments apply to triads that would intersect the
-ocean floor (\autoref{fig:bdry_triads}b). Note that both near bottom
-triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and
-$\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$
-or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is
-masked. The associated lateral fluxes (grey-black dashed line) are
-masked if \np{ln\_botmix\_triad}\forcode{ = .false.}, but left unmasked,
-giving bottom mixing, if \np{ln\_botmix\_triad}\forcode{ = .true.}.
-
-The default option \np{ln\_botmix\_triad}\forcode{ = .false.} is suitable when the
-bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}\forcode{ = 1}),
-or  for simple idealized  problems. For setups with topography without
-bbl mixing, \np{ln\_botmix\_triad}\forcode{ = .true.} may be necessary.
+The triad slope can only be defined where both the grid boxes centred at the end of the arms exist.
+Triads that would poke up through the upper ocean surface into the atmosphere,
+or down into the ocean floor, must be masked out.
+See \autoref{fig:bdry_triads}.
+Surface layer triads $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that
+require density to be specified above the ocean surface are masked (\autoref{fig:bdry_triads}a):
+this ensures that lateral tracer gradients produce no flux through the ocean surface.
+However, to prevent surface noise, it is customary to retain the $_{11}$ contributions towards
+the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$;
+this drives diapycnal tracer fluxes.
+Similar comments apply to triads that would intersect the ocean floor (\autoref{fig:bdry_triads}b).
+Note that both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and
+$\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ or $i+1,k+1$ tracer points is masked,
+i.e.\ the $i,k+1$ $u$-point is masked.
+The associated lateral fluxes (grey-black dashed line) are masked if \np{ln\_botmix\_triad}\forcode{ = .false.},
+but left unmasked, giving bottom mixing, if \np{ln\_botmix\_triad}\forcode{ = .true.}.
+
+The default option \np{ln\_botmix\_triad}\forcode{ = .false.} is suitable when the bbl mixing option is enabled
+(\key{trabbl}, with \np{nn\_bbl\_ldf}\forcode{ = 1}), or for simple idealized problems.
+For setups with topography without bbl mixing, \np{ln\_botmix\_triad}\forcode{ = .true.} may be necessary.
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[h] \begin{center}
     \includegraphics[width=0.60\textwidth]{Fig_GRIFF_bdry_triads}
     \caption{  \protect\label{fig:bdry_triads}
-      (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer
-      points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad
-      slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$
-      (blue) poking through the ocean surface are masked (faded in
-      figure). However, the lateral $_{11}$ contributions towards
-      $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$
-      (yellow line) are still applied, giving diapycnal diffusive
-      fluxes.\newline
+      (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer points (black dots),
+      and $i+1/2,1$ $u$-point (blue square).
+      Triad slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) poking through
+      the ocean surface are masked (faded in figure).
+      However, the lateral $_{11}$ contributions towards $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and
+      $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$ (yellow line) are still applied,
+      giving diapycnal diffusive fluxes.\newline
       (b) Both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and
-      $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$
-      or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point
-      is masked. The associated lateral fluxes (grey-black dashed
-      line) are masked if \protect\np{botmix\_triad}\forcode{ = .false.}, but left
-      unmasked, giving bottom mixing, if \protect\np{botmix\_triad}\forcode{ = .true.}}
+      $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ or $i+1,k+1$ tracer points is masked,
+      i.e.\ the $i,k+1$ $u$-point is masked.
+      The associated lateral fluxes (grey-black dashed line) are masked if
+      \protect\np{botmix\_triad}\forcode{ = .false.}, but left unmasked,
+      giving bottom mixing, if \protect\np{botmix\_triad}\forcode{ = .true.}}
  \end{center} \end{figure}
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 \subsection{ Limiting of the slopes within the interior}\label{sec:limit}
-As discussed in \autoref{subsec:LDF_slp_iso}, iso-neutral slopes relative to
-geopotentials must be bounded everywhere, both for consistency with the small-slope
-approximation and for numerical stability \citep{Cox1987,
-  Griffies_Bk04}. The bound chosen in \NEMO is applied to each
-component of the slope separately and has a value of $1/100$ in the ocean interior.
+As discussed in \autoref{subsec:LDF_slp_iso},
+iso-neutral slopes relative to geopotentials must be bounded everywhere,
+both for consistency with the small-slope approximation and for numerical stability \citep{Cox1987, Griffies_Bk04}.
+The bound chosen in \NEMO is applied to each component of the slope separately and
+has a value of $1/100$ in the ocean interior.
 %, ramping linearly down above 70~m depth to zero at the surface
-It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to
-geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to
-geopotentials) \autoref{eq:PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate
-surfaces, so we require
+It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to geopotentials
+(here the $\sigma_i$ are the slopes of the coordinate surfaces relative to geopotentials)
+\autoref{eq:PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate surfaces, so we require
 \begin{equation*}
   |\tilde{r}_i|\leq \tilde{r}_\mathrm{max}=0.01.
@@ -701 +652 @@
 and then recalculate the slopes $r_i$ relative to coordinates.
 Each individual triad slope
- \begin{equation}
-   \label{eq:Rtilde}
-_i^k\tilde{\mathbb{R}}_{i_p}^{k_p} = {}_i^k\mathbb{R}_{i_p}^{k_p}  + \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}
- \end{equation}
-is limited like this and then the corresponding
-$_i^k\mathbb{R}_{i_p}^{k_p} $ are recalculated and combined to form the fluxes.
-Note that where the slopes have been limited, there is now a non-zero
-iso-neutral density flux that drives dianeutral mixing.  In particular this iso-neutral density flux
-is always downwards, and so acts to reduce gravitational potential energy.
+\begin{equation}
+  \label{eq:Rtilde}
+  _i^k\tilde{\mathbb{R}}_{i_p}^{k_p} = {}_i^k\mathbb{R}_{i_p}^{k_p}  + \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}
+\end{equation}
+is limited like this and then the corresponding $_i^k\mathbb{R}_{i_p}^{k_p} $ are recalculated and
+combined to form the fluxes.
+Note that where the slopes have been limited, there is now a non-zero iso-neutral density flux that
+drives dianeutral mixing.
+In particular this iso-neutral density flux is always downwards,
+and so acts to reduce gravitational potential energy.
 
 \subsection{Tapering within the surface mixed layer}\label{sec:taper}
-Additional tapering of the iso-neutral fluxes is necessary within the
-surface mixed layer. When the Griffies triads are used, we offer two
-options for this.
+Additional tapering of the iso-neutral fluxes is necessary within the surface mixed layer.
+When the Griffies triads are used, we offer two options for this.
 
 \subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:lintaper}
-This is the option activated by the default choice
-\np{ln\_triad\_iso}\forcode{ = .false.}. Slopes $\tilde{r}_i$ relative to
-geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the
-surface, as described in option (c) of \autoref{fig:eiv_slp}, to values
+This is the option activated by the default choice \np{ln\_triad\_iso}\forcode{ = .false.}.
+Slopes $\tilde{r}_i$ relative to geopotentials are tapered linearly from their value immediately below
+the mixed layer to zero at the surface, as described in option (c) of \autoref{fig:eiv_slp}, to values
 \begin{subequations}
   \begin{equation}
-   \label{eq:rmtilde}
-     \rMLt =
-  -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h,
+    \label{eq:rmtilde}
+    \rMLt =
+    -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h,
   \end{equation}
-and then the $r_i$ relative to vertical coordinate surfaces are appropriately
-adjusted to
+  and then the $r_i$ relative to vertical coordinate surfaces are appropriately adjusted to
   \begin{equation}
-   \label{eq:rm}
- \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h.
+    \label{eq:rm}
+    \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h.
   \end{equation}
 \end{subequations}
@@ -744 +693 @@
 \end{equation}
 
-This slope tapering gives a natural connection between tracer in the
-mixed-layer and in isopycnal layers immediately below, in the
-thermocline. It is consistent with the way the $\tilde{r}_i$ are
-tapered within the mixed layer (see \autoref{sec:taperskew} below)
-so as to ensure a uniform GM eddy-induced velocity throughout the
-mixed layer. However, it gives a downwards density flux and so acts so
-as to reduce potential energy in the same way as does the slope
-limiting discussed above in \autoref{sec:limit}.
+This slope tapering gives a natural connection between tracer in the mixed-layer and
+in isopycnal layers immediately below, in the thermocline.
+It is consistent with the way the $\tilde{r}_i$ are tapered within the mixed layer
+(see \autoref{sec:taperskew} below) so as to ensure a uniform GM eddy-induced velocity throughout the mixed layer.
+However, it gives a downwards density flux and so acts so as to reduce potential energy in the same way as
+does the slope limiting discussed above in \autoref{sec:limit}.
  
-As in \autoref{sec:limit} above, the tapering
-\autoref{eq:rmtilde} is applied separately to each triad
-$_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the
-$_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume
-$z$-coordinates in the following; the conversion from
-$\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described
-above by \autoref{eq:Rtilde}.
+As in \autoref{sec:limit} above, the tapering \autoref{eq:rmtilde} is applied separately to
+each triad $_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the $_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted.
+For clarity, we assume $z$-coordinates in the following;
+the conversion from $\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as
+described above by \autoref{eq:Rtilde}.
 \begin{enumerate}
-\item Mixed-layer depth is defined so as to avoid including regions of weak
-vertical stratification in the slope definition.
- At each $i,j$ (simplified to $i$ in
-\autoref{fig:MLB_triad}), we define the mixed-layer by setting
-the vertical index of the tracer point immediately below the mixed
-layer, $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point)
-such that the potential density
-${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, where $i,k_{10}$ is
-the tracer gridbox within which the depth reaches 10~m. See the left
-side of \autoref{fig:MLB_triad}. We use the $k_{10}$-gridbox
-instead of the surface gridbox to avoid problems e.g.\ with thin
-daytime mixed-layers. Currently we use the same
-$\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is
-used to output the diagnosed mixed-layer depth
-$h_{\mathrm{ML}}=|z_{W}|_{k_{\mathrm{ML}}+1/2}$, the depth of the $w$-point
-above the $i,k_{\mathrm{ML}}$ tracer point.
-
-\item We define `basal' triad slopes
-${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ as the slopes
-of those triads whose vertical `arms' go down from the
-$i,k_{\mathrm{ML}}$ tracer point to the $i,k_{\mathrm{ML}}-1$ tracer point
-below. This is to ensure that the vertical density gradients
-associated with these basal triad slopes
-${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are
-representative of the thermocline. The four basal triads defined in the bottom part
-of \autoref{fig:MLB_triad} are then
+\item
+  Mixed-layer depth is defined so as to avoid including regions of weak vertical stratification in
+  the slope definition.
+  At each $i,j$ (simplified to $i$ in \autoref{fig:MLB_triad}),
+  we define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer,
+  $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) such that
+  the potential density ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$,
+  where $i,k_{10}$ is the tracer gridbox within which the depth reaches 10~m.
+  See the left side of \autoref{fig:MLB_triad}.
+  We use the $k_{10}$-gridbox instead of the surface gridbox to avoid problems e.g.\ with thin daytime mixed-layers.
+  Currently we use the same $\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is used to
+  output the diagnosed mixed-layer depth $h_{\mathrm{ML}}=|z_{W}|_{k_{\mathrm{ML}}+1/2}$,
+  the depth of the $w$-point above the $i,k_{\mathrm{ML}}$ tracer point.
+\item
+  We define `basal' triad slopes ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ as
+  the slopes of those triads whose vertical `arms' go down from the $i,k_{\mathrm{ML}}$ tracer point to
+  the $i,k_{\mathrm{ML}}-1$ tracer point below.
+  This is to ensure that the vertical density gradients associated with
+  these basal triad slopes ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are representative of the thermocline.
+  The four basal triads defined in the bottom part of \autoref{fig:MLB_triad} are then
 \begin{align}
   {\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p} &=
@@ -795 +734 @@
 {\:}^{k_{\mathrm{ML}}}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2}. \notag
 \end{align}
-The vertical flux associated with each of these triads passes through the $w$-point
-$i,k_{\mathrm{ML}}-1/2$ lying \emph{below} the $i,k_{\mathrm{ML}}$ tracer point,
-so it is this depth
+The vertical flux associated with each of these triads passes through
+the $w$-point $i,k_{\mathrm{ML}}-1/2$ lying \emph{below} the $i,k_{\mathrm{ML}}$ tracer point, so it is this depth
 \begin{equation}
   \label{eq:zbase}
   {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2}
 \end{equation}
-(one gridbox deeper than the
-diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper
-the slopes in \autoref{eq:rmtilde}.
-\item Finally, we calculate the adjusted triads
-${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within the mixed
-layer, by multiplying the appropriate
-${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ by the ratio of
-the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_{\mathrm{base}}}_{\,i}$. For
-instance the green triad centred on $i,k$
+one gridbox deeper than the diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper the slopes in
+\autoref{eq:rmtilde}.
+\item
+  Finally, we calculate the adjusted triads ${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within
+  the mixed layer, by multiplying the appropriate ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ by
+  the ratio of the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_{\mathrm{base}}}_{\,i}$.
+  For instance the green triad centred on $i,k$
 \begin{align}
   {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,1/2}^{-1/2} &=
@@ -824 +760 @@
 \begin{figure}[h]
 %  \fcapside {
-    \caption{\protect\label{fig:MLB_triad} Definition of
-      mixed-layer depth and calculation of linearly tapered
-      triads. The figure shows a water column at a given $i,j$
-      (simplified to $i$), with the ocean surface at the top. Tracer points are denoted by
-      bullets, and black lines the edges of the tracer cells; $k$
-      increases upwards. \newline
-      \hspace{5 em}We define the mixed-layer by setting the vertical index
-      of the tracer point immediately below the mixed layer,
-      $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point)
-      such that ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$,
-      where $i,k_{10}$ is the tracer gridbox within which the depth
-      reaches 10~m. We calculate the triad slopes within the mixed
-      layer by linearly tapering them from zero (at the surface) to
-      the `basal' slopes, the slopes of the four triads passing through the
-      $w$-point $i,k_{\mathrm{ML}}-1/2$ (blue square),
-      ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$. Triads with
-    different $i_p,k_p$, denoted by different colours, (e.g. the green
-    triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.}
+  \caption{\protect\label{fig:MLB_triad}
+    Definition of mixed-layer depth and calculation of linearly tapered triads.
+    The figure shows a water column at a given $i,j$ (simplified to $i$), with the ocean surface at the top.
+    Tracer points are denoted by bullets, and black lines the edges of the tracer cells;
+    $k$ increases upwards. \newline
+    \hspace{5 em}
+    We define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer,
+    $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) such that
+    ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$,
+    where $i,k_{10}$ is the tracer gridbox within which the depth reaches 10~m.
+    We calculate the triad slopes within the mixed layer by linearly tapering them from zero
+    (at the surface) to the `basal' slopes,
+    the slopes of the four triads passing through the $w$-point $i,k_{\mathrm{ML}}-1/2$ (blue square),
+    ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$.
+    Triads with different $i_p,k_p$, denoted by different colours,
+    (e.g. the green triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.}
 %}
   {\includegraphics[width=0.60\textwidth]{Fig_GRIFF_MLB_triads}}
@@ -849 +783 @@
 \subsubsection{Additional truncation of skew iso-neutral flux components}
 \label{subsec:Gerdes-taper}
-The alternative option is activated by setting \np{ln\_triad\_iso} =
-  true. This retains the same tapered slope $\rML$  described above for the
-calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the
-vertical tracer flux driven by vertical tracer gradients), but
-replaces the $\rML$ in the skew term by
+The alternative option is activated by setting \np{ln\_triad\_iso} = true.
+This retains the same tapered slope $\rML$  described above for the calculation of the $_{33}$ term of
+the iso-neutral diffusion tensor (the vertical tracer flux driven by vertical tracer gradients),
+but replaces the $\rML$ in the skew term by
 \begin{equation}
   \label{eq:rm*}
@@ -868 +801 @@
 \end{equation}
 This operator
-\footnote{To ensure good behaviour where horizontal density
-  gradients are weak, we in fact follow \citet{Gerdes1991} and set
-$\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.}
-then has the property it gives no vertical density flux, and so does
-not change the potential energy.
-This approach is similar to multiplying the iso-neutral  diffusion
-coefficient by $\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep
-slopes, as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}).
+\footnote{
+  To ensure good behaviour where horizontal density gradients are weak,
+  we in fact follow \citet{Gerdes1991} and
+  set $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.}
+then has the property it gives no vertical density flux, and so does not change the potential energy.
+This approach is similar to multiplying the iso-neutral diffusion coefficient by
+$\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep slopes,
+as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}).
 Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$
 
-In practice, this approach gives weak vertical tracer fluxes through
-the mixed-layer, as well as vanishing density fluxes. While it is
-theoretically advantageous that it does not change the potential
-energy, it may give a discontinuity between the
-fluxes within the mixed-layer (purely horizontal) and just below (along
-iso-neutral surfaces).
+In practice, this approach gives weak vertical tracer fluxes through the mixed-layer,
+as well as vanishing density fluxes.
+While it is theoretically advantageous that it does not change the potential energy,
+it may give a discontinuity between the fluxes within the mixed-layer (purely horizontal) and
+just below (along iso-neutral surfaces).
 % This may give strange looking results,
 % particularly where the mixed-layer depth varies strongly laterally.
@@ -893 +825 @@
 \subsection{Continuous skew flux formulation}\label{sec:continuous-skew-flux}
 
- When Gent and McWilliams's [1990] diffusion is used,
-an additional advection term is added. The associated velocity is the so called
-eddy induced velocity, 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.$ \autoref{eq:ldfslp_geo}
-is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo}
-+ \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates.
+When Gent and McWilliams's [1990] diffusion is used, an additional advection term is added.
+The associated velocity is the so called eddy induced velocity,
+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.$ \autoref{eq:ldfslp_geo} is used in $z$-coordinate,
+and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates.
 
 The eddy induced velocity is given by:
@@ -919 +850 @@
 \end{equation}
 \end{subequations}
-with $A_{e}$ the eddy induced velocity coefficient, and $\tilde{r}_1$ and $\tilde{r}_2$ the slopes between the iso-neutral and the geopotential surfaces.
-
-The traditional way to implement this additional advection is to add
-it to the Eulerian velocity prior to computing the tracer
-advection. This is implemented if \key{traldf\_eiv} is set in the
-default implementation, where \np{ln\_traldf\_triad} is set
-false. This allows us to take advantage of all the advection schemes
-offered for the tracers (see \autoref{sec:TRA_adv}) and not just a $2^{nd}$
-order advection scheme. This is particularly useful for passive
-tracers where \emph{positivity} of the advection scheme is of
-paramount importance.
-
-However, when \np{ln\_traldf\_triad} is set true, \NEMO instead
-implements eddy induced advection according to the so-called skew form
-\citep{Griffies_JPO98}. It is based on a transformation of the advective fluxes
-using the non-divergent nature of the eddy induced velocity.
-For example in the (\textbf{i},\textbf{k}) plane, the tracer advective
-fluxes per unit area in $ijk$ space can be
-transformed as follows:
+with $A_{e}$ the eddy induced velocity coefficient,
+and $\tilde{r}_1$ and $\tilde{r}_2$ the slopes between the iso-neutral and the geopotential surfaces.
+
+The traditional way to implement this additional advection is to add it to the Eulerian velocity prior to
+computing the tracer advection.
+This is implemented if \key{traldf\_eiv} is set in the default implementation,
+where \np{ln\_traldf\_triad} is set false.
+This allows us to take advantage of all the advection schemes offered for the tracers
+(see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme.
+This is particularly useful for passive tracers where
+\emph{positivity} of the advection scheme is of paramount importance.
+
+However, when \np{ln\_traldf\_triad} is set true,
+\NEMO instead implements eddy induced advection according to the so-called skew form \citep{Griffies_JPO98}.
+It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity.
+For example in the (\textbf{i},\textbf{k}) plane,
+the tracer advective fluxes per unit area in $ijk$ space can be transformed as follows:
 \begin{flalign*}
 \begin{split}
@@ -962 +891 @@
 \end{split}
 \end{flalign*}
-and since the eddy induced velocity field is non-divergent, we end up with the skew
-form of the eddy induced advective fluxes per unit area in $ijk$ space:
+and since the eddy induced velocity field is non-divergent,
+we end up with the skew form of the eddy induced advective fluxes per unit area in $ijk$ space:
 \begin{equation} \label{eq:eiv_skew_ijk}
 \textbf{F}_\mathrm{eiv}^T = \begin{pmatrix}
@@ -979 +908 @@
 \end{split}
 \end{equation}
-Note that  \autoref{eq:eiv_skew_physical} takes the same form whatever the
-vertical coordinate, though of course the slopes
-$\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:eiv_psi} are relative to geopotentials.
-The tendency associated with eddy induced velocity is then simply the convergence
-of the fluxes (\autoref{eq:eiv_skew_ijk}, \autoref{eq:eiv_skew_physical}), so
+Note that \autoref{eq:eiv_skew_physical} takes the same form whatever the vertical coordinate,
+though of course the slopes $\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:eiv_psi} are relative to
+geopotentials.
+The tendency associated with eddy induced velocity is then simply the convergence of the fluxes
+(\autoref{eq:eiv_skew_ijk}, \autoref{eq:eiv_skew_physical}), so
 \begin{equation} \label{eq:skew_eiv_conv}
 \frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[
@@ -992 +921 @@
    + e_{1} \psi_2 \partial_j T \right)  \right]
 \end{equation}
- It naturally conserves the tracer content, as it is expressed in flux
- form. Since it has the same divergence as the advective form it also
- preserves the tracer variance.
+It naturally conserves the tracer content, as it is expressed in flux form.
+Since it has the same divergence as the advective form it also preserves the tracer variance.
 
 \subsection{Discrete skew flux formulation}
-The skew fluxes in (\autoref{eq:eiv_skew_physical}, \autoref{eq:eiv_skew_ijk}), like the off-diagonal terms
-(\autoref{eq:i13c}, \autoref{eq:i31c}) of the small angle diffusion tensor, are best
-expressed in terms of the triad slopes, as in \autoref{fig:ISO_triad}
-and (\autoref{eq:i13}, \autoref{eq:i31}); but now in terms of the triad slopes
-$\tilde{\mathbb{R}}$ relative to geopotentials instead of the
-$\mathbb{R}$ relative to coordinate surfaces. The discrete form of
-\autoref{eq:eiv_skew_ijk} using the slopes \autoref{eq:R} and
+The skew fluxes in (\autoref{eq:eiv_skew_physical}, \autoref{eq:eiv_skew_ijk}),
+like the off-diagonal terms (\autoref{eq:i13c}, \autoref{eq:i31c}) of the small angle diffusion tensor,
+are best expressed in terms of the triad slopes, as in \autoref{fig:ISO_triad} and
+(\autoref{eq:i13}, \autoref{eq:i31});
+but now in terms of the triad slopes $\tilde{\mathbb{R}}$ relative to geopotentials instead of
+the $\mathbb{R}$ relative to coordinate surfaces.
+The discrete form of \autoref{eq:eiv_skew_ijk} using the slopes \autoref{eq:R} and
 defining $A_e$ at $T$-points is then given by:
 
@@ -1017 +945 @@
     \end{pmatrix},
   \end{flalign}
-  where the skew flux in the $i$-direction associated with a given
-  triad is (\autoref{eq:latflux-triad}, \autoref{eq:triadfluxu}):
+  where the skew flux in the $i$-direction associated with a given triad is (\autoref{eq:latflux-triad},
+  \autoref{eq:triadfluxu}):
   \begin{align}
     \label{eq:skewfluxu}
@@ -1034 +962 @@
 \end{subequations}
 
-Such a discretisation is consistent with the iso-neutral
-operator as it uses the same definition for the slopes.  It also
-ensures the following two key properties.
+Such a discretisation is consistent with the iso-neutral operator as it uses the same definition for the slopes.
+It also ensures the following two key properties.
 
 \subsubsection{No change in tracer variance}
-The discretization conserves tracer variance, $i.e.$ it does not
-include a diffusive component but is a `pure' advection term. This can
-be seen
-%either from Appendix \autoref{apdx:eiv_skew} or
-by considering the
-fluxes associated with a given triad slope
-$_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. For, following
-\autoref{subsec:variance} and \autoref{eq:dvar_iso_i}, the
-associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$
-drives a net rate of change of variance, summed over the two
-$T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
+The discretization conserves tracer variance, $i.e.$ it does not include a diffusive component but is a `pure' advection term.
+This can be seen %either from Appendix \autoref{apdx:eiv_skew} or
+by considering the fluxes associated with a given triad slope $_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$.
+For, following \autoref{subsec:variance} and \autoref{eq:dvar_iso_i},
+the associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ drives a net rate of change of variance,
+summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
 \begin{equation}
 \label{eq:dvar_eiv_i}
   _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k],
 \end{equation}
-while the associated vertical skew-flux gives a variance change summed over the
-$T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of
+while the associated vertical skew-flux gives a variance change summed over
+the $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of
 \begin{equation}
 \label{eq:dvar_eiv_k}
   _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
 \end{equation}
-Inspection of the definitions (\autoref{eq:skewfluxu}, \autoref{eq:skewfluxw})
-shows that these two variance changes (\autoref{eq:dvar_eiv_i}, \autoref{eq:dvar_eiv_k})
-sum to zero. Hence the two fluxes associated with each triad make no
-net contribution to the variance budget.
+Inspection of the definitions (\autoref{eq:skewfluxu}, \autoref{eq:skewfluxw}) shows that
+these two variance changes (\autoref{eq:dvar_eiv_i}, \autoref{eq:dvar_eiv_k}) sum to zero.
+Hence the two fluxes associated with each triad make no net contribution to the variance budget.
 
 \subsubsection{Reduction in gravitational PE}
-The vertical density flux associated with the vertical skew-flux
-always has the same sign as the vertical density gradient; thus, so
-long as the fluid is stable (the vertical density gradient is
-negative) the vertical density flux is negative (downward) and hence
-reduces the gravitational PE.
+The vertical density flux associated with the vertical skew-flux always has the same sign as
+the vertical density gradient;
+thus, so long as the fluid is stable (the vertical density gradient is negative)
+the vertical density flux is negative (downward) and hence reduces the gravitational PE.
 
 For the change in gravitational PE driven by the $k$-flux is
@@ -1091 +1011 @@
 \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}},
 \end{align}
-using the definition of the triad slope $\rtriad{R}$,
-\autoref{eq:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+
-\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of  $-\alpha_i^k \delta_{k+
-  k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$.
+using the definition of the triad slope $\rtriad{R}$, \autoref{eq:R} to
+express $-\alpha _i^k\delta_{i+ i_p}[T^k]+\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of
+$-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$.
 
 Where the coordinates slope, the $i$-flux gives a PE change
@@ -1108 +1027 @@
 \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}},
 \end{multline}
-(using \autoref{eq:skewfluxu}) and so the total PE change
-\autoref{eq:vert_densityPE} + \autoref{eq:lat_densityPE} associated with the triad fluxes is
+(using \autoref{eq:skewfluxu}) and so the total PE change \autoref{eq:vert_densityPE} +
+\autoref{eq:lat_densityPE} associated with the triad fluxes is
 \begin{multline}
   \label{eq:tot_densityPE}
@@ -1122 +1041 @@
 
 \subsection{Treatment of the triads at the boundaries}\label{sec:skew_bdry}
-Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes
-are masked at the boundaries in exactly the same way as are the triad
-slopes \rtriad{R} used for the iso-neutral diffusive fluxes, as
-described in \autoref{sec:iso_bdry} and
-\autoref{fig:bdry_triads}. Thus surface layer triads
-$\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are
-masked, and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$
-and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the
-$i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$
-$u$-point is masked. The namelist parameter \np{ln\_botmix\_triad} has
-no effect on the eddy-induced skew-fluxes.
+Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes are masked at the boundaries 
+in exactly the same way as are the triad slopes \rtriad{R} used for the iso-neutral diffusive fluxes, 
+as described in \autoref{sec:iso_bdry} and \autoref{fig:bdry_triads}. 
+Thus surface layer triads $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are masked, 
+and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when 
+either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is masked. 
+The namelist parameter \np{ln\_botmix\_triad} has no effect on the eddy-induced skew-fluxes.
 
 \subsection{Limiting of the slopes within the interior}\label{sec:limitskew}
-Presently, the iso-neutral slopes $\tilde{r}_i$ relative
-to geopotentials are limited to be less than $1/100$, exactly as in
-calculating the iso-neutral diffusion, \S \autoref{sec:limit}. Each
-individual triad \rtriadt{R} is so limited.
+Presently, the iso-neutral slopes $\tilde{r}_i$ relative to geopotentials are limited to be less than $1/100$, 
+exactly as in calculating the iso-neutral diffusion, \S \autoref{sec:limit}. 
+Each individual triad \rtriadt{R} is so limited.
 
 \subsection{Tapering within the surface mixed layer}\label{sec:taperskew}
-The slopes $\tilde{r}_i$ relative to
-geopotentials (and thus the individual triads \rtriadt{R}) are always tapered linearly from their value immediately below the mixed layer to zero at the
-surface \autoref{eq:rmtilde}, as described in \autoref{sec:lintaper}. This is
-option (c) of \autoref{fig:eiv_slp}. This linear tapering for the
-slopes used to calculate the eddy-induced fluxes is
-unaffected by the value of \np{ln\_triad\_iso}.
-
-The justification for this linear slope tapering is that, for $A_e$
-that is constant or varies only in the horizontal (the most commonly
-used options in \NEMO: see \autoref{sec:LDF_coef}), it is
-equivalent to a horizontal eiv (eddy-induced velocity) that is uniform
-within the mixed layer \autoref{eq:eiv_v}. This ensures that the
-eiv velocities do not restratify the mixed layer \citep{Treguier1997,
-  Danabasoglu_al_2008}. Equivantly, in terms
-of the skew-flux formulation we use here, the
-linear slope tapering within the mixed-layer gives a linearly varying
-vertical flux, and so a tracer convergence uniform in depth (the
-horizontal flux convergence is relatively insignificant within the mixed-layer).
+The slopes $\tilde{r}_i$ relative to geopotentials (and thus the individual triads \rtriadt{R}) 
+are always tapered linearly from their value immediately below the mixed layer to zero at the surface 
+\autoref{eq:rmtilde}, as described in \autoref{sec:lintaper}. 
+This is option (c) of \autoref{fig:eiv_slp}. 
+This linear tapering for the slopes used to calculate the eddy-induced fluxes is unaffected by 
+the value of \np{ln\_triad\_iso}.
+
+The justification for this linear slope tapering is that, for $A_e$ that is constant or varies only in
+the horizontal (the most commonly used options in \NEMO: see \autoref{sec:LDF_coef}),
+it is equivalent to a horizontal eiv (eddy-induced velocity) that is uniform within the mixed layer
+\autoref{eq:eiv_v}.
+This ensures that the eiv velocities do not restratify the mixed layer \citep{Treguier1997,Danabasoglu_al_2008}.
+Equivantly, in terms of the skew-flux formulation we use here,
+the linear slope tapering within the mixed-layer gives a linearly varying vertical flux,
+and so a tracer convergence uniform in depth
+(the horizontal flux convergence is relatively insignificant within the mixed-layer).
 
 \subsection{Streamfunction diagnostics}\label{sec:sfdiag}
-Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.}, diagnosed
-mean eddy-induced velocities are output. Each time step,
-streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
-$uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$
-(integer $i$, integer +1/2 $j$, integer +1/2 $k$) points (see Table
-\autoref{tab:cell}) respectively. We follow \citep{Griffies_Bk04} and
-calculate the streamfunction at a given $uw$-point from the
-surrounding four triads according to:
+Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.},
+diagnosed mean eddy-induced velocities are output.
+Each time step, streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
+$uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ (integer $i$, integer +1/2 $j$, integer +1/2 $k$)
+points (see Table \autoref{tab:cell}) respectively.
+We follow \citep{Griffies_Bk04} and calculate the streamfunction at a given $uw$-point from
+the surrounding four triads according to:
 \begin{equation}
   \label{eq:sfdiagi}
@@ -1175 +1086 @@
 \end{equation}
 The streamfunction $\psi_1$ is calculated similarly at $vw$ points.
-The eddy-induced velocities are then calculated from the
-straightforward discretisation of \autoref{eq:eiv_v}:
+The eddy-induced velocities are then calculated from the straightforward discretisation of \autoref{eq:eiv_v}:
 \begin{equation}\label{eq:eiv_v_discrete}
 \begin{split}

NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex

@@ -15 +15 @@
 $\ $\newline    % force a new line
 
-The ASM code adds the functionality to apply increments to the model variables: 
-temperature, salinity, sea surface height, velocity and sea ice concentration. 
-These are read into the model from a NetCDF file which may be produced by separate data
-assimilation code.  The code can also output model background fields which are used
-as an input to data assimilation code. This is all controlled by the namelist
-\textit{\ngn{nam\_asminc} }.  There is a brief description of all the namelist options
-provided.  To build the ASM code \key{asminc} must be set.
+The ASM code adds the functionality to apply increments to the model variables: temperature, salinity,
+sea surface height, velocity and sea ice concentration. 
+These are read into the model from a NetCDF file which may be produced by separate data assimilation code.
+The code can also output model background fields which are used as an input to data assimilation code.
+This is all controlled by the namelist \textit{\ngn{nam\_asminc} }.
+There is a brief description of all the namelist options provided.
+To build the ASM code \key{asminc} must be set.
 
 %===============================================================
@@ -28 +28 @@
 \label{sec:ASM_DI}
 
-Direct initialization (DI) refers to the instantaneous correction
-of the model background state using the analysis increment.
+Direct initialization (DI) refers to the instantaneous correction of the model background state using
+the analysis increment.
 DI is used when \np{ln\_asmdin} is set to true.
 
@@ -36 +36 @@
 
 Rather than updating the model state directly with the analysis increment,
-it may be preferable to introduce the increment gradually into the ocean
-model in order to minimize spurious adjustment processes. This technique
-is referred to as Incremental Analysis Updates (IAU) \citep{Bloom_al_MWR96}.
+it may be preferable to introduce the increment gradually into the ocean model in order to
+minimize spurious adjustment processes.
+This technique is referred to as Incremental Analysis Updates (IAU) \citep{Bloom_al_MWR96}.
 IAU is a common technique used with 3D assimilation methods such as 3D-Var or OI.
 IAU is used when \np{ln\_asmiau} is set to true.
 
-With IAU, the model state trajectory ${\bf x}$ in the assimilation window 
-($t_{0} \leq t_{i} \leq t_{N}$)
-is corrected by adding the analysis increments for temperature, salinity, horizontal velocity and SSH
-as additional tendency terms to the prognostic equations:
+With IAU, the model state trajectory ${\bf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$)
+is corrected by adding the analysis increments for temperature, salinity, horizontal velocity and SSH as
+additional tendency terms to the prognostic equations:
 \begin{eqnarray}     \label{eq:wa_traj_iau}
 {\bf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\bf x}^{b}(t_{0})] 
 \; + \; F_{i} \delta \tilde{\bf x}^{a} 
 \end{eqnarray}
-where $F_{i}$ is a weighting function for applying the increments $\delta
-\tilde{\bf x}^{a}$ defined such that $\sum_{i=1}^{N} F_{i}=1$.
-${\bf x}^b$ denotes the model initial state and ${\bf x}^a$ is the model state
-after the increments are applied. 
+where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\bf x}^{a}$ defined such that
+$\sum_{i=1}^{N} F_{i}=1$.
+${\bf x}^b$ denotes the model initial state and ${\bf x}^a$ is the model state after the increments are applied.
 To control the adjustment time of the model to the increment,
-the increment can be applied over an arbitrary sub-window,
-$t_{m} \leq t_{i} \leq t_{n}$, of the main assimilation window,
-where $t_{0} \leq t_{m} \leq t_{i}$ and $t_{i} \leq t_{n} \leq t_{N}$,
+the increment can be applied over an arbitrary sub-window, $t_{m} \leq t_{i} \leq t_{n}$,
+of the main assimilation window, where $t_{0} \leq t_{m} \leq t_{i}$ and $t_{i} \leq t_{n} \leq t_{N}$.
 Typically the increments are spread evenly over the full window.
 In addition, two different weighting functions have been implemented.
@@ -70 +67 @@
 \end{eqnarray}
 where $M = m-n$.
-The second function employs peaked hat-like weights in order to give maximum 
-weight in the centre of the sub-window, with the weighting reduced 
-linearly to a small value at the window end-points:
+The second function employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window,
+with the weighting reduced linearly to a small value at the window end-points:
 \begin{eqnarray}     \label{eq:F2_i}
 F^{(2)}_{i}
@@ -83 +79 @@
 \end{eqnarray}
 where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even. 
-The weights described by \autoref{eq:F2_i} provide a 
-smoother transition of the analysis trajectory from one assimilation cycle 
-to the next than that described by \autoref{eq:F1_i}.
+The weights described by \autoref{eq:F2_i} provide a smoother transition of the analysis trajectory from
+one assimilation cycle to the next than that described by \autoref{eq:F1_i}.
 
 %==========================================================================
@@ -92 +87 @@
 \label{sec:ASM_div_dmp}
 
-The velocity increments may be initialized by the iterative application of 
-a divergence damping operator. In iteration step $n$ new estimates of 
-velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by:
+The velocity increments may be initialized by the iterative application of a divergence damping operator.
+In iteration step $n$ new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by:
 \begin{equation} \label{eq:asm_dmp}
 \left\{ \begin{aligned}
@@ -110 +104 @@
                        +\delta _j \left[ {e_{1v}\,e_{3v}\,v^{n-1}_I} \right]} \right).
 \end{equation}
-By the application of \autoref{eq:asm_dmp} and \autoref{eq:asm_dmp} the divergence is filtered
-in each iteration, and the vorticity is left unchanged. In the presence of coastal boundaries
-with zero velocity increments perpendicular to the coast the divergence is strongly damped.
-This type of the initialisation reduces the vertical velocity magnitude  and alleviates the
-problem of the excessive unphysical vertical mixing in the first steps of the model 
-integration \citep{Talagrand_JAS72, Dobricic_al_OS07}. Diffusion coefficients are defined as 
-$A_D = \alpha e_{1t} e_{2t}$, where $\alpha = 0.2$. The divergence damping is activated by
-assigning to \np{nn\_divdmp} in the \textit{nam\_asminc} namelist a value greater than zero. 
-By choosing this value to be of the order of 100 the increments in the vertical velocity will 
-be significantly reduced.
+By the application of \autoref{eq:asm_dmp} and \autoref{eq:asm_dmp} the divergence is filtered in each iteration,
+and the vorticity is left unchanged.
+In the presence of coastal boundaries with zero velocity increments perpendicular to the coast
+the divergence is strongly damped.
+This type of the initialisation reduces the vertical velocity magnitude and
+alleviates the problem of the excessive unphysical vertical mixing in the first steps of the model integration
+\citep{Talagrand_JAS72, Dobricic_al_OS07}.
+Diffusion coefficients are defined as $A_D = \alpha e_{1t} e_{2t}$, where $\alpha = 0.2$.
+The divergence damping is activated by assigning to \np{nn\_divdmp} in the \textit{nam\_asminc} namelist
+a value greater than zero.
+By choosing this value to be of the order of 100 the increments in
+the vertical velocity will be significantly reduced.
 
 
@@ -127 +123 @@
 \label{sec:ASM_details}
 
-Here we show an example \ngn{namasm} namelist and the header of an example assimilation 
-increments file on the ORCA2 grid.
+Here we show an example \ngn{namasm} namelist and the header of an example assimilation increments file on
+the ORCA2 grid.
 
 %------------------------------------------namasm-----------------------------------------------------

NEMO/trunk/doc/latex/NEMO/subfiles/chap_CONFIG.tex

@@ -18 +18 @@
 
 
-The purpose of this part of the manual is to introduce the \NEMO reference configurations. 
-These configurations are offered as means to explore various numerical and physical options, 
-thus allowing the user to verify that the code is performing in a manner consistent with that 
-we are running. This form of verification is critical as one adopts the code for his or her particular 
-research purposes. The reference configurations also provide a sense for some of the options available 
-in the code, though by no means are all options exercised in the reference configurations.
+The purpose of this part of the manual is to introduce the \NEMO reference configurations.
+These configurations are offered as means to explore various numerical and physical options,
+thus allowing the user to verify that the code is performing in a manner consistent with that we are running.
+This form of verification is critical as one adopts the code for his or her particular research purposes.
+The reference configurations also provide a sense for some of the options available in the code,
+though by no means are all options exercised in the reference configurations.
 
 %------------------------------------------namcfg----------------------------------------------------
@@ -40 +40 @@
 $\ $\newline
 
-The 1D model option simulates a stand alone water column within the 3D \NEMO system. 
-It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers 
-or a biogeochemical model. It is set up by defining the position of the 1D water column in the grid 
+The 1D model option simulates a stand alone water column within the 3D \NEMO system.
+It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers or a biogeochemical model.
+It is set up by defining the position of the 1D water column in the grid
 (see \textit{CONFIG/SHARED/namelist\_ref} ). 
-The 1D model is a very useful tool  
-\textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes ; 
-\textit{(b)} to investigate suitable parameterisations of unresolved turbulence (surface wave
-breaking, Langmuir circulation, ...) ; 
-\textit{(c)} to compare the behaviour of different vertical mixing schemes  ; 
-\textit{(d)} to perform sensitivity studies on the vertical diffusion at a particular point of an ocean domain ; 
+The 1D model is a very useful tool
+\textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes;
+\textit{(b)} to investigate suitable parameterisations of unresolved turbulence
+(surface wave breaking, Langmuir circulation, ...);
+\textit{(c)} to compare the behaviour of different vertical mixing schemes;
+\textit{(d)} to perform sensitivity studies on the vertical diffusion at a particular point of an ocean domain;
 \textit{(d)} to produce extra diagnostics, without the large memory requirement of the full 3D model.
 
-The methodology is based on the use of the zoom functionality over the smallest possible 
-domain : a 3x3 domain centered on the grid point of interest, 
-with some extra routines. There is no need to define a new mesh, bathymetry, 
-initial state or forcing, since the 1D model will use those of the configuration it is a zoom of. 
-The chosen grid point is set in \textit{\ngn{namcfg}} namelist by setting the \np{jpizoom} and \np{jpjzoom} 
-parameters to the indices of the location of the chosen grid point.
-
-The 1D model has some specifies. First, all the horizontal derivatives are assumed to be zero, and
-second, the two components of the velocity are moved on a $T$-point. 
+The methodology is based on the use of the zoom functionality over the smallest possible domain:
+a 3x3 domain centered on the grid point of interest, with some extra routines.
+There is no need to define a new mesh, bathymetry, initial state or forcing,
+since the 1D model will use those of the configuration it is a zoom of.
+The chosen grid point is set in \textit{\ngn{namcfg}} namelist by
+setting the \np{jpizoom} and \np{jpjzoom} parameters to the indices of the location of the chosen grid point.
+
+The 1D model has some specifies. First, all the horizontal derivatives are assumed to be zero,
+and second, the two components of the velocity are moved on a $T$-point. 
 Therefore, defining \key{c1d} changes five main things in the code behaviour: 
 \begin{description}
-\item[(1)] the lateral boundary condition routine (\rou{lbc\_lnk}) set the value of the central column 
-of the 3x3 domain is imposed over the whole domain ; 
-\item[(3)] a call to \rou{lbc\_lnk} is systematically done when reading input data ($i.e.$ in \mdl{iom}) ; 
-\item[(3)] a simplified \rou{stp} routine is used (\rou{stp\_c1d}, see \mdl{step\_c1d} module) in which 
-both lateral tendancy terms and lateral physics are not called ; 
-\item[(4)] the vertical velocity is zero (so far, no attempt at introducing a Ekman pumping velocity 
-has been made) ; 
-\item[(5)] a simplified treatment of the Coriolis term is performed as $U$- and $V$-points are the same 
-(see \mdl{dyncor\_c1d}).
+\item[(1)]
+  the lateral boundary condition routine (\rou{lbc\_lnk}) set the value of the central column of
+  the 3x3 domain is imposed over the whole domain;
+\item[(3)]
+  a call to \rou{lbc\_lnk} is systematically done when reading input data ($i.e.$ in \mdl{iom});
+\item[(3)]
+  a simplified \rou{stp} routine is used (\rou{stp\_c1d}, see \mdl{step\_c1d} module) in which
+  both lateral tendancy terms and lateral physics are not called;
+\item[(4)]
+  the vertical velocity is zero
+  (so far, no attempt at introducing a Ekman pumping velocity has been made);
+\item[(5)]
+  a simplified treatment of the Coriolis term is performed as $U$- and $V$-points are the same
+  (see \mdl{dyncor\_c1d}).
 \end{description}
-All the relevant \textit{\_c1d} modules can be found in the NEMOGCM/NEMO/OPA\_SRC/C1D directory of 
+All the relevant \textit{\_c1d} modules can be found in the NEMOGCM/NEMO/OPA\_SRC/C1D directory of
 the \NEMO distribution.
 
@@ -84 +89 @@
 \label{sec:CFG_orca}
 
-The ORCA family is a series of global ocean configurations that are run together with 
-the LIM sea-ice model (ORCA-LIM) and possibly with PISCES biogeochemical model 
-(ORCA-LIM-PISCES), using various resolutions.
-An appropriate namelist is available in \path{CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg} 
-for ORCA2.
-The domain of ORCA2 configuration is defined in \ifile{ORCA\_R2\_zps\_domcfg} file, this file is available in tar file in the wiki of NEMO : \\
+The ORCA family is a series of global ocean configurations that are run together with
+the LIM sea-ice model (ORCA-LIM) and possibly with PISCES biogeochemical model (ORCA-LIM-PISCES),
+using various resolutions.
+An appropriate namelist is available in \path{CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg} for ORCA2.
+The domain of ORCA2 configuration is defined in \ifile{ORCA\_R2\_zps\_domcfg} file,
+this file is available in tar file in the wiki of NEMO: \\
 https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ReferenceConfigurations/ORCA2\_LIM3\_PISCES \\
 In this namelist\_cfg the name of domain input file is set in \ngn{namcfg} block of namelist. 
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t]   \begin{center}
-\includegraphics[width=0.98\textwidth]{Fig_ORCA_NH_mesh}
-\caption{  \protect\label{fig:MISC_ORCA_msh}     
-ORCA mesh conception. The departure from an isotropic Mercator grid start poleward of 20\degN.
-The two "north pole" are the foci of a series of embedded ellipses (blue curves) 
-which are determined analytically and form the i-lines of the ORCA mesh (pseudo latitudes). 
-Then, following \citet{Madec_Imbard_CD96}, the normal to the series of ellipses (red curves) is computed 
-which provide the j-lines of the mesh (pseudo longitudes).  }
-\end{center}   \end{figure}
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=0.98\textwidth]{Fig_ORCA_NH_mesh}
+    \caption{  \protect\label{fig:MISC_ORCA_msh}
+      ORCA mesh conception.
+      The departure from an isotropic Mercator grid start poleward of 20\degN.
+      The two "north pole" are the foci of a series of embedded ellipses (blue curves) which
+      are determined analytically and form the i-lines of the ORCA mesh (pseudo latitudes).
+      Then, following \citet{Madec_Imbard_CD96}, the normal to the series of ellipses (red curves) is computed which
+      provides the j-lines of the mesh (pseudo longitudes).  }
+  \end{center}
+\end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
@@ -111 +119 @@
 \label{subsec:CFG_orca_grid}
 
-The ORCA grid is a tripolar is based on the semi-analytical method of \citet{Madec_Imbard_CD96}. 
-It allows to construct a global orthogonal curvilinear ocean mesh which has no singularity point inside 
+The ORCA grid is a tripolar is based on the semi-analytical method of \citet{Madec_Imbard_CD96}.
+It allows to construct a global orthogonal curvilinear ocean mesh which has no singularity point inside
 the computational domain since two north mesh poles are introduced and placed on lands.
-The method involves defining an analytical set of mesh parallels in the stereographic polar plan, 
-computing the associated set of mesh meridians, and projecting the resulting mesh onto the sphere. 
-The set of mesh parallels used is a series of embedded ellipses which foci are the two mesh north 
-poles (\autoref{fig:MISC_ORCA_msh}). The resulting mesh presents no loss of continuity in 
-either the mesh lines or the scale factors, or even the scale factor derivatives over the whole 
-ocean domain, as the mesh is not a composite mesh. 
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tbp]  \begin{center}
-\includegraphics[width=1.0\textwidth]{Fig_ORCA_NH_msh05_e1_e2}
-\includegraphics[width=0.80\textwidth]{Fig_ORCA_aniso}
-\caption {  \protect\label{fig:MISC_ORCA_e1e2}
-\textit{Top}: Horizontal scale factors ($e_1$, $e_2$) and 
-\textit{Bottom}: ratio of anisotropy ($e_1 / e_2$)
-for ORCA 0.5\deg ~mesh. South of 20\degN a Mercator grid is used ($e_1 = e_2$) 
-so that the anisotropy ratio is 1. Poleward of 20\degN, the two "north pole" 
-introduce a weak anisotropy over the ocean areas ($< 1.2$) except in vicinity of Victoria Island 
-(Canadian Arctic Archipelago). }
+The method involves defining an analytical set of mesh parallels in the stereographic polar plan,
+computing the associated set of mesh meridians, and projecting the resulting mesh onto the sphere.
+The set of mesh parallels used is a series of embedded ellipses which foci are the two mesh north poles
+(\autoref{fig:MISC_ORCA_msh}).
+The resulting mesh presents no loss of continuity in either the mesh lines or the scale factors,
+or even the scale factor derivatives over the whole ocean domain, as the mesh is not a composite mesh. 
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!tbp]
+  \begin{center}
+    \includegraphics[width=1.0\textwidth]{Fig_ORCA_NH_msh05_e1_e2}
+    \includegraphics[width=0.80\textwidth]{Fig_ORCA_aniso}
+    \caption {  \protect\label{fig:MISC_ORCA_e1e2}
+      \textit{Top}: Horizontal scale factors ($e_1$, $e_2$) and
+      \textit{Bottom}: ratio of anisotropy ($e_1 / e_2$)
+      for ORCA 0.5\deg ~mesh.
+      South of 20\degN a Mercator grid is used ($e_1 = e_2$) so that the anisotropy ratio is 1.
+      Poleward of 20\degN, the two "north pole" introduce a weak anisotropy over the ocean areas ($< 1.2$) except in
+      vicinity of Victoria Island (Canadian Arctic Archipelago). }
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 
-The method is applied to Mercator grid ($i.e.$ same zonal and meridional grid spacing) poleward 
-of 20\degN, so that the Equator is a mesh line, which provides a better numerical solution 
-for equatorial dynamics. The choice of the series of embedded ellipses (position of the foci and 
-variation of the ellipses) is a compromise between maintaining  the ratio of mesh anisotropy 
-($e_1 / e_2$) close to one in the ocean (especially in area of strong eddy activities such as 
-the Gulf Stream) and keeping the smallest scale factor in the northern hemisphere larger 
-than the smallest one in the southern hemisphere.
-The resulting mesh is shown in \autoref{fig:MISC_ORCA_msh} and \autoref{fig:MISC_ORCA_e1e2} 
-for a half a degree grid (ORCA\_R05).
-The smallest ocean scale factor is found in along  Antarctica, while the ratio of anisotropy remains close to one except near the Victoria Island 
-in the Canadian Archipelago. 
+The method is applied to Mercator grid ($i.e.$ same zonal and meridional grid spacing) poleward of 20\degN,
+so that the Equator is a mesh line, which provides a better numerical solution for equatorial dynamics.
+The choice of the series of embedded ellipses (position of the foci and variation of the ellipses)
+is a compromise between maintaining the ratio of mesh anisotropy ($e_1 / e_2$) close to one in the ocean
+(especially in area of strong eddy activities such as the Gulf Stream) and keeping the smallest scale factor in
+the northern hemisphere larger than the smallest one in the southern hemisphere.
+The resulting mesh is shown in \autoref{fig:MISC_ORCA_msh} and \autoref{fig:MISC_ORCA_e1e2} for
+a half a degree grid (ORCA\_R05).
+The smallest ocean scale factor is found in along Antarctica,
+while the ratio of anisotropy remains close to one except near the Victoria Island in the Canadian Archipelago. 
 
 % -------------------------------------------------------------------------------------------------------------
@@ -154 +162 @@
 
 
-The NEMO system is provided with five built-in ORCA configurations which differ in the 
-horizontal resolution. The value of the resolution is given by the resolution at the Equator 
-expressed in degrees. Each of configuration is set through the \textit{domain\_cfg} domain configuration file, 
-which sets the grid size and configuration name parameters. The NEMO System Team provides only ORCA2 domain input file "\ifile{ORCA\_R2\_zps\_domcfg}" file  (Tab. \autoref{tab:ORCA}).
+The NEMO system is provided with five built-in ORCA configurations which differ in the horizontal resolution.
+The value of the resolution is given by the resolution at the Equator expressed in degrees.
+Each of configuration is set through the \textit{domain\_cfg} domain configuration file,
+which sets the grid size and configuration name parameters.
+The NEMO System Team provides only ORCA2 domain input file "\ifile{ORCA\_R2\_zps\_domcfg}" file
+(Tab. \autoref{tab:ORCA}).
 
 
@@ -176 +186 @@
 \hline   \hline
 \end{tabular}
-\caption{ \protect\label{tab:ORCA}   
-Domain size of ORCA family configurations.
-The flag for configurations of ORCA family need to be set in \textit{domain\_cfg} file. }
+\caption{ \protect\label{tab:ORCA}
+  Domain size of ORCA family configurations.
+  The flag for configurations of ORCA family need to be set in \textit{domain\_cfg} file. }
 \end{center}
 \end{table}
@@ -184 +194 @@
 
 
-The ORCA\_R2 configuration has the following specificity : starting from a 2\deg~ORCA mesh, 
-local mesh refinements were applied to the Mediterranean, Red, Black and Caspian Seas, 
-so that the resolution is 1\deg \time 1\deg there. A local transformation were also applied 
-with in the Tropics in order to refine the meridional resolution up to 0.5\deg at the Equator.
-
-The ORCA\_R1 configuration has only a local tropical transformation  to refine the meridional 
-resolution up to 1/3\deg~at the Equator. Note that the tropical mesh refinements in ORCA\_R2 
-and R1 strongly increases the mesh anisotropy there.
+The ORCA\_R2 configuration has the following specificity: starting from a 2\deg~ORCA mesh,
+local mesh refinements were applied to the Mediterranean, Red, Black and Caspian Seas,
+so that the resolution is 1\deg \time 1\deg there.
+A local transformation were also applied with in the Tropics in order to refine the meridional resolution up to
+0.5\deg at the Equator.
+
+The ORCA\_R1 configuration has only a local tropical transformation to refine the meridional resolution up to
+1/3\deg~at the Equator.
+Note that the tropical mesh refinements in ORCA\_R2 and R1 strongly increases the mesh anisotropy there.
 
 The ORCA\_R05 and higher global configurations do not incorporate any regional refinements.  
 
-For ORCA\_R1 and R025, setting the configuration key to 75 allows to use 75 vertical levels, 
-otherwise 46 are used. In the other ORCA configurations, 31 levels are used 
+For ORCA\_R1 and R025, setting the configuration key to 75 allows to use 75 vertical levels, otherwise 46 are used.
+In the other ORCA configurations, 31 levels are used
 (see \autoref{tab:orca_zgr} \sfcomment{HERE I need to put new table for ORCA2 values} and \autoref{fig:zgr}).
 
-Only the ORCA\_R2 is provided with all its input files in the \NEMO distribution. 
-It is very similar to that used as part of the climate model developed at IPSL for the 4th IPCC 
-assessment of climate change (Marti et al., 2009). It is also the basis for the \NEMO contribution 
-to the Coordinate Ocean-ice Reference Experiments (COREs) documented in \citet{Griffies_al_OM09}. 
-
-This version of ORCA\_R2 has 31 levels in the vertical, with the highest resolution (10m) 
-in the upper 150m (see \autoref{tab:orca_zgr} and \autoref{fig:zgr}). 
+Only the ORCA\_R2 is provided with all its input files in the \NEMO distribution.
+It is very similar to that used as part of the climate model developed at IPSL for the 4th IPCC assessment of
+climate change (Marti et al., 2009).
+It is also the basis for the \NEMO contribution to the Coordinate Ocean-ice Reference Experiments (COREs)
+documented in \citet{Griffies_al_OM09}. 
+
+This version of ORCA\_R2 has 31 levels in the vertical, with the highest resolution (10m) in the upper 150m
+(see \autoref{tab:orca_zgr} and \autoref{fig:zgr}). 
 The bottom topography and the coastlines are derived from the global atlas of Smith and Sandwell (1997). 
 The default forcing uses the boundary forcing from \citet{Large_Yeager_Rep04} (see \autoref{subsec:SBC_blk_core}), 
-which was developed for the purpose of running global coupled ocean-ice simulations 
-without an interactive atmosphere. This \citet{Large_Yeager_Rep04} dataset is available 
-through the \href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}. 
-The "normal year" of \citet{Large_Yeager_Rep04} has been chosen of the \NEMO distribution 
-since release v3.3. 
-
-ORCA\_R2 pre-defined configuration can also be run with an AGRIF zoom over the Agulhas 
-current area ( \key{agrif}  defined) and, by setting the appropriate variables, see \path{CONFIG/SHARED/namelist_ref}
-a regional Arctic or peri-Antarctic configuration is extracted from an ORCA\_R2 or R05 configurations
-using sponge layers at open boundaries. 
+which was developed for the purpose of running global coupled ocean-ice simulations without
+an interactive atmosphere.
+This \citet{Large_Yeager_Rep04} dataset is available through
+the \href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}.
+The "normal year" of \citet{Large_Yeager_Rep04} has been chosen of the \NEMO distribution since release v3.3. 
+
+ORCA\_R2 pre-defined configuration can also be run with an AGRIF zoom over the Agulhas current area
+(\key{agrif} defined) and, by setting the appropriate variables, see \path{CONFIG/SHARED/namelist_ref}.
+A regional Arctic or peri-Antarctic configuration is extracted from an ORCA\_R2 or R05 configurations using
+sponge layers at open boundaries. 
 
 % -------------------------------------------------------------------------------------------------------------
@@ -225 +237 @@
 \label{sec:CFG_gyre}
 
-The GYRE configuration \citep{Levy_al_OM10} has been built to simulate
-the seasonal cycle of a double-gyre box model. It consists in an idealized domain 
-similar to that used in the studies of \citet{Drijfhout_JPO94} and \citet{Hazeleger_Drijfhout_JPO98, 
-Hazeleger_Drijfhout_JPO99, Hazeleger_Drijfhout_JGR00, Hazeleger_Drijfhout_JPO00}, 
-over which an analytical seasonal forcing is applied. This allows to investigate the 
-spontaneous generation of a large number of interacting, transient mesoscale eddies 
+The GYRE configuration \citep{Levy_al_OM10} has been built to
+simulate the seasonal cycle of a double-gyre box model.
+It consists in an idealized domain similar to that used in the studies of \citet{Drijfhout_JPO94} and
+\citet{Hazeleger_Drijfhout_JPO98, Hazeleger_Drijfhout_JPO99, Hazeleger_Drijfhout_JGR00, Hazeleger_Drijfhout_JPO00},
+over which an analytical seasonal forcing is applied.
+This allows to investigate the spontaneous generation of a large number of interacting, transient mesoscale eddies 
 and their contribution to the large scale circulation. 
 
-The domain geometry is a closed rectangular basin on the $\beta$-plane centred 
-at $\sim$ 30\degN and rotated by 45\deg, 3180~km long, 2120~km wide 
-and 4~km deep (\autoref{fig:MISC_strait_hand}). 
-The domain is bounded by vertical walls and by a flat bottom. The configuration is 
-meant to represent an idealized North Atlantic or North Pacific basin. 
-The circulation is forced by analytical profiles of wind and buoyancy fluxes. 
-The applied forcings vary seasonally in a sinusoidal manner between winter 
-and summer extrema \citep{Levy_al_OM10}. 
-The wind stress is zonal and its curl changes sign at 22\degN and 36\degN. 
-It forces a subpolar gyre in the north, a subtropical gyre in the wider part of the domain 
-and a small recirculation gyre in the southern corner. 
-The net heat flux takes the form of a restoring toward a zonal apparent air 
-temperature profile. A portion of the net heat flux which comes from the solar radiation
-is allowed to penetrate within the water column. 
-The fresh water flux is also prescribed and varies zonally. 
-It is determined such as, at each time step, the basin-integrated flux is zero. 
-The basin is initialised at rest with vertical profiles of temperature and salinity 
-uniformly applied to the whole domain.
-
-The GYRE configuration is set like an analytical configuration. Through \np{ln\_read\_cfg}\forcode{ = .false.} in \textit{namcfg} namelist defined in the reference configuration \path{CONFIG/GYRE/EXP00/namelist_cfg} anaylitical definition of grid in GYRE is done in usrdef\_hrg, usrdef\_zgr routines. Its horizontal resolution 
-(and thus the size of the domain) is determined by setting \np{nn\_GYRE} in  \ngn{namusr\_def}: \\
+The domain geometry is a closed rectangular basin on the $\beta$-plane centred at $\sim$ 30\degN and
+rotated by 45\deg, 3180~km long, 2120~km wide and 4~km deep (\autoref{fig:MISC_strait_hand}).
+The domain is bounded by vertical walls and by a flat bottom.
+The configuration is meant to represent an idealized North Atlantic or North Pacific basin.
+The circulation is forced by analytical profiles of wind and buoyancy fluxes.
+The applied forcings vary seasonally in a sinusoidal manner between winter and summer extrema \citep{Levy_al_OM10}. 
+The wind stress is zonal and its curl changes sign at 22\degN and 36\degN.
+It forces a subpolar gyre in the north, a subtropical gyre in the wider part of the domain and
+a small recirculation gyre in the southern corner.
+The net heat flux takes the form of a restoring toward a zonal apparent air temperature profile.
+A portion of the net heat flux which comes from the solar radiation is allowed to penetrate within the water column.
+The fresh water flux is also prescribed and varies zonally.
+It is determined such as, at each time step, the basin-integrated flux is zero.
+The basin is initialised at rest with vertical profiles of temperature and salinity uniformly applied to
+the whole domain.
+
+The GYRE configuration is set like an analytical configuration.
+Through \np{ln\_read\_cfg}\forcode{ = .false.} in \textit{namcfg} namelist defined in
+the reference configuration \path{CONFIG/GYRE/EXP00/namelist_cfg}
+analytical definition of grid in GYRE is done in usrdef\_hrg, usrdef\_zgr routines.
+Its horizontal resolution (and thus the size of the domain) is determined by
+setting \np{nn\_GYRE} in \ngn{namusr\_def}: \\
 \np{jpiglo} $= 30 \times$ \np{nn\_GYRE} + 2   \\
 \np{jpjglo} $= 20 \times$ \np{nn\_GYRE} + 2   \\
-Obviously, the namelist parameters have to be adjusted to the chosen resolution, see the Configurations 
-pages on the NEMO web site (Using NEMO\/Configurations) .
+Obviously, the namelist parameters have to be adjusted to the chosen resolution,
+see the Configurations pages on the NEMO web site (Using NEMO\/Configurations).
 In the vertical, GYRE uses the default 30 ocean levels (\jp{jpk}\forcode{ = 31}) (\autoref{fig:zgr}).
 
-The GYRE configuration is also used in benchmark test as it is very simple to increase 
-its resolution and as it does not requires any input file. For example, keeping a same model size 
-on each processor while increasing the number of processor used is very easy, even though the 
-physical integrity of the solution can be compromised. Benchmark is activate via \np{ln\_bench}\forcode{ = .true.} in \ngn{namusr\_def} in namelist \path{CONFIG/GYRE/EXP00/namelist_cfg}.
-
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t]   \begin{center}
-\includegraphics[width=1.0\textwidth]{Fig_GYRE}
-\caption{  \protect\label{fig:GYRE}   
-Snapshot of relative vorticity at the surface of the model domain 
-in GYRE R9, R27 and R54. From \citet{Levy_al_OM10}.}
-\end{center}   \end{figure}
+The GYRE configuration is also used in benchmark test as it is very simple to increase its resolution and
+as it does not requires any input file.
+For example, keeping a same model size on each processor while increasing the number of processor used is very easy,
+even though the physical integrity of the solution can be compromised.
+Benchmark is activate via \np{ln\_bench}\forcode{ = .true.} in \ngn{namusr\_def} in
+namelist \path{CONFIG/GYRE/EXP00/namelist_cfg}.
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=1.0\textwidth]{Fig_GYRE}
+    \caption{  \protect\label{fig:GYRE}
+      Snapshot of relative vorticity at the surface of the model domain in GYRE R9, R27 and R54.
+      From \citet{Levy_al_OM10}.}
+  \end{center}
+\end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
@@ -280 +297 @@
 \label{sec:MISC_config_AMM}
 
-The AMM, Atlantic Margins Model, is a regional model covering the
-Northwest European Shelf domain on a regular lat-lon grid at
-approximately 12km horizontal resolution. The appropriate 
-\textit{\&namcfg} namelist  is available in \textit{CONFIG/AMM12/EXP00/namelist\_cfg}.
+The AMM, Atlantic Margins Model, is a regional model covering the Northwest European Shelf domain on
+a regular lat-lon grid at approximately 12km horizontal resolution.
+The appropriate \textit{\&namcfg} namelist  is available in \textit{CONFIG/AMM12/EXP00/namelist\_cfg}.
 It is used to build the correct dimensions of the AMM domain.
 
-This configuration tests several features of NEMO functionality specific
-to the shelf seas.
-In particular, the AMM uses $S$-coordinates in the vertical rather than
-$z$-coordinates and is forced with tidal lateral boundary conditions
-using a flather boundary condition from the BDY module.
-The AMM configuration  uses the GLS (\key{zdfgls}) turbulence scheme, the
-VVL non-linear free surface(\key{vvl}) and time-splitting
-(\key{dynspg\_ts}).
-
-In addition to the tidal boundary condition the model may also take
-open boundary conditions from a North Atlantic model. Boundaries may be
-completely omitted by setting \np{ln\_bdy} to false.
-Sample surface fluxes, river forcing and a sample initial restart file
-are included to test a realistic model run. The Baltic boundary is
-included within the river input file and is specified as a river source.
-Unlike ordinary river points the Baltic inputs also include salinity and
-temperature data.
+This configuration tests several features of NEMO functionality specific to the shelf seas.
+In particular, the AMM uses $S$-coordinates in the vertical rather than $z$-coordinates and
+is forced with tidal lateral boundary conditions using a flather boundary condition from the BDY module.
+The AMM configuration uses the GLS (\key{zdfgls}) turbulence scheme,
+the VVL non-linear free surface(\key{vvl}) and time-splitting (\key{dynspg\_ts}).
+
+In addition to the tidal boundary condition the model may also take open boundary conditions from
+a North Atlantic model.
+Boundaries may be completely omitted by setting \np{ln\_bdy} to false.
+Sample surface fluxes, river forcing and a sample initial restart file are included to test a realistic model run.
+The Baltic boundary is included within the river input file and is specified as a river source.
+Unlike ordinary river points the Baltic inputs also include salinity and temperature data.
 
 \end{document}

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex

@@ -17 +17 @@
 \label{sec:DIA_io_old}
 
-The model outputs are of three types: the restart file, the output listing, and 
-the diagnostic output file(s).
-The restart file is used internally by the code when the user wants to start the model with 
+The model outputs are of three types: the restart file, the output listing, and the diagnostic output file(s).
+The restart file is used internally by the code when the user wants to start the model with
 initial conditions defined by a previous simulation.
-It contains all the information that is necessary in order for there to be no changes in 
-the model results (even at the computer precision) between a run performed with several restarts and 
+It contains all the information that is necessary in order for there to be no changes in the model results
+(even at the computer precision) between a run performed with several restarts and
 the same run performed in one step.
-It should be noted that this requires that the restart file contains two consecutive time steps for 
-all the prognostic variables, and that it is saved in the same binary format as the one used by 
-the computer that is to read it (in particular, 32 bits binary IEEE format must not be used for this file).
-
-The output listing and file(s) are predefined but should be checked and eventually adapted to 
-the user's needs.
+It should be noted that this requires that the restart file contains two consecutive time steps for
+all the prognostic variables, and that it is saved in the same binary format as the one used by the computer that
+is to read it (in particular, 32 bits binary IEEE format must not be used for this file).
+
+The output listing and file(s) are predefined but should be checked and eventually adapted to the user's needs.
 The output listing is stored in the $ocean.output$ file.
 The information is printed from within the code on the logical unit $numout$.
@@ -35 +33 @@
 
 By default, diagnostic output files are written in NetCDF format.
-Since version 3.2, when defining \key{iomput}, an I/O server has been added which 
-provides more flexibility in the choice of the fields to be written as well as 
-how the writing work is distributed over the processors in massively parallel computing. 
-A complete description of the use of this I/O server is presented in the next section. 
-
-By default, \key{iomput} is not defined, NEMO produces NetCDF with the old IOIPSL library which 
-has been kept for compatibility and its easy installation.
-However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2, 
-many diagnostic options have been added presuming the use of \key{iomput}. 
-The usefulness of the default IOIPSL-based option is expected to reduce with each new release. 
-If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module and 
-contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period of 
+Since version 3.2, when defining \key{iomput}, an I/O server has been added which
+provides more flexibility in the choice of the fields to be written as well as how
+the writing work is distributed over the processors in massively parallel computing.
+A complete description of the use of this I/O server is presented in the next section.
+
+By default, \key{iomput} is not defined,
+NEMO produces NetCDF with the old IOIPSL library which has been kept for compatibility and its easy installation.
+However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2,
+many diagnostic options have been added presuming the use of \key{iomput}.
+The usefulness of the default IOIPSL-based option is expected to reduce with each new release.
+If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module and
+contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period of
 nn\_write time-steps (namelist parameter). 
 
@@ -57 +55 @@
 \label{sec:DIA_iom}
 
-Since version 3.2, iomput is the NEMO output interface of choice. 
-It has been designed to be simple to use, flexible and efficient. 
+Since version 3.2, iomput is the NEMO output interface of choice.
+It has been designed to be simple to use, flexible and efficient.
 The two main purposes of iomput are: 
 
 \begin{enumerate}
-   \item The complete and flexible control of the output files through external XML files adapted by 
-         the user from standard templates.
-   \item To achieve high performance and scalable output through the optional distribution of 
-         all diagnostic output related tasks to dedicated processes.
+\item
+  The complete and flexible control of the output files through external XML files adapted by
+  the user from standard templates.
+\item
+  To achieve high performance and scalable output through the optional distribution of
+  all diagnostic output related tasks to dedicated processes.
 \end{enumerate}
 
@@ -72 +72 @@
 
 \begin{itemize}
-   \item The choice of output frequencies that can be different for each file 
-         (including real months and years).
-   \item The choice of file contents; includes complete flexibility over which data are written in 
-         which files (the same data can be written in different files).
-   \item The possibility to split output files at a chosen frequency.
-   \item The possibility to extract a vertical or an horizontal subdomain.
-   \item The choice of the temporal operation to perform, 
-         e.g.: average, accumulate, instantaneous, min, max and once.
-   \item Control over metadata via a large XML "database" of possible output fields.
+\item
+  The choice of output frequencies that can be different for each file (including real months and years).
+\item
+  The choice of file contents; includes complete flexibility over which data are written in which files
+  (the same data can be written in different files).
+\item
+  The possibility to split output files at a chosen frequency.
+\item
+  The possibility to extract a vertical or an horizontal subdomain.
+\item
+  The choice of the temporal operation to perform, $e.g.$: average, accumulate, instantaneous, min, max and once.
+\item
+  Control over metadata via a large XML "database" of possible output fields.
 \end{itemize}
 
-In addition, iomput allows the user to add in the code the output of any new variable (scalar, 2D or 
-3D) in a very easy way.
+In addition, iomput allows the user to add in the code the output of any new variable (scalar, 2D or 3D)
+in a very easy way.
 All details of iomput functionalities are listed in the following subsections.
-Examples of the XML files that control the outputs can be found in: 
-\path{NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml}, \path{NEMOGCM/CONFIG/SHARED/field_def.xml} and 
-\path{NEMOGCM/CONFIG/SHARED/domain_def.xml}. \\
+Examples of the XML files that control the outputs can be found in: \path{NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml},
+\path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml}. \\
 
 The second functionality targets output performance when running in parallel (\key{mpp\_mpi}).
-Iomput provides the possibility to specify N dedicated I/O processes (in addition to 
-the NEMO processes) to collect and write the outputs.
-With an appropriate choice of N by the user, the bottleneck associated with the writing of 
+Iomput provides the possibility to specify N dedicated I/O processes (in addition to the NEMO processes)
+to collect and write the outputs.
+With an appropriate choice of N by the user, the bottleneck associated with the writing of
 the output files can be greatly reduced.
 
-In version 3.6, the iom\_put interface depends on an external code called 
-\href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-1.0}{XIOS-1.0} 
+In version 3.6, the iom\_put interface depends on
+an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-1.0}{XIOS-1.0} 
 (use of revision 618 or higher is required).
-This new IO server can take advantage of the parallel I/O functionality of NetCDF4 to 
+This new IO server can take advantage of the parallel I/O functionality of NetCDF4 to
 create a single output file and therefore to bypass the rebuilding phase.
-Note that writing in parallel into the same NetCDF files requires that 
-your NetCDF4 library is linked to an HDF5 library that has been correctly compiled ($i.e.$ with 
-the configure option $--$enable-parallel).
+Note that writing in parallel into the same NetCDF files requires that your NetCDF4 library is linked to
+an HDF5 library that has been correctly compiled ($i.e.$ with the configure option $--$enable-parallel).
 Note that the files created by iomput through XIOS are incompatible with NetCDF3.
-All post-processsing and visualization tools must therefore be compatible with 
-NetCDF4 and not only NetCDF3.
-
-Even if not using the parallel I/O functionality of NetCDF4, using N dedicated I/O servers, 
-where N is typically much less than the number of NEMO processors, 
-will reduce the number of output files created.
-This can greatly reduce the post-processing burden usually associated with using 
-large numbers of NEMO processors.
-Note that for smaller configurations, the rebuilding phase can be avoided, even without 
-a parallel-enabled NetCDF4 library, simply by employing only one dedicated I/O server.
+All post-processsing and visualization tools must therefore be compatible with NetCDF4 and not only NetCDF3.
+
+Even if not using the parallel I/O functionality of NetCDF4, using N dedicated I/O servers,
+where N is typically much less than the number of NEMO processors, will reduce the number of output files created.
+This can greatly reduce the post-processing burden usually associated with using large numbers of NEMO processors.
+Note that for smaller configurations, the rebuilding phase can be avoided,
+even without a parallel-enabled NetCDF4 library, simply by employing only one dedicated I/O server.
 
 \subsection{XIOS: XML Inputs-Outputs Server}
@@ -120 +119 @@
 \subsubsection{Attached or detached mode?}
 
-Iomput is based on \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS}, 
+Iomput is based on \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS},
 the io\_server developed by Yann Meurdesoif from IPSL.
 The behaviour of the I/O subsystem is controlled by settings in the external XML files listed above.
@@ -127 +126 @@
 \xmlline|<variable id="using_server" type="bool"></variable>|
 
-The {\tt using\_server} setting determines whether or not the server will be used in 
-\textit{attached mode} (as a library) [{\tt> false <}] or in \textit{detached mode} (as 
-an external executable on N additional, dedicated cpus) [{\tt > true <}].
+The {\tt using\_server} setting determines whether or not the server will be used in \textit{attached mode}
+(as a library) [{\tt> false <}] or in \textit{detached mode}
+(as an external executable on N additional, dedicated cpus) [{\tt > true <}].
 The \textit{attached mode} is simpler to use but much less efficient for massively parallel applications.
 The type of each file can be either ''multiple\_file'' or ''one\_file''.
 
-In \textit{attached mode} and if the type of file is ''multiple\_file'', 
+In \textit{attached mode} and if the type of file is ''multiple\_file'',
 then each NEMO process will also act as an IO server and produce its own set of output files.
 Superficially, this emulates the standard behaviour in previous versions.
-However, the subdomain written out by each process does not correspond to 
+However, the subdomain written out by each process does not correspond to
 the \forcode{jpi x jpj x jpk} domain actually computed by the process (although it may if \forcode{jpni=1}).
 Instead each process will have collected and written out a number of complete longitudinal strips.
-If the ''one\_file'' option is chosen then all processes will collect their longitudinal strips and 
+If the ''one\_file'' option is chosen then all processes will collect their longitudinal strips and
 write (in parallel) to a single output file.
 
-In \textit{detached mode} and if the type of file is ''multiple\_file'', then 
-each stand-alone XIOS process will collect data for a range of complete longitudinal strips and 
+In \textit{detached mode} and if the type of file is ''multiple\_file'',
+then each stand-alone XIOS process will collect data for a range of complete longitudinal strips and
 write to its own set of output files.
-If the ''one\_file'' option is chosen then all XIOS processes will collect their longitudinal strips and 
+If the ''one\_file'' option is chosen then all XIOS processes will collect their longitudinal strips and
 write (in parallel) to a single output file. 
 Note running in detached mode requires launching a Multiple Process Multiple Data (MPMD) parallel job.
-The following subsection provides a typical example but the syntax will vary in 
-different MPP environments.
+The following subsection provides a typical example but the syntax will vary in different MPP environments.
 
 \subsubsection{Number of cpu used by XIOS in detached mode}
@@ -156 +154 @@
 The number of cores dedicated to XIOS should be from \texttildelow1/10 to \texttildelow1/50 of the number of 
 cores dedicated to NEMO.
-Some manufacturers suggest using O($\sqrt{N}$) dedicated IO processors for N processors but 
+Some manufacturers suggest using O($\sqrt{N}$) dedicated IO processors for N processors but
 this is a general recommendation and not specific to NEMO.
-It is difficult to provide precise recommendations because the optimal choice will depend on 
+It is difficult to provide precise recommendations because the optimal choice will depend on
 the particular hardware properties of the target system 
-(parallel filesystem performance, available memory, memory bandwidth etc.) and 
-the volume and frequency of data to be created.
+(parallel filesystem performance, available memory, memory bandwidth etc.)
+and the volume and frequency of data to be created.
 Here is an example of 2 cpus for the io\_server and 62 cpu for nemo using mpirun:
 \cmd|mpirun -np 62 ./nemo.exe : -np 2 ./xios_server.exe|
@@ -167 +165 @@
 \subsubsection{Control of XIOS: the context in iodef.xml}
 
-As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in 
-the XIOS context in iodef.xml.
+As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml.
 See the XML basics section below for more details on XML syntax and rules.
 
@@ -205 +202 @@
 \subsubsection{Installation}
 
-As mentioned, XIOS is supported separately and must be downloaded and compiled before 
-it can be used with NEMO.
-See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for 
-help and guidance.
+As mentioned, XIOS is supported separately and must be downloaded and compiled before it can be used with NEMO.
+See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for help and guidance.
 NEMO will need to link to the compiled XIOS library.
-The 
-\href{https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ModelInterfacing/InputsOutputs#Inputs-OutputsusingXIOS}
-{XIOS with NEMO}
-guide provides an example illustration of how this can be achieved.
+The \href{https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ModelInterfacing/InputsOutputs#Inputs-OutputsusingXIOS}
+{XIOS with NEMO} guide provides an example illustration of how this can be achieved.
 
 \subsubsection{Add your own outputs}
 
 It is very easy to add your own outputs with iomput.
-Many standard fields and diagnostics are already prepared ($i.e.$, steps 1 to 3 below have been done) and 
+Many standard fields and diagnostics are already prepared ($i.e.$, steps 1 to 3 below have been done) and
 simply need to be activated by including the required output in a file definition in iodef.xml (step 4).
 To add new output variables, all 4 of the following steps must be taken.
 
 \begin{enumerate}
-   \item[1.] in NEMO code, add a \forcode{CALL iom\_put( 'identifier', array )} where you want to 
-               output a 2D or 3D array.
-   \item[2.] If necessary, add \forcode{USE iom ! I/O manager library} to the list of used modules in 
-               the upper part of your module.
-   \item[3.] in the field\_def.xml file, add the definition of your variable using 
-               the same identifier you used in the f90 code (see subsequent sections for 
-               a details of the XML syntax and rules).
-For example:
+\item[1.]
+  in NEMO code, add a \forcode{CALL iom\_put( 'identifier', array )} where you want to output a 2D or 3D array.
+\item[2.]
+  If necessary, add \forcode{USE iom ! I/O manager library} to the list of used modules in
+  the upper part of your module.
+\item[3.]
+  in the field\_def.xml file, add the definition of your variable using the same identifier you used in the f90 code
+  (see subsequent sections for a details of the XML syntax and rules).
+  For example:
 
 \begin{xmllines}
@@ -241 +235 @@
 \end{xmllines}
 
-Note your definition must be added to the field\_group whose reference grid is consistent with 
-the size of the array passed to iomput.
-The grid\_ref attribute refers to definitions set in iodef.xml which, in turn, reference grids and 
-axes either defined in the code (iom\_set\_domain\_attr and iom\_set\_axis\_attr in \mdl{iom}) or 
-defined in the domain\_def.xml file.
+Note your definition must be added to the field\_group whose reference grid is consistent with the size of
+the array passed to iomput.
+The grid\_ref attribute refers to definitions set in iodef.xml which, in turn,
+reference grids and axes either defined in the code
+(iom\_set\_domain\_attr and iom\_set\_axis\_attr in \mdl{iom}) or defined in the domain\_def.xml file.
 $e.g.$:
 
@@ -252 +246 @@
 \end{xmllines}
 
-Note, if your array is computed within the surface module each \np{nn\_fsbc} time\_step, 
+Note, if your array is computed within the surface module each \np{nn\_fsbc} time\_step,
 add the field definition within the field\_group defined with the id "SBC":
-\xmlcode{<field_group id="SBC" ...>} which has been defined with the correct frequency of 
-operations (iom\_set\_field\_attr in \mdl{iom})
-   \item[4.] add your field in one of the output files defined in iodef.xml 
-(again see subsequent sections for syntax and rules)
+\xmlcode{<field_group id="SBC" ...>} which has been defined with the correct frequency of operations
+(iom\_set\_field\_attr in \mdl{iom})
+\item[4.]
+  add your field in one of the output files defined in iodef.xml
+  (again see subsequent sections for syntax and rules)
 
 \begin{xmllines}
@@ -274 +269 @@
 
 XML tags begin with the less-than character ("$<$") and end with the greater-than character ("$>$").
-You use tags to mark the start and end of elements, which are the logical units of information in 
-an XML document.
-In addition to marking the beginning of an element, XML start tags also provide a place to 
-specify attributes.
-An attribute specifies a single property for an element, using a name/value pair, for example: 
+You use tags to mark the start and end of elements, which are the logical units of information in an XML document.
+In addition to marking the beginning of an element, XML start tags also provide a place to specify attributes.
+An attribute specifies a single property for an element, using a name/value pair, for example:
 \xmlcode{<a b="x" c="y" d="z"> ... </a>}.
 See \href{http://www.xmlnews.org/docs/xml-basics.html}{here} for more details.
@@ -284 +277 @@
 \subsubsection{Structure of the XML file used in NEMO}
 
-The XML file used in XIOS is structured by 7 families of tags: 
+The XML file used in XIOS is structured by 7 families of tags:
 context, axis, domain, grid, field, file and variable.
 Each tag family has hierarchy of three flavors (except for context):
@@ -302 +295 @@
 
 Each element may have several attributes.
-Some attributes are mandatory, other are optional but have a default value and 
-other are completely optional.
+Some attributes are mandatory, other are optional but have a default value and other are completely optional.
 Id is a special attribute used to identify an element or a group of elements.
 It must be unique for a kind of element.
 It is optional, but no reference to the corresponding element can be done if it is not defined.
 
-The XML file is split into context tags that are used to isolate IO definition from different codes or 
-different parts of a code.
+The XML file is split into context tags that are used to isolate IO definition from
+different codes or different parts of a code.
 No interference is possible between 2 different contexts.
 Each context has its own calendar and an associated timestep.
@@ -370 +362 @@
  
 \noindent In NEMO, by default, the field and domain definition is done in 2 separate files:
-\path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml} 
-that are included in the main iodef.xml file through the following commands:
+\path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml} that
+are included in the main iodef.xml file through the following commands:
 \begin{xmllines}
 <field_definition src="./field_def.xml" />
@@ -380 +372 @@
 
 XML extensively uses the concept of inheritance.
-XML has a tree based structure with a parent-child oriented relation: 
-all children inherit attributes from parent, but an attribute defined in a child replace 
-the inherited attribute value. 
-Note that the special attribute ''id'' is never inherited. \\ \\
+XML has a tree based structure with a parent-child oriented relation: all children inherit attributes from parent,
+but an attribute defined in a child replace the inherited attribute value.
+Note that the special attribute ''id'' is never inherited.
+\\
+\\
 example 1: Direct inheritance.
 
@@ -393 +386 @@
 \end{xmllines}
 
-The field ''sst'' which is part (or a child) of the field\_definition will inherit the value ''average'' 
-of the attribute ''operation'' from its parent.
+The field ''sst'' which is part (or a child) of the field\_definition will inherit the value ''average'' of
+the attribute ''operation'' from its parent.
 Note that a child can overwrite the attribute definition inherited from its parents.
-In the example above, the field ''sss'' will for example output instantaneous values instead of 
-average values. \\ \\
+In the example above, the field ''sss'' will for example output instantaneous values instead of average values.
+\\
+\\
 example 2: Inheritance by reference.
 
@@ -418 +412 @@
 
 Groups can be used for 2 purposes.
-Firstly, the group can be used to define common attributes to be shared by the elements of 
+Firstly, the group can be used to define common attributes to be shared by the elements of
 the group through inheritance.
 In the following example, we define a group of field that will share a common grid ''grid\_T\_2D''.
-Note that for the field ''toce'', we overwrite the grid definition inherited from the group by 
-''grid\_T\_3D''.
+Note that for the field ''toce'', we overwrite the grid definition inherited from the group by ''grid\_T\_3D''.
 
 \begin{xmllines}
@@ -456 +449 @@
 \subsection{Detailed functionalities}
 
-The file \path{NEMOGCM/CONFIG/ORCA2_LIM/iodef_demo.xml} provides several examples of the use of 
+The file \path{NEMOGCM/CONFIG/ORCA2_LIM/iodef_demo.xml} provides several examples of the use of
 the new functionalities offered by the XML interface of XIOS.
 
 \subsubsection{Define horizontal subdomains}
 
-Horizontal subdomains are defined through the attributs zoom\_ibegin, zoom\_jbegin, zoom\_ni, zoom\_nj of 
+Horizontal subdomains are defined through the attributs zoom\_ibegin, zoom\_jbegin, zoom\_ni, zoom\_nj of
 the tag family domain.
 It must therefore be done in the domain part of the XML file. 
@@ -472 +465 @@
 \end{xmllines}
 
-The use of this subdomain is done through the redefinition of the attribute domain\_ref of 
-the tag family field.
+The use of this subdomain is done through the redefinition of the attribute domain\_ref of the tag family field.
 For example:
 
@@ -483 +475 @@
 
 Moorings are seen as an extrem case corresponding to a 1 by 1 subdomain. 
-The Equatorial section, the TAO, RAMA and PIRATA moorings are alredy registered in the code and 
+The Equatorial section, the TAO, RAMA and PIRATA moorings are already registered in the code and
 can therefore be outputted without taking care of their (i,j) position in the grid.
 These predefined domains can be activated by the use of specific domain\_ref:
-''EqT'', ''EqU'' or ''EqW'' for the equatorial sections and the mooring position for 
-TAO, RAMA and PIRATA followed by ''T'' (for example: ''8s137eT'', ''1.5s80.5eT'' ...)
+''EqT'', ''EqU'' or ''EqW'' for the equatorial sections and
+the mooring position for TAO, RAMA and PIRATA followed by ''T'' (for example: ''8s137eT'', ''1.5s80.5eT'' ...)
 
 \begin{xmllines}
@@ -503 +495 @@
 \subsubsection{Define vertical zooms}
 
-Vertical zooms are defined through the attributs zoom\_begin and zoom\_end of the tag family axis. 
-It must therefore be done in the axis part of the XML file. 
+Vertical zooms are defined through the attributs zoom\_begin and zoom\_end of the tag family axis.
+It must therefore be done in the axis part of the XML file.
 For example, in \path{NEMOGCM/CONFIG/ORCA2_LIM/iodef_demo.xml}, we provide the following example:
 
@@ -513 +505 @@
 \end{xmllines}
 
-The use of this vertical zoom is done through the redefinition of the attribute axis\_ref of 
-the tag family field.
+The use of this vertical zoom is done through the redefinition of the attribute axis\_ref of the tag family field.
 For example:
 
@@ -539 +530 @@
 \end{xmllines}
 
-However it is often very convienent to define the file name with the name of the experiment, 
-the output file frequency and the date of the beginning and the end of the simulation 
+However it is often very convienent to define the file name with the name of the experiment,
+the output file frequency and the date of the beginning and the end of the simulation
 (which are informations stored either in the namelist or in the XML file).
-To do so, we added the following rule: if the id of the tag file is ''fileN'' (where N = 1 to 999 on 
-1 to 3 digits) or one of the predefined sections or moorings (see next subsection), 
-the following part of the name and the name\_suffix (that can be inherited) will be automatically 
-replaced by:
+To do so, we added the following rule:
+if the id of the tag file is ''fileN'' (where N = 1 to 999 on 1 to 3 digits) or
+one of the predefined sections or moorings (see next subsection),
+the following part of the name and the name\_suffix (that can be inherited) will be automatically replaced by:
 
 \begin{table} \scriptsize
@@ -580 +571 @@
 \end{forlines}
 
-\noindent will give the following file name radical:
-\ifile{myfile\_ORCA2\_19891231\_freq1d}
+\noindent will give the following file name radical: \ifile{myfile\_ORCA2\_19891231\_freq1d}
 
 \subsubsection{Other controls of the XML attributes from NEMO}
 
 The values of some attributes are defined by subroutine calls within NEMO 
-(calls to iom\_set\_domain\_attr, iom\_set\_axis\_attr and iom\_set\_field\_attr in \mdl{iom}). 
-Any definition given in the XML file will be overwritten. 
-By convention, these attributes are defined to ''auto'' (for string) or ''0000'' (for integer) in 
-the XML file (but this is not necessary). \\
-
-Here is the list of these attributes: \\
+(calls to iom\_set\_domain\_attr, iom\_set\_axis\_attr and iom\_set\_field\_attr in \mdl{iom}).
+Any definition given in the XML file will be overwritten.
+By convention, these attributes are defined to ''auto'' (for string) or ''0000'' (for integer) in the XML file
+(but this is not necessary).
+\\
+
+Here is the list of these attributes:
+\\
 
 \begin{table} \scriptsize
@@ -631 +623 @@
 
 \begin{enumerate}
-\item Simple computation: directly define the computation when refering to the variable in 
-the file definition.
+\item
+  Simple computation: directly define the computation when refering to the variable in the file definition.
 
 \begin{xmllines}
@@ -640 +632 @@
 \end{xmllines}
 
-\item Simple computation: define a new variable and use it in the file definition.
+\item
+  Simple computation: define a new variable and use it in the file definition.
 
 in field\_definition:
@@ -654 +647 @@
 \end{xmllines}
 
-Note that in this case, the following syntaxe \xmlcode{<field field_ref="sst2" />} is not working as 
+Note that in this case, the following syntaxe \xmlcode{<field field_ref="sst2" />} is not working as
 sst2 won't be evaluated.
 
-\item Change of variable precision:
+\item
+  Change of variable precision:
 
 \begin{xmllines}
@@ -667 +661 @@
 
 Note that, then the code is crashing, writting real4 variables forces a numerical convection from 
-real8 to real4 which will create an internal error in NetCDF and will avoid the creation of 
-the output files.
-Forcing double precision outputs with prec="8" (for example in the field\_definition) will 
-avoid this problem.
-
-\item add user defined attributes:
+real8 to real4 which will create an internal error in NetCDF and will avoid the creation of the output files.
+Forcing double precision outputs with prec="8" (for example in the field\_definition) will avoid this problem.
+
+\item
+  add user defined attributes:
 
 \begin{xmllines}
@@ -687 +680 @@
 \end{xmllines}
 
-\item use of the ``@'' function: example 1, weighted temporal average
+\item
+  use of the ``@'' function: example 1, weighted temporal average
 
  - define a new variable in field\_definition
@@ -706 +700 @@
 
 The freq\_op="5d" attribute is used to define the operation frequency of the ``@'' function: here 5 day.
-The temporal operation done by the ``@'' is the one defined in the field definition: 
+The temporal operation done by the ``@'' is the one defined in the field definition:
 here we use the default, average.
 So, in the above case, @toce\_e3t will do the 5-day mean of toce*e3t.
-Operation="instant" refers to the temporal operation to be performed on the field''@toce\_e3t / @e3t'': 
-here the temporal average is alreday done by the ``@'' function so we just use instant to do the ratio of 
+Operation="instant" refers to the temporal operation to be performed on the field''@toce\_e3t / @e3t'':
+here the temporal average is alreday done by the ``@'' function so we just use instant to do the ratio of
 the 2 mean values.
 field\_ref="toce" means that attributes not explicitely defined, are inherited from toce field.
 Note that in this case, freq\_op must be equal to the file output\_freq.
 
-\item use of the ``@'' function: example 2, monthly SSH standard deviation
+\item
+  use of the ``@'' function: example 2, monthly SSH standard deviation
 
  - define a new variable in field\_definition
@@ -737 +732 @@
 
 The freq\_op="1m" attribute is used to define the operation frequency of the ``@'' function: here 1 month.
-The temporal operation done by the ``@'' is the one defined in the field definition: 
+The temporal operation done by the ``@'' is the one defined in the field definition:
 here we use the default, average.
 So, in the above case, @ssh2 will do the monthly mean of ssh*ssh.
-Operation="instant" refers to the temporal operation to be performed on 
-the field ''sqrt( @ssh2 - @ssh * @ssh )'': 
+Operation="instant" refers to the temporal operation to be performed on the field ''sqrt( @ssh2 - @ssh * @ssh )'': 
 here the temporal average is alreday done by the ``@'' function so we just use instant.
 field\_ref="ssh" means that attributes not explicitely defined, are inherited from ssh field.
 Note that in this case, freq\_op must be equal to the file output\_freq.
 
-\item use of the ``@'' function: example 3, monthly average of SST diurnal cycle
+\item
+  use of the ``@'' function: example 3, monthly average of SST diurnal cycle
 
  - define 2 new variables in field\_definition
@@ -770 +765 @@
 
 The freq\_op="1d" attribute is used to define the operation frequency of the ``@'' function: here 1 day.
-The temporal operation done by the ``@'' is the one defined in the field definition: here maximum for 
-sstmax and minimum for sstmin.
+The temporal operation done by the ``@'' is the one defined in the field definition:
+here maximum for sstmax and minimum for sstmin.
 So, in the above case, @sstmax will do the daily max and @sstmin the daily min.
-Operation="average" refers to the temporal operation to be performed on the field ``@sstmax - @sstmin'': 
+Operation="average" refers to the temporal operation to be performed on the field ``@sstmax - @sstmin'':
 here monthly mean (of daily max - daily min of the sst).
 field\_ref="sst" means that attributes not explicitely defined, are inherited from sst field.
@@ -1126 +1121 @@
 
 Output from the XIOS-1.0 IO server is compliant with 
-\href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} 
-of the CF metadata standard. 
-Therefore while a user may wish to add their own metadata to the output files (as demonstrated in 
-example 4 of section \autoref{subsec:IOM_xmlref}) the metadata should, for the most part, comply with 
-the CF-1.5 standard.
-
-Some metadata that may significantly increase the file size (horizontal cell areas and 
-vertices) are controlled by the namelist parameter \np{ln\_cfmeta} in the \ngn{namrun} namelist.
+\href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} of
+the CF metadata standard. 
+Therefore while a user may wish to add their own metadata to the output files (as demonstrated in example 4 of
+section \autoref{subsec:IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard.
+
+Some metadata that may significantly increase the file size (horizontal cell areas and vertices) are controlled by
+the namelist parameter \np{ln\_cfmeta} in the \ngn{namrun} namelist.
 This must be set to true if these metadata are to be included in the output files.
 
@@ -1144 +1138 @@
 
 Since version 3.3, support for NetCDF4 chunking and (loss-less) compression has been included.
-These options build on the standard NetCDF output and allow the user control over the size of 
-the chunks via namelist settings.
+These options build on the standard NetCDF output and allow the user control over the size of the chunks via
+namelist settings.
 Chunking and compression can lead to significant reductions in file sizes for a small runtime overhead.
-For a fuller discussion on chunking and other performance issues the reader is referred to 
-the NetCDF4 documentation found 
-\href{http://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#Chunking}{here}.
-
-The new features are only available when the code has been linked with a NetCDF4 library 
-(version 4.1 onwards, recommended) which has been built with HDF5 support 
-(version 1.8.4 onwards, recommended).
-Datasets created with chunking and compression are not backwards compatible with 
-NetCDF3 "classic" format but most analysis codes can be relinked simply with the new libraries and 
-will then read both NetCDF3 and NetCDF4 files.
-NEMO executables linked with NetCDF4 libraries can be made to produce NetCDF3 files by 
+For a fuller discussion on chunking and other performance issues the reader is referred to
+the NetCDF4 documentation found \href{http://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#Chunking}{here}.
+
+The new features are only available when the code has been linked with a NetCDF4 library
+(version 4.1 onwards, recommended) which has been built with HDF5 support (version 1.8.4 onwards, recommended).
+Datasets created with chunking and compression are not backwards compatible with NetCDF3 "classic" format but
+most analysis codes can be relinked simply with the new libraries and will then read both NetCDF3 and NetCDF4 files.
+NEMO executables linked with NetCDF4 libraries can be made to produce NetCDF3 files by
 setting the \np{ln\_nc4zip} logical to false in the \textit{namnc4} namelist:
 
@@ -1166 +1157 @@
 
 If \key{netcdf4} has not been defined, these namelist parameters are not read.
-In this case, \np{ln\_nc4zip} is set false and dummy routines for 
-a few NetCDF4-specific functions are defined.
-These functions will not be used but need to be included so that compilation is possible with 
-NetCDF3 libraries.
-
-When using NetCDF4 libraries, \key{netcdf4} should be defined even if the intention is to 
+In this case, \np{ln\_nc4zip} is set false and dummy routines for a few NetCDF4-specific functions are defined.
+These functions will not be used but need to be included so that compilation is possible with NetCDF3 libraries.
+
+When using NetCDF4 libraries, \key{netcdf4} should be defined even if the intention is to
 create only NetCDF3-compatible files.
-This is necessary to avoid duplication between the dummy routines and 
-the actual routines present in the library.
+This is necessary to avoid duplication between the dummy routines and the actual routines present in the library.
 Most compilers will fail at compile time when faced with such duplication.
-Thus when linking with NetCDF4 libraries the user must define \key{netcdf4} and 
+Thus when linking with NetCDF4 libraries the user must define \key{netcdf4} and
 control the type of NetCDF file produced via the namelist parameter.
 
-Chunking and compression is applied only to 4D fields and there is no advantage in 
-chunking across more than one time dimension since previously written chunks would have to 
-be read back and decompressed before being added to.
+Chunking and compression is applied only to 4D fields and
+there is no advantage in chunking across more than one time dimension since
+previously written chunks would have to be read back and decompressed before being added to.
 Therefore, user control over chunk sizes is provided only for the three space dimensions.
 The user sets an approximate number of chunks along each spatial axis.
-The actual size of the chunks will depend on global domain size for mono-processors or, 
-more likely, the local processor domain size for distributed processing.
-The derived values are subject to practical minimum values (to avoid wastefully small chunk sizes) and 
+The actual size of the chunks will depend on global domain size for mono-processors or, more likely,
+the local processor domain size for distributed processing.
+The derived values are subject to practical minimum values (to avoid wastefully small chunk sizes) and
 cannot be greater than the domain size in any dimension.
 The algorithm used is:
@@ -1203 +1191 @@
 \end{forlines}
 
-\noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\tt 46x38x1} respectively in 
+\noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\tt 46x38x1} respectively in
 the mono-processor case (i.e. global domain of {\small\tt 182x149x31}).
 An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in 
-table \autoref{tab:NC4} which compares the results of two short runs of 
-the ORCA2\_LIM reference configuration with a 4x2 mpi partitioning.
-Note the variation in the compression ratio achieved which reflects chiefly the dry to 
-wet volume ratio of each processing region.
+table \autoref{tab:NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with
+a 4x2 mpi partitioning.
+Note the variation in the compression ratio achieved which reflects chiefly the dry to wet volume ratio of
+each processing region.
 
 %------------------------------------------TABLE----------------------------------------------------
@@ -1249 +1237 @@
 %----------------------------------------------------------------------------------------------------
 
-When \key{iomput} is activated with \key{netcdf4} chunking and compression parameters for 
-fields produced via \np{iom\_put} calls are set via an equivalent and identically named namelist to 
-\textit{namnc4} in \np{xmlio\_server.def}.
-Typically this namelist serves the mean files whilst the \ngn{ namnc4} in 
-the main namelist file continues to serve the restart files.
-This duplication is unfortunate but appropriate since, if using io\_servers, the domain sizes of 
-the individual files produced by the io\_server processes may be different to those produced by 
+When \key{iomput} is activated with \key{netcdf4} chunking and compression parameters for fields produced via
+\np{iom\_put} calls are set via an equivalent and identically named namelist to \textit{namnc4} in
+\np{xmlio\_server.def}.
+Typically this namelist serves the mean files whilst the \ngn{ namnc4} in the main namelist file continues to
+serve the restart files.
+This duplication is unfortunate but appropriate since, if using io\_servers, the domain sizes of
+the individual files produced by the io\_server processes may be different to those produced by
 the invidual processing regions and different chunking choices may be desired.
  
@@ -1269 +1257 @@
 %-------------------------------------------------------------------------------------------------------------
 
-Each trend of the dynamics and/or temperature and salinity time evolution equations can be send to 
-\mdl{trddyn} and/or \mdl{trdtra} modules (see TRD directory) just after their computation ($i.e.$ at 
-the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines).
+Each trend of the dynamics and/or temperature and salinity time evolution equations can be send to
+\mdl{trddyn} and/or \mdl{trdtra} modules (see TRD directory) just after their computation
+($i.e.$ at the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines).
 This capability is controlled by options offered in \ngn{namtrd} namelist.
 Note that the output are done with xIOS, and therefore the \key{IOM} is required.
@@ -1278 +1266 @@
 
 \begin{description}
-   \item[\np{ln\_glo\_trd}]: at each \np{nn\_trd} time-step a check of the basin averaged properties of 
-                              the momentum and tracer equations is performed.
-                              This also includes a check of $T^2$, $S^2$, $\tfrac{1}{2} (u^2+v2)$, and 
-                              potential energy time evolution equations properties;
-   \item[\np{ln\_dyn\_trd}]: each 3D trend of the evolution of the two momentum components is output;
-   \item[\np{ln\_dyn\_mxl}]: each 3D trend of the evolution of the two momentum components averaged over 
-                              the mixed layer is output;
-   \item[\np{ln\_vor\_trd}]: a vertical summation of the moment tendencies is performed, then 
-                              the curl is computed to obtain the barotropic vorticity tendencies which 
-                              are output;
-   \item[\np{ln\_KE\_trd}] : each 3D trend of the Kinetic Energy equation is output ;
-   \item[\np{ln\_tra\_trd}]: each 3D trend of the evolution of temperature and salinity is output ;
-   \item[\np{ln\_tra\_mxl}]: each 2D trend of the evolution of temperature and salinity averaged over 
-                              the mixed layer is output;
+\item[\np{ln\_glo\_trd}]:
+  at each \np{nn\_trd} time-step a check of the basin averaged properties of
+  the momentum and tracer equations is performed.
+  This also includes a check of $T^2$, $S^2$, $\tfrac{1}{2} (u^2+v2)$,
+  and potential energy time evolution equations properties;
+\item[\np{ln\_dyn\_trd}]:
+  each 3D trend of the evolution of the two momentum components is output;
+\item[\np{ln\_dyn\_mxl}]:
+  each 3D trend of the evolution of the two momentum components averaged over the mixed layer is output;
+\item[\np{ln\_vor\_trd}]:
+  a vertical summation of the moment tendencies is performed,
+  then the curl is computed to obtain the barotropic vorticity tendencies which are output;
+\item[\np{ln\_KE\_trd}] :
+  each 3D trend of the Kinetic Energy equation is output;
+\item[\np{ln\_tra\_trd}]:
+  each 3D trend of the evolution of temperature and salinity is output;
+\item[\np{ln\_tra\_mxl}]:
+  each 2D trend of the evolution of temperature and salinity averaged over the mixed layer is output;
 \end{description}
 
@@ -1298 +1290 @@
 
 \textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested.
-In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl} are 
-not working, and none of the options have been tested with variable volume ($i.e.$ \key{vvl} defined).
+In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl} are not working,
+and none of the options have been tested with variable volume ($i.e.$ \key{vvl} defined).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1311 +1303 @@
 %--------------------------------------------------------------------------------------------------------------
 
-The on-line computation of floats advected either by the three dimensional velocity field or 
-constraint to remain at a given depth ($w = 0$ in the computation) have been introduced in 
-the system during the CLIPPER project.
+The on-line computation of floats advected either by the three dimensional velocity field or constraint to
+remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project.
 Options are defined by \ngn{namflo} namelis variables.
-The algorithm used is based either on the work of \cite{Blanke_Raynaud_JPO97} (default option), or 
-on a $4^th$ Runge-Hutta algorithm (\np{ln\_flork4}\forcode{ = .true.}).
-Note that the \cite{Blanke_Raynaud_JPO97} algorithm have the advantage of providing trajectories which 
+The algorithm used is based either on the work of \cite{Blanke_Raynaud_JPO97} (default option),
+or on a $4^th$ Runge-Hutta algorithm (\np{ln\_flork4}\forcode{ = .true.}).
+Note that the \cite{Blanke_Raynaud_JPO97} algorithm have the advantage of providing trajectories which
 are consistent with the numeric of the code, so that the trajectories never intercept the bathymetry.
 
 \subsubsection{Input data: initial coordinates}
 
-Initial coordinates can be given with Ariane Tools convention (IJK coordinates,
-(\np{ln\_ariane}\forcode{ = .true.}) ) or with longitude and latitude.
+Initial coordinates can be given with Ariane Tools convention
+(IJK coordinates, (\np{ln\_ariane}\forcode{ = .true.}) ) or with longitude and latitude.
 
 In case of Ariane convention, input filename is \np{init\_float\_ariane}.
@@ -1368 +1359 @@
 
 \np{jpnfl} is the total number of floats during the run.
-When initial positions are read in a restart file (\np{ln\_rstflo}\forcode{ = .true.} ), 
+When initial positions are read in a restart file (\np{ln\_rstflo}\forcode{ = .true.} ),
 \np{jpnflnewflo} can be added in the initialization file.
 
 \subsubsection{Output data}
 
-\np{nn\_writefl} is the frequency of writing in float output file and \np{nn\_stockfl} is 
-the frequency of creation of the float restart file.
+\np{nn\_writefl} is the frequency of writing in float output file and \np{nn\_stockfl} is the frequency of
+creation of the float restart file.
 
 Output data can be written in ascii files (\np{ln\_flo\_ascii}\forcode{ = .true.}).
@@ -1382 +1373 @@
 There are 2 possibilities:
 
- - if (\key{iomput}) is used, outputs are selected in  iodef.xml.
-   Here it is an example of specification to put in files description section:
+- if (\key{iomput}) is used, outputs are selected in  iodef.xml.
+Here it is an example of specification to put in files description section:
 
 \begin{xmllines}
@@ -1401 +1392 @@
  -  if (\key{iomput}) is not used, a file called \ifile{trajec\_float} will be created by IOIPSL library.
 
-See also \href{http://stockage.univ-brest.fr/~grima/Ariane/}{here} the web site describing 
-the off-line use of this marvellous diagnostic tool.
+ See also \href{http://stockage.univ-brest.fr/~grima/Ariane/}{here} the web site describing the off-line use of
+ this marvellous diagnostic tool.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1418 +1409 @@
 This on-line Harmonic analysis is actived with \key{diaharm}.
 
-Some parameters are available in namelist \ngn{namdia\_harm} :
+Some parameters are available in namelist \ngn{namdia\_harm}:
 
  - \np{nit000\_han} is the first time step used for harmonic analysis
@@ -1430 +1421 @@
  - \np{tname}       is an array with names of tidal constituents to analyse
 
-\np{nit000\_han} and \np{nitend\_han} must be between \np{nit000} and \np{nitend} of the simulation.
-The restart capability is not implemented.
-
-The Harmonic analysis solve the following equation:
+ \np{nit000\_han} and \np{nitend\_han} must be between \np{nit000} and \np{nitend} of the simulation.
+ The restart capability is not implemented.
+
+ The Harmonic analysis solve the following equation:
 
 \[h_{i} - A_{0} + \sum^{nb\_ana}_{j=1}[A_{j}cos(\nu_{j}t_{j}-\phi_{j})] = e_{i}\]
 
-With $A_{j}$, $\nu_{j}$, $\phi_{j}$, the amplitude, frequency and phase for each wave and 
-$e_{i}$ the error.
+With $A_{j}$, $\nu_{j}$, $\phi_{j}$, the amplitude, frequency and phase for each wave and $e_{i}$ the error.
 $h_{i}$ is the sea level for the time $t_{i}$ and $A_{0}$ is the mean sea level. \\
 We can rewrite this equation:
@@ -1459 +1449 @@
 %-------------------------------------------------------------------------------------------------------------
 
-A module is available to compute the transport of volume, heat and salt through sections. 
+A module is available to compute the transport of volume, heat and salt through sections.
 This diagnostic is actived with \key{diadct}.
 
 Each section is defined by the coordinates of its 2 extremities.
-The pathways between them are contructed using tools which can be found in 
-\texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT} 
-and are written in a binary file 
-\texttt{section\_ijglobal.diadct\_ORCA2\_LIM} 
-which is later read in by NEMO to compute on-line transports.
+The pathways between them are contructed using tools which can be found in \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT}
+and are written in a binary file \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is later read in by
+NEMO to compute on-line transports.
 
 The on-line transports module creates three output ascii files:
@@ -1477 +1465 @@
 - \texttt{salt\_transport}   for   salt transports (unit: $10^{9}Kg s^{-1}$) \\
 
-Namelist variables in \ngn{namdct} control how frequently the flows are summed and the time scales over 
-which they are averaged, as well as the level of output for debugging:
+Namelist variables in \ngn{namdct} control how frequently the flows are summed and the time scales over which
+they are averaged, as well as the level of output for debugging:
 \np{nn\_dct}   : frequency of instantaneous transports computing
 \np{nn\_dctwri}: frequency of writing ( mean of instantaneous transports )
@@ -1485 +1473 @@
 \subsubsection{Creating a binary file containing the pathway of each section}
 
-In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, 
-the file \textit{ {list\_sections.ascii\_global}} contains a list of all the sections that are to 
-be computed (this list of sections is based on MERSEA project metrics).
+In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run},
+the file \textit{ {list\_sections.ascii\_global}} contains a list of all the sections that are to be computed
+(this list of sections is based on MERSEA project metrics).
 
 Another file is available for the GYRE configuration (\texttt{ {list\_sections.ascii\_GYRE}}).
@@ -1505 +1493 @@
  - \texttt{ice}       to compute surface and volume ice transports, \texttt{noice}     if no. \\
 
-\noindent The results of the computing of transports, and the directions of positive and 
-negative flow do not depend on the order of the 2 extremities in this file. \\
-
-\noindent If nclass $\neq$ 0,the next lines contain the class type and the nclass bounds: \\
+ \noindent The results of the computing of transports, and the directions of positive and
+ negative flow do not depend on the order of the 2 extremities in this file. \\
+
+\noindent If nclass $\neq$ 0, the next lines contain the class type and the nclass bounds: \\
 {\scriptsize \texttt{
 long1 lat1 long2 lat2 nclass (ok/no)strpond (no)ice section\_name \\
@@ -1531 +1519 @@
  - \texttt{zsigp} for potential density classes \\
   
-The script \texttt{job.ksh} computes the pathway for each section and 
-creates a binary file \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is read by NEMO. \\
-
-It is possible to use this tools for new configuations: \texttt{job.ksh} has to be updated with 
-the coordinates file name and path. \\
-
-Examples of two sections, the ACC\_Drake\_Passage with no classes, and the ATL\_Cuba\_Florida with 
-4 temperature clases (5 class bounds), are shown: \\
+ The script \texttt{job.ksh} computes the pathway for each section and creates a binary file
+ \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is read by NEMO. \\
+
+ It is possible to use this tools for new configuations: \texttt{job.ksh} has to be updated with
+ the coordinates file name and path. \\
+
+ Examples of two sections, the ACC\_Drake\_Passage with no classes,
+ and the ATL\_Cuba\_Florida with 4 temperature clases (5 class bounds), are shown: \\
 \noindent {\scriptsize \texttt{
 -68.    -54.5   -60.    -64.7  00 okstrpond noice ACC\_Drake\_Passage \\
@@ -1559 +1547 @@
 transport\_total}}                                     \\
 
-For sections with classes, the first \texttt{nclass-1} lines correspond to the transport for 
-each class and the last line corresponds to the total transport summed over all classes.
-For sections with no classes, class number \texttt{1} corresponds to \texttt{total class} and 
+For sections with classes, the first \texttt{nclass-1} lines correspond to the transport for each class and
+the last line corresponds to the total transport summed over all classes.
+For sections with no classes, class number \texttt{1} corresponds to \texttt{total class} and
 this class is called \texttt{N}, meaning \texttt{none}.
 
@@ -1568 +1556 @@
 - \texttt{transport\_direction2} is the negative part of the transport ($\leq$ 0). \\
 
-\noindent The \texttt{section slope coefficient} gives information about the significance of 
-transports signs and direction: \\
+\noindent The \texttt{section slope coefficient} gives information about the significance of transports signs and
+direction: \\
 
 \begin{table} \scriptsize
@@ -1591 +1579 @@
 
 
-Changes in steric sea level are caused when changes in the density of the water column imply an 
-expansion or contraction of the column.
-It is essentially produced through surface heating/cooling and to a lesser extent through 
-non-linear effects of the equation of state (cabbeling, thermobaricity...).
+Changes in steric sea level are caused when changes in the density of the water column imply an expansion or
+contraction of the column.
+It is essentially produced through surface heating/cooling and to a lesser extent through non-linear effects of
+the equation of state (cabbeling, thermobaricity...).
 Non-Boussinesq models contain all ocean effects within the ocean acting on the sea level.
 In particular, they include the steric effect.
-In contrast, Boussinesq models, such as \NEMO, conserve volume, rather than mass, 
+In contrast, Boussinesq models, such as \NEMO, conserve volume, rather than mass,
 and so do not properly represent expansion or contraction.
 The steric effect is therefore not explicitely represented.
-This approximation does not represent a serious error with respect to the flow field calculated by 
-the model \citep{Greatbatch_JGR94}, but extra attention is required when investigating sea level, 
-as steric changes are an important contribution to local changes in sea level on seasonal and 
-climatic time scales.
-This is especially true for investigation into sea level rise due to global warming. 
-
-Fortunately, the steric contribution to the sea level consists of a spatially uniform component that 
+This approximation does not represent a serious error with respect to the flow field calculated by the model
+\citep{Greatbatch_JGR94}, but extra attention is required when investigating sea level,
+as steric changes are an important contribution to local changes in sea level on seasonal and climatic time scales.
+This is especially true for investigation into sea level rise due to global warming.
+
+Fortunately, the steric contribution to the sea level consists of a spatially uniform component that
 can be diagnosed by considering the mass budget of the world ocean \citep{Greatbatch_JGR94}.
-In order to better understand how global mean sea level evolves and thus how the steric sea level can 
-be diagnosed, we compare, in the following, the non-Boussinesq and Boussinesq cases.
+In order to better understand how global mean sea level evolves and thus how the steric sea level can be diagnosed,
+we compare, in the following, the non-Boussinesq and Boussinesq cases.
 
 Let denote
@@ -1628 +1615 @@
 \]
 
-Temporal changes in total mass is obtained from the density conservation equation :
+Temporal changes in total mass is obtained from the density conservation equation:
 
 \[ \frac{1}{e_3} \partial_t ( e_3\,\rho) + \nabla( \rho \, \textbf{U} ) 
@@ -1634 +1621 @@
  \label{eq:Co_nBq} \]
 
-where $\rho$ is the \textit{in situ} density, and \textit{emp} the surface mass 
-exchanges with the other media of the Earth system (atmosphere, sea-ice, land). 
+where $\rho$ is the \textit{in situ} density, and \textit{emp} the surface mass exchanges with the other media of
+the Earth system (atmosphere, sea-ice, land).
 Its global averaged leads to the total mass change 
 
@@ -1642 +1629 @@
 
 where $\overline{\textit{emp}} = \int_S \textit{emp}\,ds$ is the net mass flux through the ocean surface.
-Bringing \autoref{eq:Mass_nBq} and the time derivative of \autoref{eq:MV_nBq} together leads to 
+Bringing \autoref{eq:Mass_nBq} and the time derivative of \autoref{eq:MV_nBq} together leads to
 the evolution equation of the mean sea level
 
@@ -1649 +1636 @@
  \label{eq:ssh_nBq} \]
 
-The first term in equation \autoref{eq:ssh_nBq} alters sea level by adding or 
-subtracting mass from the ocean. 
-The second term arises from temporal changes in the global mean density; $i.e.$ from steric effects. 
-
-In a Boussinesq fluid, $\rho$ is replaced by $\rho_o$ in all the equation except when 
-$\rho$ appears multiplied by the gravity ($i.e.$ in the hydrostatic balance of the primitive Equations).
-In particular, the mass conservation equation, \autoref{eq:Co_nBq}, degenerates into 
-the incompressibility equation:
+The first term in equation \autoref{eq:ssh_nBq} alters sea level by adding or subtracting mass from the ocean. 
+The second term arises from temporal changes in the global mean density; $i.e.$ from steric effects.
+
+In a Boussinesq fluid, $\rho$ is replaced by $\rho_o$ in all the equation except when $\rho$ appears multiplied by
+the gravity ($i.e.$ in the hydrostatic balance of the primitive Equations).
+In particular, the mass conservation equation, \autoref{eq:Co_nBq}, degenerates into the incompressibility equation:
 
 \[ \frac{1}{e_3} \partial_t ( e_3 ) + \nabla( \textbf{U} )
@@ -1667 +1652 @@
  \label{eq:V_Bq} \]
 
-Only the volume is conserved, not mass, or, more precisely, the mass which is conserved is 
-the Boussinesq mass, $\mathcal{M}_o = \rho_o \mathcal{V}$.
-The total volume (or equivalently the global mean sea level) is altered only by net volume fluxes across 
-the ocean surface, not by changes in mean mass of the ocean: the steric effect is missing in 
-a Boussinesq fluid.
-
-Nevertheless, following \citep{Greatbatch_JGR94}, the steric effect on the volume can be diagnosed by 
+Only the volume is conserved, not mass, or, more precisely, the mass which is conserved is the Boussinesq mass,
+$\mathcal{M}_o = \rho_o \mathcal{V}$.
+The total volume (or equivalently the global mean sea level) is altered only by net volume fluxes across
+the ocean surface, not by changes in mean mass of the ocean: the steric effect is missing in a Boussinesq fluid.
+
+Nevertheless, following \citep{Greatbatch_JGR94}, the steric effect on the volume can be diagnosed by
 considering the mass budget of the ocean. 
-The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by 
-surface mass flux must be compensated by a spatially uniform change in the mean sea level due to 
-expansion/contraction of the ocean \citep{Greatbatch_JGR94}.
-In others words, the Boussinesq mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$, 
-the  total mass of the ocean seen by the Boussinesq model, via the steric contribution to the sea level, 
+The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by surface mass flux
+must be compensated by a spatially uniform change in the mean sea level due to expansion/contraction of the ocean
+\citep{Greatbatch_JGR94}.
+In others words, the Boussinesq mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$,
+the total mass of the ocean seen by the Boussinesq model, via the steric contribution to the sea level,
 $\eta_s$, a spatially uniform variable, as follows:
 
@@ -1685 +1669 @@
  \label{eq:M_Bq} \]
 
-Any change in $\mathcal{M}$ which cannot be explained by the net mass flux through 
-the ocean surface is converted into a mean change in sea level.
-Introducing the total density anomaly, $\mathcal{D}= \int_D d_a \,dv$, where 
-$d_a = (\rho -\rho_o ) / \rho_o$ is the density anomaly used in \NEMO (cf. \autoref{subsec:TRA_eos}) in 
-\autoref{eq:M_Bq} leads to a very simple form for the steric height:
+Any change in $\mathcal{M}$ which cannot be explained by the net mass flux through the ocean surface
+is converted into a mean change in sea level.
+Introducing the total density anomaly, $\mathcal{D}= \int_D d_a \,dv$,
+where $d_a = (\rho -\rho_o ) / \rho_o$ is the density anomaly used in \NEMO (cf. \autoref{subsec:TRA_eos})
+in \autoref{eq:M_Bq} leads to a very simple form for the steric height:
 
 \[ \eta_s = - \frac{1}{\mathcal{A}} \mathcal{D}
@@ -1699 +1683 @@
 We do not recommend that.
 Indeed, in this case $\rho_o$ depends on the initial state of the ocean.
-Since $\rho_o$ has a direct effect on the dynamics of the ocean (it appears in 
-the pressure gradient term of the momentum equation) it is definitively not a good idea when 
-inter-comparing experiments.
+Since $\rho_o$ has a direct effect on the dynamics of the ocean
+(it appears in the pressure gradient term of the momentum equation)
+it is definitively not a good idea when inter-comparing experiments.
 We better recommend to fixe once for all $\rho_o$ to $1035\;Kg\,m^{-3}$.
-This value is a sensible choice for the reference density used in a Boussinesq ocean climate model since, 
-with the exception of only a small percentage of the ocean, density in the World Ocean varies by 
-no more than 2$\%$ from this value (\cite{Gill1982}, page 47).
-
-Second, we have assumed here that the total ocean surface, $\mathcal{A}$, does not change when 
-the sea level is changing as it is the case in all global ocean GCMs 
+This value is a sensible choice for the reference density used in a Boussinesq ocean climate model since,
+with the exception of only a small percentage of the ocean, density in the World Ocean varies by no more than
+2$\%$ from this value (\cite{Gill1982}, page 47).
+
+Second, we have assumed here that the total ocean surface, $\mathcal{A}$,
+does not change when the sea level is changing as it is the case in all global ocean GCMs
 (wetting and drying of grid point is not allowed).
   
-Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of 
-free surface which is considered.
+Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of free surface which is considered.
 In the non linear free surface case, $i.e.$ \key{vvl} defined, it is given by
 
@@ -1719 +1702 @@
  \label{eq:discrete_steric_Bq_nfs} \]
 
-whereas in the linear free surface, the volume above the \textit{z=0} surface must be explicitly 
-taken into account to better approximate the total ocean mass and thus the steric sea level:
+whereas in the linear free surface,
+the volume above the \textit{z=0} surface must be explicitly taken into account to
+better approximate the total ocean mass and thus the steric sea level:
 
 \[ \eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t}e_{2t}e_{3t} + \sum_{i,\,j} d_a\; e_{1t}e_{2t} \eta }
@@ -1729 +1713 @@
 In the real ocean, sea ice (and snow above it)  depresses the liquid seawater through its mass loading.
 This depression is a result of the mass of sea ice/snow system acting on the liquid ocean.
-There is, however, no dynamical effect associated with these depressions in the liquid ocean sea level, 
+There is, however, no dynamical effect associated with these depressions in the liquid ocean sea level,
 so that there are no associated ocean currents.
-Hence, the dynamically relevant sea level is the effective sea level, $i.e.$ the sea level as if 
-sea ice (and snow) were converted to liquid seawater \citep{Campin_al_OM08}.
-However, in the current version of \NEMO the sea-ice is levitating above the ocean without 
-mass exchanges between ice and ocean.
-Therefore the model effective sea level is always given by $\eta + \eta_s$, 
-whether or not there is sea ice present.
+Hence, the dynamically relevant sea level is the effective sea level,
+$i.e.$ the sea level as if sea ice (and snow) were converted to liquid seawater \citep{Campin_al_OM08}.
+However, in the current version of \NEMO the sea-ice is levitating above the ocean without mass exchanges between
+ice and ocean.
+Therefore the model effective sea level is always given by $\eta + \eta_s$, whether or not there is sea ice present.
 
 In AR5 outputs, the thermosteric sea level is demanded.
@@ -1747 +1730 @@
 where $S_o$ and $p_o$ are the initial salinity and pressure, respectively.
 
-Both steric and thermosteric sea level are computed in \mdl{diaar5} which needs
-the \key{diaar5} defined to be called.
+Both steric and thermosteric sea level are computed in \mdl{diaar5} which needs the \key{diaar5} defined to
+be called.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1761 +1744 @@
 \subsection{Depth of various quantities (\protect\mdl{diahth})}
 
-Among the available diagnostics the following ones are obtained when defining the \key{diahth} CPP key: 
+Among the available diagnostics the following ones are obtained when defining the \key{diahth} CPP key:
 
 - the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
@@ -1782 +1765 @@
 %-----------------------------------------------------------------------------------------
 
-The poleward heat and salt transports, their advective and diffusive component, and 
-the meriodional stream function can be computed on-line in \mdl{diaptr} \np{ln\_diaptr} to true 
+The poleward heat and salt transports, their advective and diffusive component,
+and the meriodional stream function can be computed on-line in \mdl{diaptr} \np{ln\_diaptr} to true
 (see the \textit{\ngn{namptr} } namelist below).
-When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, 
-Indian, Pacific and Indo-Pacific Oceans (defined north of 30\deg S) as well as for the World Ocean.
-The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains 
-three 2D mask arrays, the Indo-Pacific mask been deduced from the sum of the Indian and 
-Pacific mask (\autoref{fig:mask_subasins}).
+When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, Indian,
+Pacific and Indo-Pacific Oceans (defined north of 30\deg S) as well as for the World Ocean.
+The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays,
+the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}).
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \begin{center}
-   \includegraphics[width=1.0\textwidth]{Fig_mask_subasins}
-   \caption{   \protect\label{fig:mask_subasins}
-      Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute the heat and 
-      salt transports as well as the meridional stream-function: Atlantic basin (red), 
-      Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green).
-      Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay are removed from 
-      the sub-basins. Note also that the Arctic Ocean has been split into Atlantic and 
-      Pacific basins along the North fold line.}
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=1.0\textwidth]{Fig_mask_subasins}
+    \caption{  \protect\label{fig:mask_subasins}
+      Decomposition of the World Ocean (here ORCA2) into sub-basin used in to
+      compute the heat and salt transports as well as the meridional stream-function:
+      Atlantic basin (red), Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green).
+      Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay are removed from the sub-basins.
+      Note also that the Arctic Ocean has been split into Atlantic and Pacific basins along the North fold line.}
 \end{center} \end{figure}  
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -1827 +1809 @@
 The 25 hour mean is available for daily runs by summing up the 25 hourly instantananeous hourly values from
 midnight at the start of the day to midight at the day end.
-This diagnostic is actived with the logical  $ln\_dia25h$
+This diagnostic is actived with the logical $ln\_dia25h$.
 
 % -----------------------------------------------------------
@@ -1839 +1821 @@
 %----------------------------------------------------------------------------------------------------------
 
-A module is available to output the surface (top), mid water and bed diagnostics of a set of 
-standard variables.
-This can be a useful diagnostic when hourly or sub-hourly output is required in 
-high resolution tidal outputs.
+A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables.
+This can be a useful diagnostic when hourly or sub-hourly output is required in high resolution tidal outputs.
 The tidal signal is retained but the overall data usage is cut to just three vertical levels.
 Also the bottom level is calculated for each cell.
-This diagnostic is actived with the logical  $ln\_diatmb$
+This diagnostic is actived with the logical $ln\_diatmb$.
 
 % -----------------------------------------------------------
@@ -1859 +1839 @@
 
 in the zonal, meridional and vertical directions respectively.
-The vertical component is included although it is not strictly valid as the vertical velocity is 
-calculated from the continuity equation rather than as a prognostic variable.
+The vertical component is included although it is not strictly valid as the vertical velocity is calculated from
+the continuity equation rather than as a prognostic variable.
 Physically this represents the rate at which information is propogated across a grid cell.
-Values greater than 1 indicate that information is propagated across more than one grid cell in 
-a single time step.
-
-The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in 
-the \ngn{namctl} namelist.
+Values greater than 1 indicate that information is propagated across more than one grid cell in a single time step.
+
+The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist.
 The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii.
-In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with 
-the coordinates of where the maximum value occurs.
-At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is 
-printed along with the coordinates of each.
+In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of
+where the maximum value occurs.
+At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along
+with the coordinates of each.
 The maximum values from the run are also copied to the ocean.output file. 
 

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex

@@ -18 +18 @@
 The skin temperature can be split into three parts:
 \begin{itemize}
-\item A foundation SST which is free from diurnal warming.
-\item A warm layer, typically ~3\,m thick, where heating from solar radiation can
-cause a warm stably stratified layer during the daytime
-\item A cool skin, a thin layer, approximately ~1\,mm thick, where long wave cooling
-is dominant and cools the immediate ocean surface.
+\item
+  A foundation SST which is free from diurnal warming.
+\item
+  A warm layer, typically ~3\,m thick,
+  where heating from solar radiation can cause a warm stably stratified layer during the daytime
+\item
+  A cool skin, a thin layer, approximately ~1\, mm thick,
+  where long wave cooling is dominant and cools the immediate ocean surface.
 \end{itemize}
 
-Models are provided for both the warm layer, \mdfl{diurnal_bulk}, and the cool skin,
-\mdl{cool_skin}.  Foundation SST is not considered as it can be obtained
-either from the main NEMO model ($i.e.$ from the temperature of the top few model levels)
-or from some other source.  
-It must be noted that both the cool skin and warm layer models produce estimates of 
-the change in temperature ($\Delta T_{\rm{cs}}$ and $\Delta T_{\rm{wl}}$) 
-and both must be added to a foundation SST to obtain the true skin temperature.
+Models are provided for both the warm layer, \mdfl{diurnal_bulk}, and the cool skin, \mdl{cool_skin}.
+Foundation SST is not considered as it can be obtained either from the main NEMO model
+($i.e.$ from the