Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex	(revision 10501)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex	(revision 10502)
@@ -3,5 +3,5 @@
 \begin{document}
 % ================================================================
-% Chapter 2 Space and Time Domain (DOM)
+% Chapter 2 ——— Space and Time Domain (DOM)
 % ================================================================
 \chapter{Space Domain (DOM)}
@@ -40,5 +40,5 @@
 \begin{figure}[!tb]
   \begin{center}
-    \includegraphics[width=0.90\textwidth]{Fig_cell}
+    \includegraphics[]{Fig_cell}
     \caption{
       \protect\label{fig:cell}
@@ -46,6 +46,6 @@
       $t$ indicates scalar points where temperature, salinity, density, pressure and
       horizontal divergence are defined.
-      ($u$,$v$,$w$) indicates vector points,
-      and $f$ indicates vorticity points where both relative and planetary vorticities are defined
+      $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and
+      planetary vorticities are defined.
     }
   \end{center}
@@ -64,23 +64,21 @@
 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points.
 
-The ocean mesh (\ie the position of all the scalar and vector points) is defined by
-the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$.
+The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that
+gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}.
 In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of
 the grid-point where the scale factors are defined.
 Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}.
-As a result,
-the mesh on which partial derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$,
-and $\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity.
-Discrete partial derivatives are formulated by the traditional,
-centred second order finite difference approximation while
-the scale factors are chosen equal to their local analytical value.
+As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and
+$\pd[]{z}$ are evaluated in a uniform mesh with a grid size of unity.
+Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation
+while the scale factors are chosen equal to their local analytical value.
 An important point here is that the partial derivative of the scale factors must be evaluated by
 centred finite difference approximation, not from their analytical expression.
-This preserves the symmetry of the discrete set of equations and
-therefore satisfies many of the continuous properties (see \autoref{apdx:C}).
+This preserves the symmetry of the discrete set of equations and therefore satisfies many of
+the continuous properties (see \autoref{apdx:C}).
 A similar, related remark can be made about the domain size:
 when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors
-(see \autoref{eq:DOM_bar}) in the next section).
+(see \autoref{eq:DOM_bar} in the next section).
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -89,17 +87,25 @@
     \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|}
       \hline
-      T	&$i$ 		& $j$		& $k$		 \\ \hline
-      u	& $i+1/2$	& $j$		& $k$	 	\\ \hline
-      v	& $i$	 	& $j+1/2$	& $k$	 	\\ \hline
-      w	& $i$	 	& $j$		& $k+1/2$ 	\\ \hline
-      f	& $i+1/2$ 	& $j+1/2$	& $k$	 	\\ \hline
-      uw	& $i+1/2$ 	& $j$		& $k+1/2$ 	\\ \hline
-      vw	& $i$	 	& $j+1/2$	& $k+1/2$ 	\\ \hline
-      fw	& $i+1/2$ 	& $j+1/2$	& $k+1/2$ 	\\ \hline
+      T	& $i 		 $ & $j      $ & $k      $ \\
+      \hline
+      u	& $i + 1/2$ & $j      $ & $k      $ \\
+      \hline
+      v	& $i	 	 $ & $j + 1/2$ & $k      $ \\
+      \hline
+      w	& $i	 	 $ & $j      $ & $k + 1/2$ \\
+      \hline
+      f	& $i + 1/2$ & $j + 1/2$ & $k      $ \\
+      \hline
+      uw	& $i + 1/2$ & $j      $ & $k + 1/2$ \\
+      \hline
+      vw	& $i	 	 $ & $j + 1/2$ & $k + 1/2$ \\
+      \hline
+      fw	& $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\
+      \hline
     \end{tabular}
     \caption{
       \protect\label{tab:cell}
       Location of grid-points as a function of integer or integer and a half value of the column, line or level.
-      This indexing is only used for the writing of the semi-discrete equation.
+      This indexing is only used for the writing of the semi -discrete equation.
       In the code, the indexing uses integer values only and has a reverse direction in the vertical
       (see \autoref{subsec:DOM_Num_Index})
@@ -115,92 +121,89 @@
 \label{subsec:DOM_operators}
 
-Given the values of a variable $q$ at adjacent points,
-the differencing and averaging operators at the midpoint between them are:
-\[
+Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at
+the midpoint between them are:
+\begin{alignat*}{2}
   % \label{eq:di_mi}
-  \begin{split}
-    \delta_i [q]       &=  \  \    q(i+1/2)  - q(i-1/2)		\\
-    \overline q^{\,i} &= \left\{ q(i+1/2) + q(i-1/2) \right\} \; / \; 2
-  \end{split}
-\]
-
-Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and $k+1/2$.
+  \delta_i [q]      &= &       &q (i + 1/2) - q (i - 1/2) \\
+  \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2
+\end{alignat*}
+
+Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$.
 Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at
 a $t$-point has its three components defined at $u$-, $v$- and $w$-points while
-its Laplacien is defined at $t$-point.
-These operators have the following discrete forms in the curvilinear $s$-coordinate system:
+its Laplacian is defined at $t$-point.
+These operators have the following discrete forms in the curvilinear $s$-coordinates system:
 \[
   % \label{eq:DOM_grad}
-  \nabla q\equiv 	\frac{1}{e_{1u} } \delta_{i+1/2 } [q] \;\,\mathbf{i}
-  +	\frac{1}{e_{2v} } \delta_{j+1/2 } [q] \;\,\mathbf{j}
-  +	\frac{1}{e_{3w}} \delta_{k+1/2} [q] \;\,\mathbf{k}
+  \nabla q \equiv   \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i
+                  + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j
+                  + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k
 \]
 \begin{multline*}
   % \label{eq:DOM_lap}
-  \Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
-  \;\left(          \delta_i  \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right]
-    +                        \delta_j  \left[ \frac{e_{1v}\,e_{3v}}  {e_{2v}} \;\delta_{j+1/2} [q] \right] \;  \right)		\\
-  +\frac{1}{e_{3t}} \delta_k \left[ \frac{1}{e_{3w} }                     \;\delta_{k+1/2} [q] \right]
+  \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}}
+                    \; \lt[   \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt)
+                            + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] \\
+                  + \frac{1}{e_{3t}}
+                              \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt]
 \end{multline*}
 
-Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$
-defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points,
-and its divergence defined at $t$-points:
-\begin{align*}
-  % \label{eq:DOM_curl}
-  \nabla \times {\rm{\bf A}}\equiv &
-                                     \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right)  &\ \mathbf{i} \\
-  +& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1  \right] -\delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right)  &\ \mathbf{j} \\
-  +& \frac{1}{e_{1f} \,e_{2f}    } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2  \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right)  &\ \mathbf{k}
-\end{align*}
-\begin{align*}
-  % \label{eq:DOM_div}
-  \nabla \cdot \rm{\bf A} \equiv
-  \frac{1}{e_{1t}\,e_{2t}\,e_{3t}} \left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right]
-  +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right]
-\end{align*}
+Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at
+vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and
+its divergence defined at $t$-points:
+\begin{multline}
+% \label{eq:DOM_curl}
+  \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}}
+                                 \Big[   \delta_{j + 1/2} (e_{3w} \, a_3)
+                                       - \delta_{k + 1/2} (e_{2v} \, a_2) \Big] \vect i \\
+                               + \frac{1}{e_{2u} \, e_{3uw}}
+                                 \Big[   \delta_{k + 1/2} (e_{1u} \, a_1)
+                                       - \delta_{i + 1/2} (e_{3w} \, a_3) \Big] \vect j \\
+                               + \frac{1}{e_{1f} \, e_{2f}}
+                                 \Big[   \delta_{i + 1/2} (e_{2v} \, a_2)
+                                       - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k
+\end{multline}
+\begin{equation}
+% \label{eq:DOM_div}
+  \nabla \cdot \vect A \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}}
+                                \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big]
+                              + \frac{1}{e_{3t}} \delta_k (a_3)
+\end{equation}
 
 The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which
-is a masked field (\ie equal to zero inside solid area):
+is a masked field (i.e. equal to zero inside solid area):
 \begin{equation}
   \label{eq:DOM_bar}
-  \bar q 	=         \frac{1}{H}    \int_{k^b}^{k^o} {q\;e_{3q} \,dk}
-  \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} }
+  \bar q = \frac{1}{H} \int_{k^b}^{k^o} q \; e_{3q} \, dk \equiv \frac{1}{H_q} \sum \limits_k q \; e_{3q}
 \end{equation}
 where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points,
-$k^b$ and $k^o$ are the bottom and surface $k$-indices,
-and the symbol $k^o$ refers to a summation over all grid points of the same type in the direction indicated by
-the subscript (here $k$). 
+$k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over
+all grid points of the same type in the direction indicated by the subscript (here $k$).
 
 In continuous form, the following properties are satisfied:
-\begin{equation}
+\begin{gather}
   \label{eq:DOM_curl_grad}
-  \nabla \times \nabla q ={\rm {\bf {0}}}
-\end{equation}
-\begin{equation}
+  \nabla \times \nabla q = \vect 0 \\
   \label{eq:DOM_div_curl}
-  \nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0
-\end{equation}
+  \nabla \cdot (\nabla \times \vect A) = 0
+\end{gather}
 
 It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as
-the scalar $q$ is taken at $t$-points and
-the vector \textbf{A} has its components defined at vector points $(u,v,w)$.
+the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at
+vector points $(u,v,w)$.
 
 Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area.
-Using integration by parts it can be shown that
-the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators,
-and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$, $\overline{\,\cdot\,}^{\,k}$ and
-$\overline{\,\cdot\,}^{\,k}$) are symmetric linear operators, \ie
-\begin{align}
+Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$)
+are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$,
+$\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie
+\begin{alignat}{4}
   \label{eq:DOM_di_adj}
-  \sum\limits_i { a_i \;\delta_i \left[ b \right]}
-  &\equiv -\sum\limits_i {\delta_{i+1/2} \left[ a \right]\;b_{i+1/2} }      \\
+  &\sum \limits_i a_i \; \delta_i [b]      &\equiv &- &&\sum \limits_i \delta      _{   i + 1/2} [a] &b_{i + 1/2} \\
   \label{eq:DOM_mi_adj}
-  \sum\limits_i { a_i \;\overline b^{\,i}}
-  & \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} }
-\end{align}
-
-In other words, the adjoint of the differencing and averaging operators are $\delta_i^*=\delta_{i+1/2}$ and 
-${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively. 
+  &\sum \limits_i a_i \; \overline b^{\, i} &\equiv &  &&\sum \limits_i \overline a ^{\, i + 1/2}     &b_{i + 1/2}
+\end{alignat}
+
+In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 
+$(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively.
 These two properties will be used extensively in the \autoref{apdx:C} to
 demonstrate integral conservative properties of the discrete formulation chosen.
@@ -215,5 +218,5 @@
 \begin{figure}[!tb]
   \begin{center}
-    \includegraphics[width=0.90\textwidth]{Fig_index_hor}
+    \includegraphics[]{Fig_index_hor}
     \caption{
       \protect\label{fig:index_hor}
@@ -230,5 +233,5 @@
 Therefore a specific integer indexing must be defined for points other than $t$-points
 (\ie velocity and vorticity grid-points).
-Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k=1$.
+Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k = 1$.
 
 % -----------------------------------
@@ -250,10 +253,10 @@
 \label{subsec:DOM_Num_Index_vertical}
 
-In the vertical, the chosen indexing requires special attention since
-the $k$-axis is re-orientated downward in the \fortran code compared to
-the indexing used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}.
-The sea surface corresponds to the $w$-level $k=1$ which is the same index as $t$-level just below
+In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in
+the \fortran code compared to the indexing used in the semi -discrete equations and
+given in \autoref{subsec:DOM_cell}.
+The sea surface corresponds to the $w$-level $k = 1$ which is the same index as $t$-level just below
 (\autoref{fig:index_vert}).
-The last $w$-level ($k=jpk$) either corresponds to the ocean floor or is inside the bathymetry while
+The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while
 the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}).
 Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index,
@@ -262,11 +265,12 @@
 have the same $i$ or $j$ index
 (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}).
-Since the scale factors are chosen to be strictly positive, a \emph{minus sign} appears in the \fortran 
-code \emph{before all the vertical derivatives} of the discrete equations given in this documentation.
+Since the scale factors are chosen to be strictly positive,
+a \textit{minus sign} appears in the \fortran code \textit{before all the vertical derivatives} of
+the discrete equations given in this documentation.
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!pt]
   \begin{center}
-    \includegraphics[width=.90\textwidth]{Fig_index_vert}
+    \includegraphics[]{Fig_index_vert}
     \caption{
       \protect\label{fig:index_vert}
@@ -287,7 +291,4 @@
 The total size of the computational domain is set by the parameters \np{jpiglo},
 \np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively.
-%%%
-%%%
-%%%
 Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when
 the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined,
@@ -299,37 +300,39 @@
 \section{Needed fields}
 \label{sec:DOM_fields}
-The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
+The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that
+gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}.
 The associated scale factors are defined using the analytical first derivative of the transformation
 \autoref{eq:scale_factors}.
-Necessary fields for configuration definition are: \\
-Geographic position :
-
-longitude: glamt, glamu, glamv and glamf (at T, U, V and F point)
-
-latitude: gphit, gphiu, gphiv and gphif (at T, U, V and F point)\\
-Coriolis parameter (if domain not on the sphere): 
-
- ff\_f  and  ff\_t (at T and F point)\\
-Scale factors : 
+Necessary fields for configuration definition are:
+
+\begin{itemize}
+\item
+  Geographic position:
+  longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and
+  latitude  with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif}
+  (all respectively at T, U, V and F point)
+\item
+  Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t}
+  (at T and F point)
+\item
+  Scale factors:
+  \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction),
+  \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and
+  \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\
+  \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits), 
+  \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed.
+\end{itemize}
  
- e1t, e1u, e1v and e1f (on i direction),
-
- e2t, e2u, e2v and e2f (on j direction) and
-
- ie1e2u\_v, e1e2u , e1e2v   
- 
-e1e2u , e1e2v are u and v surfaces (if gridsize reduction in some straits)\\
-ie1e2u\_v is a flag to flag set u and  v surfaces are neither read nor computed.\\
- 
-These fields can be read in an domain input file which name is setted in
-\np{cn\_domcfg} parameter specified in \ngn{namcfg}.
+These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in
+\ngn{namcfg}.
 
 \nlst{namcfg}
-or they can be defined in an analytical way in MY\_SRC directory of the configuration.
+
+Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration.
 For Reference Configurations of NEMO input domain files are supplied by NEMO System Team.
-For analytical definition of input fields two routines are supplied: \mdl{userdef\_hgr} and \mdl{userdef\_zgr}.
-They are an example of GYRE configuration parameters, and they are available in NEMO/OPA\_SRC/USR directory,
-they provide the horizontal and vertical mesh. 
+For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}.
+They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory,
+they provide the horizontal and vertical mesh.
 % -------------------------------------------------------------------------------------------------------------
 %        Needed fields 
@@ -366,26 +369,28 @@
 ($i$ and $j$, respectively) (geographical configuration of the mesh),
 the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$,
-and their derivatives $\lambda'(i)$ $\varphi'(j)$ in the \mdl{domhgr} module.
+and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module.
 The model computes the grid-point positions and scale factors in the horizontal plane as follows:
-\begin{flalign*}
-  \lambda_t &\equiv \text{glamt}= \lambda(i)	  & \varphi_t &\equiv \text{gphit} = \varphi(j)\\
-  \lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\
-  \lambda_v &\equiv \text{glamv}= \lambda(i)       & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\
-  \lambda_f &\equiv \text{glamf }= \lambda(i+1/2)& \varphi_f &\equiv \text{gphif }= \varphi(j+1/2)
-\end{flalign*}
-\begin{flalign*}
-  e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i)		\; \cos\varphi(j)  |&
-  e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\
-  e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)	\; \cos\varphi(j)  |&
-  e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\
-  e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i)		\; \cos\varphi(j+1/2)  |&
-  e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|\\
-  e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)\; \cos\varphi(j+1/2)  |&
-  e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|
-\end{flalign*}
+\begin{align*}
+   \lambda_t &\equiv \text{glamt} =      \lambda (i      )
+  &\varphi_t &\equiv \text{gphit} =      \varphi (j      ) \\
+   \lambda_u &\equiv \text{glamu} =      \lambda (i + 1/2)
+  &\varphi_u &\equiv \text{gphiu} =      \varphi (j      ) \\
+   \lambda_v &\equiv \text{glamv} =      \lambda (i      )
+  &\varphi_v &\equiv \text{gphiv} =      \varphi (j + 1/2) \\
+   \lambda_f &\equiv \text{glamf} =      \lambda (i + 1/2)
+  &\varphi_f &\equiv \text{gphif} =      \varphi (j + 1/2) \\
+   e_{1t}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j      ) |
+  &e_{2t}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\
+   e_{1u}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j      ) |
+  &e_{2u}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\
+   e_{1v}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j + 1/2) |
+  &e_{2v}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | \\
+   e_{1f}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) |
+  &e_{2f}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         |
+\end{align*}
 where the last letter of each computational name indicates the grid point considered and
 $r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants).
 Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points,
-thus no specific arrays are defined at $w$-points. 
+thus no specific arrays are defined at $w$-points.
 
 Note that the definition of the scale factors
@@ -405,12 +410,13 @@
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[width=0.90\textwidth]{Fig_zgr_e3}
+    \includegraphics[]{Fig_zgr_e3}
     \caption{
       \protect\label{fig:zgr_e3}
       Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,
       and (b) analytically derived grid-point position and scale factors.
-      For both grids here,
-      the same $w$-point depth has been chosen but in (a) the $t$-points are set half way between $w$-points while
-      in (b) they are defined from an analytical function: $z(k)=5\,(k-1/2)^3 - 45\,(k-1/2)^2 + 140\,(k-1/2) - 150$.
+      For both grids here, the same $w$-point depth has been chosen but
+      in (a) the $t$-points are set half way between $w$-points while
+      in (b) they are defined from an analytical function:
+      $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$.
       Note the resulting difference between the value of the grid-size $\Delta_k$ and
       those of the scale factor $e_k$.
@@ -426,5 +432,4 @@
 \label{subsec:DOM_hgr_msh_choice}
 
-
 % -------------------------------------------------------------------------------------------------------------
 %        Grid files
@@ -434,5 +439,5 @@
 
 All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks)
-can be saved in files if \np{nn\_msh} $\not= 0$ (namelist variable in \ngn{namdom}).
+can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}).
 This can be particularly useful for plots and off-line diagnostics.
 In some cases, the user may choose to make a local modification of a scale factor in the code.
@@ -441,5 +446,5 @@
 An example is Gibraltar Strait in the ORCA2 configuration.
 When such modifications are done,
-the output grid written when \np{nn\_msh} $\not= 0$ is no more equal to the input grid.
+the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid.
 
 % ================================================================
@@ -466,5 +471,5 @@
 \begin{figure}[!tb]
   \begin{center}
-    \includegraphics[width=1.0\textwidth]{Fig_z_zps_s_sps}
+    \includegraphics[]{Fig_z_zps_s_sps}
     \caption{
       \protect\label{fig:z_zps_s_sps}
@@ -475,5 +480,5 @@
       (d) hybrid $s-z$ coordinate,
       (e) hybrid $s-z$ coordinate with partial step, and
-      (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}\forcode{ = .false.}).
+      (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}~\forcode{= .false.}).
       Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).
     }
@@ -485,15 +490,15 @@
 must be done once of all at the beginning of an experiment.
 It is not intended as an option which can be enabled or disabled in the middle of an experiment.
-Three main choices are offered (\autoref{fig:z_zps_s_sps}a to c):
-$z$-coordinate with full step bathymetry (\np{ln\_zco}\forcode{ = .true.}),
-$z$-coordinate with partial step bathymetry (\np{ln\_zps}\forcode{ = .true.}),
-or generalized, $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}).
+Three main choices are offered (\autoref{fig:z_zps_s_sps}):
+$z$-coordinate with full step bathymetry (\np{ln\_zco}~\forcode{= .true.}),
+$z$-coordinate with partial step bathymetry (\np{ln\_zps}~\forcode{= .true.}),
+or generalized, $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}).
 Hybridation of the three main coordinates are available:
-$s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}e).
+$s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}).
 By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that
-the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}f).
-When a linear free surface is assumed (\np{ln\_linssh}\forcode{ = .true.}),
-the vertical coordinate are fixed in time, but the seawater can move up and down across the z=0 surface
-(in other words, the top of the ocean in not a rigid-lid). 
+the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}).
+When a linear free surface is assumed (\np{ln\_linssh}~\forcode{= .true.}),
+the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface
+(in other words, the top of the ocean in not a rigid-lid).
 The last choice in terms of vertical coordinate concerns the presence (or not) in
 the model domain of ocean cavities beneath ice shelves.
@@ -502,27 +507,25 @@
 and partial step are also applied at the ocean/ice shelf interface.
 
-Contrary to the horizontal grid, the vertical grid is computed in the code and
-no provision is made for reading it from a file.
+Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for
+reading it from a file.
 The only input file is the bathymetry (in meters) (\ifile{bathy\_meter})
 \footnote{
   N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file,
-  so that the computation of the number of wet ocean point in each water column is by-passed
-}. 
-If \np{ln\_isfcav}\forcode{ = .true.},
-an extra file input file describing the ice shelf draft (in meters) (\ifile{isf\_draft\_meter}) is needed.
+  so that the computation of the number of wet ocean point in each water column is by-passed}.
+If \np{ln\_isfcav}~\forcode{= .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing
+the ice shelf draft (in meters) is needed.
 
 After reading the bathymetry, the algorithm for vertical grid definition differs between the different options:
 \begin{description}
 \item[\textit{zco}]
-  set a reference coordinate transformation $z_0 (k)$, and set $z(i,j,k,t)=z_0 (k)$.
+  set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$.
 \item[\textit{zps}]
-  set a reference coordinate transformation $z_0 (k)$,
-  and calculate the thickness of the deepest level at each $(i,j)$ point using the bathymetry,
-  to obtain the final three-dimensional depth and scale factor arrays.
+  set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at
+  each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays.
 \item[\textit{sco}]
-  smooth the bathymetry to fulfil the hydrostatic consistency criteria and
+  smooth the bathymetry to fulfill the hydrostatic consistency criteria and
   set the three-dimensional transformation.
 \item[\textit{s-z} and \textit{s-zps}]
-  smooth the bathymetry to fulfil the hydrostatic consistency criteria and
+  smooth the bathymetry to fulfill the hydrostatic consistency criteria and
   set the three-dimensional transformation $z(i,j,k)$,
   and possibly introduce masking of extra land points to better fit the original bathymetry file.
@@ -532,13 +535,12 @@
 %%%
 
-Unless a linear free surface is used (\np{ln\_linssh}\forcode{ = .false.}),
+Unless a linear free surface is used (\np{ln\_linssh}~\forcode{= .false.}),
 the arrays describing the grid point depths and vertical scale factors are three set of
 three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step.
-The time at which they are defined is indicated by a suffix:$\_b$, $\_n$, or $\_a$, respectively.
+The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively.
 They are updated at each model time step using a fixed reference coordinate system which
 computer names have a $\_0$ suffix.
-When the linear free surface option is used (\np{ln\_linssh}\forcode{ = .true.}),
-\textit{before}, \textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart. 
-
+When the linear free surface option is used (\np{ln\_linssh}~\forcode{= .true.}), \textit{before},
+\textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -551,12 +553,12 @@
 (found in \ngn{namdom} namelist): 
 \begin{description}
-\item[\np{nn\_bathy}\forcode{ = 0}]:
+\item[\np{nn\_bathy}~\forcode{= 0}]:
   a flat-bottom domain is defined.
   The total depth $z_w (jpk)$ is given by the coordinate transformation.
-  The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}. 
-\item[\np{nn\_bathy}\forcode{ = -1}]:
+  The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}.
+\item[\np{nn\_bathy}~\forcode{= -1}]:
   a domain with a bump of topography one third of the domain width at the central latitude.
-  This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. 
-\item[\np{nn\_bathy}\forcode{ = 1}]:
+  This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount.
+\item[\np{nn\_bathy}~\forcode{= 1}]:
   read a bathymetry and ice shelf draft (if needed).
   The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at
@@ -569,21 +571,21 @@
   The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at
   each grid point of the model grid.
-  This file is only needed if \np{ln\_isfcav}\forcode{ = .true.}.
+  This file is only needed if \np{ln\_isfcav}~\forcode{= .true.}.
   Defining the ice shelf draft will also define the ice shelf edge and the grounding line position.
 \end{description}
 
 When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies
-(e.g, great lakes, Caspian sea...)
-even if the model resolution does not allow their communication with the rest of the ocean.
+(\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with
+the rest of the ocean.
 This is unnecessary when the ocean is forced by fixed atmospheric conditions,
 so these seas can be removed from the ocean domain.
 The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}),
-but the code has to be adapted to the user's configuration. 
+but the code has to be adapted to the user's configuration.
 
 % -------------------------------------------------------------------------------------------------------------
 %        z-coordinate  and reference coordinate transformation
 % -------------------------------------------------------------------------------------------------------------
-\subsection[$Z$-coordinate (\protect\np{ln\_zco}\forcode{ = .true.}) and ref. coordinate]
-				{$Z$-coordinate (\protect\np{ln\_zco}\forcode{ = .true.}) and reference coordinate}
+\subsection[$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and ref. coordinate]
+				{$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and reference coordinate}
 \label{subsec:DOM_zco}
 
@@ -591,5 +593,5 @@
 \begin{figure}[!tb]
   \begin{center}
-    \includegraphics[width=0.90\textwidth]{Fig_zgr}
+    \includegraphics[]{Fig_zgr}
     \caption{
       \protect\label{fig:zgr}
@@ -602,9 +604,10 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ and $gdepw_0$ for
-$t$- and $w$-points, respectively.
-As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the ocean surface.
+The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points,
+respectively.
+As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels.
+$gdepw_0(1)$ is the ocean surface.
 There are at most \jp{jpk}-1 $t$-points inside the ocean,
-the additional $t$-point at $jk=jpk$ is below the sea floor and is not used.
+the additional $t$-point at $jk = jpk$ is below the sea floor and is not used.
 The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose
 analytical derivative with respect to $k$ provides the vertical scale factors.
@@ -613,67 +616,62 @@
 using parameters provided in the \ngn{namcfg} namelist.
 
-It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr=0}).
-In that case,
-the parameters \jp{jpk} (number of $w$-levels) and \np{pphmax} (total ocean depth in meters) fully define the grid. 
+It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}~\forcode{= 0}).
+In that case, the parameters \jp{jpk} (number of $w$-levels) and
+\np{pphmax} (total ocean depth in meters) fully define the grid.
 
 For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface.
 The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps): 
-\begin{equation}
+\begin{gather}
   \label{eq:DOM_zgr_ana_1}
-  \begin{split}
-    z_0 (k) 	&= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\
-    e_3^0 (k) 	&= \left| -h_0 -h_1 \;\tanh \left( {{(k-h_{th} )} / {h_{cr} }} \right) \right|
-  \end{split}
-\end{equation}
-where $k=1$ to \jp{jpk} for $w$-levels and $k=1$ to $k=1$ for $T-$levels.
+    z_0  (k) = h_{sur} - h_0 \; k - \; h_1 \; \log  \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\
+    e_3^0(k) = \lt|    - h_0      -    h_1 \; \tanh \big[        (k - h_{th}) / h_{cr}  \big] \rt|
+\end{gather}
+where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels.
 Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with
 a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}).
 
-If the ice shelf cavities are opened (\np{ln\_isfcav}\forcode{ = .true.}), the definition of $z_0$ is the same.
+If the ice shelf cavities are opened (\np{ln\_isfcav}~\forcode{= .true.}), the definition of $z_0$ is the same.
 However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to:
 \begin{equation}
   \label{eq:DOM_zgr_ana_2}
   \begin{split}
-    e_3^T(k) &= z_W (k+1) - z_W (k)   \\
-    e_3^W(k) &= z_T (k)   - z_T (k-1) \\
+    e_3^T(k) &= z_W (k + 1) - z_W (k    ) \\
+    e_3^W(k) &= z_T (k    ) - z_T (k - 1)
   \end{split}
 \end{equation}
 This formulation decrease the self-generated circulation into the ice shelf cavity 
 (which can, in extreme case, leads to blow up).\\
-
  
-The most used vertical grid for ORCA2 has $10~m$ ($500~m)$ resolution in the surface (bottom) layers and
+The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and
 a depth which varies from 0 at the sea surface to a minimum of $-5000~m$.
 This leads to the following conditions:
 \begin{equation}
   \label{eq:DOM_zgr_coef}
-  \begin{split}
-    e_3 (1+1/2)		&=10. \\
-    e_3 (jpk-1/2)	&=500. \\
-    z(1)			&=0. \\
-    z(jpk)			&=-5000. \\
-  \end{split}
+  \begin{array}{ll}
+    e_3 (1   + 1/2) =  10. & z(1  ) =     0. \\
+    e_3 (jpk - 1/2) = 500. & z(jpk) = -5000.
+  \end{array}
 \end{equation}
 
-With the choice of the stretching $h_{cr} =3$ and the number of levels \jp{jpk}=$31$,
-the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in
+With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$,
+the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in
 \autoref{eq:DOM_zgr_ana_2} have been determined such that
 \autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method.
 For the first standard ORCA2 vertical grid this led to the following values:
-$h_{sur} =4762.96$, $h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$.
+$h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$.
 The resulting depths and scale factors as a function of the model levels are shown in
 \autoref{fig:zgr} and given in \autoref{tab:orca_zgr}.
-Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. 
-
-Rather than entering parameters $h_{sur}$, $h_{0}$, and $h_{1}$ directly, it is possible to recalculate them.
-In that case the user sets \np{ppsur}\forcode{ = }\np{ppa0}\forcode{ = }\np{ppa1}\forcode{ = 999999}.,
+Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist.
+
+Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them.
+In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$.,
 in \ngn{namcfg} namelist, and specifies instead the four following parameters:
 \begin{itemize}
 \item
-  \np{ppacr}=$h_{cr} $: stretching factor (nondimensional).
+  \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional).
   The larger \np{ppacr}, the smaller the stretching.
   Values from $3$ to $10$ are usual.
 \item
-  \np{ppkth}=$h_{th} $: is approximately the model level at which maximum stretching occurs
+  \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs
   (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk})
 \item
@@ -683,6 +681,5 @@
 \end{itemize}
 As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are:
-\jp{jpk}\forcode{ = 46}, \np{ppacr}\forcode{ = 9}, \np{ppkth}\forcode{ = 23.563},
-\np{ppdzmin}\forcode{ = 6}m, \np{pphmax}\forcode{ = 5750}m.
+\jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$.
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -691,36 +688,68 @@
     \begin{tabular}{c||r|r|r|r}
       \hline
-      \textbf{LEVEL}& \textbf{gdept\_1d}& \textbf{gdepw\_1d}& \textbf{e3t\_1d }& \textbf{e3w\_1d  } \\ \hline
-      1	&	\textbf{  5.00}	&       0.00 &	\textbf{ 10.00} &	 10.00 \\	\hline
-      2	&	\textbf{15.00}	& 	  10.00 &	\textbf{ 10.00} &	 10.00 \\	\hline
-      3	&	\textbf{25.00}	&	  20.00 &	\textbf{ 10.00} & 	 10.00 \\	\hline
-      4	&	\textbf{35.01}	&	  30.00 & 	\textbf{ 10.01} & 	 10.00 \\	\hline
-      5	&	\textbf{45.01}	&	  40.01 &	\textbf{ 10.01} &	 10.01 \\	\hline
-      6	&	\textbf{55.03}	&	  50.02 &	\textbf{ 10.02} & 	 10.02 \\	\hline
-      7	&	\textbf{65.06}	&	  60.04 &	\textbf{ 10.04} &	 10.03 \\	\hline
-      8	&	\textbf{75.13}	&	  70.09 &	\textbf{ 10.09} &	 10.06 \\	\hline
-      9	&	\textbf{85.25}	&	  80.18 &	\textbf{ 10.17} &	 10.12 \\	\hline
-      10	& 	\textbf{95.49}	& 	  90.35 &	\textbf{ 10.33} &	 10.24 \\	\hline
-      11	& 	\textbf{105.97}	& 	 100.69 &	\textbf{ 10.65} &	 10.47 \\	\hline
-      12	& 	\textbf{116.90}	& 	 111.36 &	\textbf{ 11.27} &	 10.91 \\	\hline
-      13	& 	\textbf{128.70}	& 	 122.65 &	\textbf{ 12.47} &	 11.77 \\	\hline
-      14	& 	\textbf{142.20}	& 	 135.16 &	\textbf{ 14.78} &	 13.43 \\	\hline
-      15	& 	\textbf{158.96}	& 	 150.03 &	\textbf{ 19.23} &	 16.65 \\	\hline
-      16	& 	\textbf{181.96}	& 	 169.42 &	\textbf{ 27.66} &	 22.78 \\	\hline
-      17	& 	\textbf{216.65}	& 	 197.37 & 	\textbf{ 43.26} &	 34.30 \\ \hline
-      18	& 	\textbf{272.48}	& 	 241.13 & 	\textbf{ 70.88} &	 55.21 \\ \hline
-      19	& 	\textbf{364.30}	& 	 312.74 & 	\textbf{116.11} &	 90.99 \\ \hline
-      20	& 	\textbf{511.53}	& 	 429.72 & 	\textbf{181.55} & 	146.43 \\ \hline
-      21	& 	\textbf{732.20}	& 	 611.89 & 	\textbf{261.03} & 	220.35 \\ \hline
-      22	& 	\textbf{1033.22}&	 872.87 & 	\textbf{339.39} & 	301.42 \\ \hline
-      23	& 	\textbf{1405.70}&	1211.59 & \textbf{402.26} & 	373.31 \\ \hline
-      24	& 	\textbf{1830.89}&	1612.98 & \textbf{444.87} & 	426.00 \\ \hline
-      25	& 	\textbf{2289.77}&	2057.13 & \textbf{470.55} & 	459.47 \\ \hline
-      26	& 	\textbf{2768.24}&	2527.22 & \textbf{484.95} & 	478.83 \\ \hline
-      27	& 	\textbf{3257.48}&	3011.90 & \textbf{492.70} & 	489.44 \\ \hline
-      28	& 	\textbf{3752.44}&	3504.46 & \textbf{496.78} & 	495.07 \\ \hline
-      29	& 	\textbf{4250.40}&	4001.16 & \textbf{498.90} & 	498.02 \\ \hline
-      30	& 	\textbf{4749.91}&	4500.02 & \textbf{500.00} &	499.54 \\ \hline
-      31	& 	\textbf{5250.23}&	5000.00 &	\textbf{500.56} &	500.33 \\ \hline
+      \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\
+      \hline
+      1              & \textbf{     5.00} &               0.00 & \textbf{   10.00} &            10.00 \\
+      \hline
+      2              & \textbf{    15.00} &	            10.00 & \textbf{   10.00} &            10.00 \\
+      \hline
+      3              & \textbf{    25.00} &              20.00 & \textbf{   10.00} &            10.00 \\
+      \hline
+      4              & \textbf{    35.01} &	            30.00 & \textbf{   10.01} &            10.00 \\
+      \hline
+      5              & \textbf{    45.01} &              40.01 & \textbf{   10.01} &            10.01 \\
+      \hline
+      6              & \textbf{    55.03} &              50.02 & \textbf{   10.02} &            10.02 \\
+      \hline
+      7              & \textbf{    65.06} &              60.04 & \textbf{   10.04} &            10.03 \\
+      \hline
+      8              & \textbf{    75.13} &              70.09 & \textbf{   10.09} &            10.06 \\
+      \hline
+      9              & \textbf{    85.25} &              80.18 & \textbf{   10.17} &            10.12 \\
+      \hline
+      10             & \textbf{    95.49} &              90.35 & \textbf{   10.33} &            10.24 \\
+      \hline
+      11             & \textbf{   105.97} &             100.69 & \textbf{   10.65} &            10.47 \\
+      \hline
+      12             & \textbf{   116.90} &             111.36 & \textbf{   11.27} &            10.91 \\
+      \hline
+      13             & \textbf{   128.70} &             122.65 & \textbf{   12.47} &            11.77 \\
+      \hline
+      14             & \textbf{   142.20} &             135.16 & \textbf{   14.78} &            13.43 \\
+      \hline
+      15             & \textbf{   158.96} &             150.03 & \textbf{   19.23} &            16.65 \\
+      \hline
+      16             & \textbf{   181.96} &             169.42 & \textbf{   27.66} &            22.78 \\
+      \hline
+      17             & \textbf{   216.65} &             197.37 & \textbf{   43.26} &            34.30 \\
+      \hline
+      18             & \textbf{   272.48} &             241.13 & \textbf{   70.88} &            55.21 \\
+      \hline
+      19             & \textbf{   364.30} &             312.74 & \textbf{  116.11} &            90.99 \\
+      \hline
+      20             & \textbf{   511.53} &             429.72 & \textbf{  181.55} &           146.43 \\
+      \hline
+      21             & \textbf{   732.20} &             611.89 & \textbf{  261.03} &           220.35 \\
+      \hline
+      22             & \textbf{  1033.22} &             872.87 & \textbf{  339.39} &           301.42 \\
+      \hline
+      23             & \textbf{  1405.70} &            1211.59 & \textbf{  402.26} &           373.31 \\
+      \hline
+      24             & \textbf{  1830.89} &            1612.98 & \textbf{  444.87} &           426.00 \\
+      \hline
+      25             & \textbf{  2289.77} &            2057.13 & \textbf{  470.55} &           459.47 \\
+      \hline
+      26             & \textbf{  2768.24} &            2527.22 & \textbf{  484.95} &           478.83 \\
+      \hline
+      27             & \textbf{  3257.48} &            3011.90 & \textbf{  492.70} &           489.44 \\
+      \hline
+      28             & \textbf{  3752.44} &            3504.46 & \textbf{  496.78} &           495.07 \\
+      \hline
+      29             & \textbf{  4250.40} &            4001.16 & \textbf{  498.90} &           498.02 \\
+      \hline
+      30             & \textbf{  4749.91} &            4500.02 & \textbf{  500.00} &           499.54 \\
+      \hline
+      31             & \textbf{  5250.23} &            5000.00 & \textbf{  500.56} &           500.33 \\
+      \hline
     \end{tabular}
   \end{center}
@@ -736,5 +765,5 @@
 %        z-coordinate with partial step
 % -------------------------------------------------------------------------------------------------------------
-\subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}\forcode{ = .true.})}
+\subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}~\forcode{= .true.})}
 \label{subsec:DOM_zps}
 %--------------------------------------------namdom-------------------------------------------------------
@@ -744,6 +773,6 @@
 
 In $z$-coordinate partial step,
-the depths of the model levels are defined by the reference analytical function $z_0 (k)$ as described in
-the previous section, \emph{except} in the bottom layer.
+the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in
+the previous section, \textit{except} in the bottom layer.
 The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to
 allow a better representation of the bathymetry, especially in the case of small slopes
@@ -752,9 +781,9 @@
 With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$.
 The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$:
-the maximum thickness allowed is $2*e_{3t}(jpk-1)$.
+the maximum thickness allowed is $2*e_{3t}(jpk - 1)$.
 This has to be kept in mind when specifying values in \ngn{namdom} namelist,
 as the maximum depth \np{pphmax} in partial steps:
-for example, with \np{pphmax}$=5750~m$ for the DRAKKAR 45 layer grid,
-the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk-1)$ being $250~m$).
+for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid,
+the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$).
 Two variables in the namdom namelist are used to define the partial step vertical grid.
 The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is
@@ -767,5 +796,5 @@
 %        s-coordinate
 % -------------------------------------------------------------------------------------------------------------
-\subsection{$S$-coordinate (\protect\np{ln\_sco}\forcode{ = .true.})}
+\subsection{$S$-coordinate (\protect\np{ln\_sco}~\forcode{= .true.})}
 \label{subsec:DOM_sco}
 %------------------------------------------nam_zgr_sco---------------------------------------------------
@@ -774,14 +803,12 @@
 %--------------------------------------------------------------------------------------------------------------
 Options are defined in \ngn{namzgr\_sco}.
-In $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}), the depth and thickness of the model levels are defined from
+In $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}), the depth and thickness of the model levels are defined from
 the product of a depth field and either a stretching function or its derivative, respectively:
 
-\[
+\begin{align*}
   % \label{eq:DOM_sco_ana}
-  \begin{split}
-    z(k) 		&= h(i,j) \; z_0(k)	\\
-    e_3(k)	&= h(i,j) \; z_0'(k)
-  \end{split}
-\]
+  z(k)   &= h(i,j) \; z_0 (k) \\
+  e_3(k) &= h(i,j) \; z_0'(k)
+\end{align*}
 
 where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and
@@ -789,7 +816,7 @@
 The depth field $h$ is not necessary the ocean depth,
 since a mixed step-like and bottom-following representation of the topography can be used
-(\autoref{fig:z_zps_s_sps}d-e) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}f).
+(\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}).
 The namelist parameter \np{rn\_rmax} determines the slope at which
-the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 
+the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate.
 The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as
 the minimum and maximum depths at which the terrain-following vertical coordinate is calculated.
@@ -799,9 +826,9 @@
 
 The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true
-(\np{ln\_s\_SH94}\forcode{ = .false.} and \np{ln\_s\_SF12}\forcode{ = .false.}). 
+(\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}).
 This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}:
 
 \[
-  z = s_{min}+C\left(s\right)\left(H-s_{min}\right)
+  z = s_{min} + C (s) (H - s_{min})
   % \label{eq:SH94_1}
 \]
@@ -810,26 +837,20 @@
 allows a $z$-coordinate to placed on top of the stretched coordinate,
 and $z$ is the depth (negative down from the asea surface).
+\begin{gather*}
+  s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1
+  % \label{eq:DOM_s}
+ \\
+  % \label{eq:DOM_sco_function}
+  C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)}
+\end{gather*}
+
+A stretching function,
+modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}),
+is also available and is more commonly used for shelf seas modelling:
 
 \[
-  s = -\frac{k}{n-1} \quad \text{ and } \quad 0 \leq k \leq n-1
-  % \label{eq:DOM_s}
-\]
-
-\[
-  % \label{eq:DOM_sco_function}
-  \begin{split}
-    C(s)	&=  \frac{ \left[	  \tanh{ \left( \theta \, (s+b) \right)}
-        - \tanh{ \left(  \theta \, b      \right)}  \right]}
-    {2\;\sinh \left( \theta \right)}
-  \end{split}
-\]
-
-A stretching function,
-modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}\forcode{ = .true.}),
-is also available and is more commonly used for shelf seas modelling:
-
-\[
-  C\left(s\right) =   \left(1 - b \right)\frac{ \sinh\left( \theta s\right)}{\sinh\left(\theta\right)} +      \\
-  b\frac{ \tanh \left[ \theta \left(s + \frac{1}{2} \right)\right] - \tanh\left(\frac{\theta}{2}\right)}{ 2\tanh\left (\frac{\theta}{2}\right)}
+  C(s) =   (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)}
+         + b       \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] -   \tanh \lt( \frac{\theta}{2} \rt)}
+                        {                                                  2 \tanh \lt( \frac{\theta}{2} \rt)}
   % \label{eq:SH94_2}
 \]
@@ -838,5 +859,5 @@
 \begin{figure}[!ht]
   \begin{center}
-    \includegraphics[width=1.0\textwidth]{Fig_sco_function}
+    \includegraphics[]{Fig_sco_function}
     \caption{
       \protect\label{fig:sco_function}
@@ -848,8 +869,7 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-where $H_c$ is the critical depth (\np{rn\_hc}) at which
-the coordinate transitions from pure $\sigma$ to the stretched coordinate,
-and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and bottom control parameters such that
-$0\leqslant \theta \leqslant 20$, and $0\leqslant b\leqslant 1$.
+where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to
+the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and
+bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$.
 $b$ has been designed to allow surface and/or bottom increase of the vertical resolution
 (\autoref{fig:sco_function}).
@@ -859,21 +879,20 @@
 In this case the a stretching function $\gamma$ is defined such that:
 
-\[
-  z = -\gamma h \quad \text{ with } \quad 0 \leq \gamma \leq 1
+\begin{equation}
+  z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1
   % \label{eq:z}
-\]
+\end{equation}
 
 The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate:
 
-\[
+\begin{gather*}
   % \label{eq:DOM_gamma_deriv}
-  \gamma= A\left(\sigma-\frac{1}{2}\left(\sigma^{2}+f\left(\sigma\right)\right)\right)+B\left(\sigma^{3}-f\left(\sigma\right)\right)+f\left(\sigma\right)
-\]
-
-Where:
-\[
+  \gamma =   A \lt( \sigma   - \frac{1}{2} (\sigma^2     + f (\sigma)) \rt)
+           + B \lt( \sigma^3 - f           (\sigma) \rt) + f (\sigma)       \\
+  \intertext{Where:}
   % \label{eq:DOM_gamma}
-  f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2} \quad \text{ and } \quad \sigma = \frac{k}{n-1}
-\]
+  f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2}
+  \quad \text{and} \quad \sigma = \frac{k}{n - 1}
+\end{gather*}
 
 This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of
@@ -892,15 +911,15 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht]
-   \includegraphics[width=1.0\textwidth]{Fig_DOM_compare_coordinates_surface}
-   \caption{
-     A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines),
-     a 50 level $Z$-coordinate (contoured surfaces) and
-     the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in
-     the surface 100m for a idealised bathymetry that goes from 50m to 5500m depth.
-     For clarity every third coordinate surface is shown.
-   }
-   \label{fig:fig_compare_coordinates_surface}
+  \includegraphics[]{Fig_DOM_compare_coordinates_surface}
+  \caption{
+    A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines),
+    a 50 level $Z$-coordinate (contoured surfaces) and
+    the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface $100~m$ for
+    a idealised bathymetry that goes from $50~m$ to $5500~m$ depth.
+    For clarity every third coordinate surface is shown.
+  }
+  \label{fig:fig_compare_coordinates_surface}
 \end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 This gives a smooth analytical stretching in computational space that is constrained to
@@ -925,10 +944,10 @@
 
 % -------------------------------------------------------------------------------------------------------------
-%        \zstar- or \sstar-coordinate
-% -------------------------------------------------------------------------------------------------------------
-\subsection{$Z^*$- or $S^*$-coordinate (\protect\np{ln\_linssh}\forcode{ = .false.}) }
+%        z*- or s*-coordinate
+% -------------------------------------------------------------------------------------------------------------
+\subsection{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}~\forcode{= .false.})}
 \label{subsec:DOM_zgr_star}
 
-This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 
+This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site.
 
 %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances
@@ -940,11 +959,10 @@
 \label{subsec:DOM_msk}
 
-Whatever the vertical coordinate used,
-the model offers the possibility of representing the bottom topography with steps that
-follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}.
-The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy,
-which gives the number of ocean levels (\ie those that are not masked) at each $t$-point.
-mbathy is computed from the meter bathymetry using the definiton of gdept as
-the number of $t$-points which gdept $\leq$ bathy.
+Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with
+steps that follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}.
+The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which
+gives the number of ocean levels (\ie those that are not masked) at each $t$-point.
+mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which
+gdept $\leq$ bathy.
 
 Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after
@@ -954,38 +972,42 @@
 As for the representation of bathymetry, a 2D integer array, misfdep, is created.
 misfdep defines the level of the first wet $t$-point.
-All the cells between $k=1$ and $misfdep(i,j)-1$ are masked.
-By default, misfdep(:,:)=1 and no cells are masked.
+All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked.
+By default, $misfdep(:,:) = 1$ and no cells are masked.
 
 In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into 
 the cavities are performed in the \textit{zgr\_isf} routine.
-The compatibility between ice shelf draft and bathymetry is checked. 
-All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked). 
+The compatibility between ice shelf draft and bathymetry is checked.
+All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked).
 If only one cell on the water column is opened at $t$-, $u$- or $v$-points,
 the bathymetry or the ice shelf draft is dug to fit this constrain.
-If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked.\\ 
+If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked.
 
 From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows:
-\begin{align*}
-  tmask(i,j,k) &= \begin{cases}   \; 0&   \text{ if $k < misfdep(i,j) $ } \\
-    \; 1&   \text{ if $misfdep(i,j) \leq k\leq mbathy(i,j)$  }    \\
-    \; 0&   \text{ if $k > mbathy(i,j)$  }    \end{cases}     \\
-  umask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)	\\
-  vmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j+1,k)	\\
-  fmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)	\\
-               &    \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\
-  wmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j,k-1) \text{ with } wmask(i,j,1) = tmask(i,j,1)
-\end{align*}
+\begin{alignat*}{2}
+  tmask(i,j,k) &= &  &
+    \begin{cases}
+                  0 &\text{if $                  k  <    misfdep(i,j)$} \\
+                  1 &\text{if $misfdep(i,j) \leq k \leq   mbathy(i,j)$} \\
+                  0 &\text{if $                  k  >     mbathy(i,j)$}
+    \end{cases}
+  \\
+  umask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\
+  vmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j + 1,k) \\
+  fmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\
+               &  &* &tmask(i,j,k) * tmask(i + 1,j,    k) \\
+  wmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j,k - 1) \\
+  \text{with~} wmask(i,j,1) &= & &tmask(i,j,1)
+\end{alignat*}
 
 Note that, without ice shelves cavities,
 masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}).
 Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 
-exactly in the same way as for the bottom boundary. 
+exactly in the same way as for the bottom boundary.
 
 The specification of closed lateral boundaries requires that at least
 the first and last rows and columns of the \textit{mbathy} array are set to zero.
-In the particular case of an east-west cyclical boundary condition,
-\textit{mbathy} has its last column equal to the second one and its first column equal to the last but one 
-(and so too the mask arrays) (see \autoref{fig:LBC_jperio}).
-
+In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to
+the second one and its first column equal to the last but one (and so too the mask arrays)
+(see \autoref{fig:LBC_jperio}).
 
 % ================================================================
@@ -1000,7 +1022,8 @@
 
 Options are defined in \ngn{namtsd}.
-By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter.
+By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and
+salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter.
 \begin{description}
-\item[\np{ln\_tsd\_init}\forcode{ = .true.}]
+\item[\np{ln\_tsd\_init}~\forcode{= .true.}]
   use a T and S input files that can be given on the model grid itself or on their native input data grid.
   In the latter case,
@@ -1009,7 +1032,7 @@
   The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures.
   The computation is done in the \mdl{dtatsd} module.
-\item[\np{ln\_tsd\_init}\forcode{ = .false.}]
-  use constant salinity value of 35.5 psu and an analytical profile of temperature (typical of the tropical ocean),
-  see \rou{istate\_t\_s} subroutine called from \mdl{istate} module.
+\item[\np{ln\_tsd\_init}~\forcode{= .false.}]
+  use constant salinity value of $35.5~psu$ and an analytical profile of temperature
+  (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module.
 \end{description}
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 10501)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 10502)
@@ -15,19 +15,15 @@
 %STEVEN :  is the use of the word "positive" to describe a scheme enough, or should it be "positive definite"? I added a comment to this effect on some instances of this below
 
-%\newpage
-
-Using the representation described in \autoref{chap:DOM},
-several semi-discrete space forms of the tracer equations are available depending on
-the vertical coordinate used and on the physics used.
+Using the representation described in \autoref{chap:DOM}, several semi -discrete space forms of
+the tracer equations are available depending on the vertical coordinate used and on the physics used.
 In all the equations presented here, the masking has been omitted for simplicity.
-One must be aware that all the quantities are masked fields and
-that each time a mean or difference operator is used,
-the resulting field is multiplied by a mask.
+One must be aware that all the quantities are masked fields and that each time a mean or
+difference operator is used, the resulting field is multiplied by a mask.
 
 The two active tracers are potential temperature and salinity.
 Their prognostic equations can be summarized as follows:
 \[
-  \text{NXT} = \text{ADV}+\text{LDF}+\text{ZDF}+\text{SBC}
-  \ (+\text{QSR})\ (+\text{BBC})\ (+\text{BBL})\ (+\text{DMP})
+  \text{NXT} =     \text{ADV} + \text{LDF} + \text{ZDF} + \text{SBC}
+               + \{\text{QSR},  \text{BBC},  \text{BBL},  \text{DMP}\}
 \]
 
@@ -39,8 +35,10 @@
 The terms QSR, BBC, BBL and DMP are optional.
 The external forcings and parameterisations require complex inputs and complex calculations
-(\eg bulk formulae, estimation of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and 
-described in \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively.
-Note that \mdl{tranpc}, the non-penetrative convection module, although located in the NEMO/OPA/TRA directory as
-it directly modifies the tracer fields, is described with the model vertical physics (ZDF) together with
+(\eg bulk formulae, estimation of mixing coefficients) that are carried out in the SBC,
+LDF and ZDF modules and described in \autoref{chap:SBC}, \autoref{chap:LDF} and
+\autoref{chap:ZDF}, respectively.
+Note that \mdl{tranpc}, the non-penetrative convection module, although located in
+the \path{./src/OCE/TRA} directory as it directly modifies the tracer fields,
+is described with the model vertical physics (ZDF) together with
 other available parameterization of convection.
 
@@ -50,12 +48,12 @@
 
 The different options available to the user are managed by namelist logicals or CPP keys.
-For each equation term  \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx},
+For each equation term \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx},
 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme.
 The CPP key (when it exists) is \key{traTTT}.
 The equivalent code can be found in the \textit{traTTT} or \textit{traTTT\_xxx} module,
-in the NEMO/OPA/TRA directory.
+in the \path{./src/OCE/TRA} directory.
 
 The user has the option of extracting each tendency term on the RHS of the tracer equation for output
-(\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}\forcode{ = .true.}), as described in \autoref{chap:DIA}.
+(\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~\forcode{= .true.}), as described in \autoref{chap:DIA}.
 
 % ================================================================
@@ -75,15 +73,13 @@
 \begin{equation}
   \label{eq:tra_adv}
-  ADV_\tau =-\frac{1}{b_t} \left(
-    \;\delta_i \left[ e_{2u}\,e_{3u} \;  u\; \tau_u  \right]
-    +\delta_j \left[ e_{1v}\,e_{3v}  \;  v\; \tau_v  \right] \; \right)
-  -\frac{1}{e_{3t}} \;\delta_k \left[ w\; \tau_w \right]
+  ADV_\tau = - \frac{1}{b_t} \Big(   \delta_i [ e_{2u} \, e_{3u} \; u \; \tau_u]
+                                   + \delta_j [ e_{1v} \, e_{3v} \; v \; \tau_v] \Big)
+             - \frac{1}{e_{3t}} \delta_k [w \; \tau_w]
 \end{equation}
-where $\tau$ is either T or S, and $b_t= e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells.
+where $\tau$ is either T or S, and $b_t = e_{1t} \, e_{2t} \, e_{3t}$ is the volume of $T$-cells.
 The flux form in \autoref{eq:tra_adv} implicitly requires the use of the continuity equation.
-Indeed, it is obtained by using the following equality:
-$\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ which
-results from the use of the continuity equation,  $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$
-(which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, \ie \np{ln\_linssh}\forcode{ = .true.}).
+Indeed, it is obtained by using the following equality: $\nabla \cdot (\vect U \, T) = \vect U \cdot \nabla T$ which
+results from the use of the continuity equation, $\partial_t e_3 + e_3 \; \nabla \cdot \vect U = 0$
+(which reduces to $\nabla \cdot \vect U = 0$ in linear free surface, \ie \np{ln\_linssh}~\forcode{= .true.}).
 Therefore it is of paramount importance to design the discrete analogue of the advection tendency so that
 it is consistent with the continuity equation in order to enforce the conservation properties of
@@ -94,5 +90,5 @@
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[width=0.9\textwidth]{Fig_adv_scheme}
+    \includegraphics[]{Fig_adv_scheme}
     \caption{
       \protect\label{fig:adv_scheme}
@@ -120,42 +116,42 @@
 since the normal velocity is zero there.
 At the sea surface the boundary condition depends on the type of sea surface chosen:
+
 \begin{description}
 \item[linear free surface:]
-  (\np{ln\_linssh}\forcode{ = .true.})
+  (\np{ln\_linssh}~\forcode{= .true.})
   the first level thickness is constant in time:
-  the vertical boundary condition is applied at the fixed surface $z=0$ rather than on the moving surface $z=\eta$.
+  the vertical boundary condition is applied at the fixed surface $z = 0$ rather than on
+  the moving surface $z = \eta$.
   There is a non-zero advective flux which is set for all advection schemes as
-  $\left. {\tau_w } \right|_{k=1/2} =T_{k=1} $,
-  \ie the product of surface velocity (at $z=0$) by the first level tracer value.
+  $\tau_w|_{k = 1/2} = T_{k = 1}$, \ie the product of surface velocity (at $z = 0$) by
+  the first level tracer value.
 \item[non-linear free surface:]
-  (\np{ln\_linssh}\forcode{ = .false.})
+  (\np{ln\_linssh}~\forcode{= .false.})
   convergence/divergence in the first ocean level moves the free surface up/down.
   There is no tracer advection through it so that the advective fluxes through the surface are also zero.
 \end{description}
+
 In all cases, this boundary condition retains local conservation of tracer.
 Global conservation is obtained in non-linear free surface case, but \textit{not} in the linear free surface case.
 Nevertheless, in the latter case, it is achieved to a good approximation since
 the non-conservative term is the product of the time derivative of the tracer and the free surface height,
-two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}.
-
-The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco})
-is the centred (\textit{now}) \textit{effective} ocean velocity,
-\ie the \textit{eulerian} velocity (see \autoref{chap:DYN}) plus
-the eddy induced velocity (\textit{eiv}) and/or
-the mixed layer eddy induced velocity (\textit{eiv}) when
-those parameterisations are used (see \autoref{chap:LDF}).
+two quantities that are not correlated \citep{Roullet_Madec_JGR00, Griffies_al_MWR01, Campin2004}.
+
+The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco}) is
+the centred (\textit{now}) \textit{effective} ocean velocity, \ie the \textit{eulerian} velocity
+(see \autoref{chap:DYN}) plus the eddy induced velocity (\textit{eiv}) and/or
+the mixed layer eddy induced velocity (\textit{eiv}) when those parameterisations are used
+(see \autoref{chap:LDF}).
 
 Several tracer advection scheme are proposed, namely a $2^{nd}$ or $4^{th}$ order centred schemes (CEN),
-a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT),
-a Monotone Upstream Scheme for Conservative Laws scheme (MUSCL),
-a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3),
-and a Quadratic Upstream Interpolation for Convective Kinematics with
-Estimated Streaming Terms scheme (QUICKEST).
-The choice is made in the \textit{\ngn{namtra\_adv}} namelist,
-by setting to \forcode{.true.} one of the logicals \textit{ln\_traadv\_xxx}.
-The corresponding code can be found in the \mdl{traadv\_xxx} module,
-where \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme.
-By default (\ie in the reference namelist, \ngn{namelist\_ref}), all the logicals are set to \forcode{.false.}.
-If the user does not select an advection scheme in the configuration namelist (\ngn{namelist\_cfg}),
+a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT), a Monotone Upstream Scheme for
+Conservative Laws scheme (MUSCL), a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3),
+and a Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms scheme (QUICKEST).
+The choice is made in the \ngn{namtra\_adv} namelist, by setting to \forcode{.true.} one of
+the logicals \textit{ln\_traadv\_xxx}.
+The corresponding code can be found in the \textit{traadv\_xxx.F90} module, where
+\textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme.
+By default (\ie in the reference namelist, \textit{namelist\_ref}), all the logicals are set to \forcode{.false.}.
+If the user does not select an advection scheme in the configuration namelist (\textit{namelist\_cfg}),
 the tracers will \textit{not} be advected!
 
@@ -163,28 +159,34 @@
 The choosing an advection scheme is a complex matter which depends on the model physics, model resolution,
 type of tracer, as well as the issue of numerical cost. In particular, we note that
-(1) CEN and FCT schemes require an explicit diffusion operator while the other schemes are diffusive enough so that
-they do not necessarily need additional diffusion;
-(2) CEN and UBS are not \textit{positive} schemes
-\footnote{negative values can appear in an initially strictly positive tracer field which is advected},
-implying that false extrema are permitted.
-Their use is not recommended on passive tracers;
-(3) It is recommended that the same advection-diffusion scheme is used on both active and passive tracers.
-Indeed, if a source or sink of a passive tracer depends on an active one,
-the difference of treatment of active and passive tracers can create very nice-looking frontal structures that
-are pure numerical artefacts.
+
+\begin{enumerate}
+\item
+  CEN and FCT schemes require an explicit diffusion operator while the other schemes are diffusive enough so that
+  they do not necessarily need additional diffusion;
+\item
+  CEN and UBS are not \textit{positive} schemes
+  \footnote{negative values can appear in an initially strictly positive tracer field which is advected},
+  implying that false extrema are permitted.
+  Their use is not recommended on passive tracers;
+\item
+  It is recommended that the same advection-diffusion scheme is used on both active and passive tracers.
+\end{enumerate}
+
+Indeed, if a source or sink of a passive tracer depends on an active one, the difference of treatment of active and
+passive tracers can create very nice-looking frontal structures that are pure numerical artefacts.
 Nevertheless, most of our users set a different treatment on passive and active tracers,
 that's the reason why this possibility is offered.
-We strongly suggest them to perform a sensitivity experiment using a same treatment to
-assess the robustness of their results.
+We strongly suggest them to perform a sensitivity experiment using a same treatment to assess the robustness of
+their results.
 
 % -------------------------------------------------------------------------------------------------------------
 %        2nd and 4th order centred schemes
 % -------------------------------------------------------------------------------------------------------------
-\subsection{CEN: Centred scheme (\protect\np{ln\_traadv\_cen}\forcode{ = .true.})}
+\subsection{CEN: Centred scheme (\protect\np{ln\_traadv\_cen}~\forcode{= .true.})}
 \label{subsec:TRA_adv_cen}
 
 %        2nd order centred scheme  
 
-The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}\forcode{ = .true.}.
+The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}~\forcode{= .true.}.
 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by
 setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$.
@@ -196,13 +198,13 @@
 \begin{equation}
   \label{eq:tra_adv_cen2}
-  \tau_u^{cen2} =\overline T ^{i+1/2}
+  \tau_u^{cen2} = \overline T ^{i + 1/2}
 \end{equation}
 
-CEN2 is non diffusive (\ie it conserves the tracer variance, $\tau^2)$ but dispersive
+CEN2 is non diffusive (\ie it conserves the tracer variance, $\tau^2$) but dispersive
 (\ie it may create false extrema).
 It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to
 produce a sensible solution.
 The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter,
-so $T$ in (\autoref{eq:tra_adv_cen2}) is the \textit{now} tracer value. 
+so $T$ in (\autoref{eq:tra_adv_cen2}) is the \textit{now} tracer value.
 
 Note that using the CEN2, the overall tracer advection is of second order accuracy since
@@ -216,11 +218,10 @@
 \begin{equation}
   \label{eq:tra_adv_cen4}
-  \tau_u^{cen4} =\overline{   T - \frac{1}{6}\,\delta_i \left[ \delta_{i+1/2}[T] \,\right]   }^{\,i+1/2}
+  \tau_u^{cen4} = \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \, \Big]}^{\,i + 1/2}
 \end{equation}
-In the vertical direction (\np{nn\_cen\_v}\forcode{ = 4}),
+In the vertical direction (\np{nn\_cen\_v}~\forcode{= 4}),
 a $4^{th}$ COMPACT interpolation has been prefered \citep{Demange_PhD2014}.
 In the COMPACT scheme, both the field and its derivative are interpolated, which leads, after a matrix inversion,
-spectral characteristics similar to schemes of higher order \citep{Lele_JCP1992}.
- 
+spectral characteristics similar to schemes of higher order \citep{Lele_JCP1992}. 
 
 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme but
@@ -249,8 +250,8 @@
 %        FCT scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}\forcode{ = .true.})}
+\subsection{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}~\forcode{= .true.})}
 \label{subsec:TRA_adv_tvd}
 
-The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}\forcode{ = .true.}.
+The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}~\forcode{= .true.}.
 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by
 setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$.
@@ -263,11 +264,11 @@
   \label{eq:tra_adv_fct}
   \begin{split}
-    \tau_u^{ups}&=
+    \tau_u^{ups} &=
     \begin{cases}
-      T_{i+1} 	& \text{if $\ u_{i+1/2} <     0$} \hfill \\
-      T_i   		& \text{if $\ u_{i+1/2} \geq 0$} \hfill \\
+                     T_{i + 1} & \text{if~} u_{i + 1/2} <    0 \\
+                     T_i       & \text{if~} u_{i + 1/2} \geq 0 \\
     \end{cases}
-    \\ \\
-    \tau_u^{fct}&=\tau_u^{ups} +c_u \;\left( {\tau_u^{cen} -\tau_u^{ups} } \right)
+    \\
+    \tau_u^{fct} &= \tau_u^{ups} + c_u \, \big( \tau_u^{cen} - \tau_u^{ups} \big)
   \end{split}
 \end{equation}
@@ -278,5 +279,5 @@
 The one chosen in \NEMO is described in \citet{Zalesak_JCP79}.
 $c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field.
-The resulting scheme is quite expensive but \emph{positive}.
+The resulting scheme is quite expensive but \textit{positive}.
 It can be used on both active and passive tracers.
 A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}.
@@ -294,13 +295,13 @@
 $\tau_u^{ups}$ is evaluated using the \textit{before} tracer.
 In other words, the advective part of the scheme is time stepped with a leap-frog scheme
-while a forward scheme is used for the diffusive part. 
+while a forward scheme is used for the diffusive part.
 
 % -------------------------------------------------------------------------------------------------------------
 %        MUSCL scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}\forcode{ = .true.})}
+\subsection{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}~\forcode{= .true.})}
 \label{subsec:TRA_adv_mus}
 
-The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}\forcode{ = .true.}.
+The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}~\forcode{= .true.}.
 MUSCL implementation can be found in the \mdl{traadv\_mus} module.
 
@@ -309,20 +310,20 @@
 two $T$-points (\autoref{fig:adv_scheme}).
 For example, in the $i$-direction :
-\[
+\begin{equation}
   % \label{eq:tra_adv_mus}
-  \tau_u^{mus} = \left\{
-    \begin{aligned}
-      &\tau_i  &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right)
-      &\ \widetilde{\partial _i \tau}  & \quad \text{if }\;u_{i+1/2} \geqslant 0      \\
-      &\tau_{i+1/2} &+\frac{1}{2}\;\left( 1+\frac{u_{i+1/2} \;\rdt}{e_{1u} } \right)
-      &\ \widetilde{\partial_{i+1/2} \tau } & \text{if }\;u_{i+1/2} <0
-    \end{aligned}
-  \right.
-\]
-where $\widetilde{\partial _i \tau}$ is the slope of the tracer on which a limitation is imposed to
+  \tau_u^{mus} = \lt\{
+  \begin{split}
+                       \tau_i         &+ \frac{1}{2} \lt( 1 - \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt)
+                       \widetilde{\partial_i         \tau} & \text{if~} u_{i + 1/2} \geqslant 0 \\
+                       \tau_{i + 1/2} &+ \frac{1}{2} \lt( 1 + \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt)
+                       \widetilde{\partial_{i + 1/2} \tau} & \text{if~} u_{i + 1/2} <         0
+  \end{split}
+                                                                                                      \rt.
+\end{equation}
+where $\widetilde{\partial_i \tau}$ is the slope of the tracer on which a limitation is imposed to
 ensure the \textit{positive} character of the scheme.
 
-The time stepping is performed using a forward scheme,
-that is the \textit{before} tracer field is used to evaluate $\tau_u^{mus}$.
+The time stepping is performed using a forward scheme, that is the \textit{before} tracer field is used to
+evaluate $\tau_u^{mus}$.
 
 For an ocean grid point adjacent to land and where the ocean velocity is directed toward land,
@@ -330,13 +331,13 @@
 This choice ensure the \textit{positive} character of the scheme.
 In addition, fluxes round a grid-point where a runoff is applied can optionally be computed using upstream fluxes
-(\np{ln\_mus\_ups}\forcode{ = .true.}).
+(\np{ln\_mus\_ups}~\forcode{= .true.}).
 
 % -------------------------------------------------------------------------------------------------------------
 %        UBS scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}\forcode{ = .true.})}
+\subsection{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}~\forcode{= .true.})}
 \label{subsec:TRA_adv_ubs}
 
-The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}\forcode{ = .true.}.
+The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}~\forcode{= .true.}.
 UBS implementation can be found in the \mdl{traadv\_mus} module.
 
@@ -347,18 +348,18 @@
 \begin{equation}
   \label{eq:tra_adv_ubs}
-  \tau_u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{
-    \begin{aligned}
-      &\tau"_i        	& \quad \text{if }\ u_{i+1/2} \geqslant 0      \\
-      &\tau"_{i+1}	& \quad \text{if }\ u_{i+1/2}       <       0
-    \end{aligned}
-  \right.
+  \tau_u^{ubs} = \overline T ^{i + 1/2} - \frac{1}{6}
+    \begin{cases}
+                                                      \tau"_i       & \text{if~} u_{i + 1/2} \geqslant 0 \\
+                                                      \tau"_{i + 1} & \text{if~} u_{i + 1/2} <         0
+    \end{cases}
+  \quad
+  \text{where~} \tau"_i = \delta_i \lt[ \delta_{i + 1/2} [\tau] \rt]
 \end{equation}
-where $\tau "_i =\delta_i \left[ {\delta_{i+1/2} \left[ \tau \right]} \right]$.
-
-This results in a dissipatively dominant (\ie hyper-diffusive) truncation error
+
+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.
-Nevertheless the scheme is not \emph{positive}, meaning that false extrema are permitted,
+Nevertheless the scheme is not \textit{positive}, meaning that false extrema are permitted,
 but the amplitude of such are significantly reduced over the centred second or fourth order method.
 Therefore it is not recommended that it should be applied to a passive tracer that requires positivity.
@@ -368,5 +369,5 @@
 \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}.
 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or a $4^th$ order COMPACT scheme
-(\np{nn\_cen\_v}\forcode{ = 2 or 4}).
+(\np{nn\_cen\_v}~\forcode{= 2 or 4}).
 
 For stability reasons (see \autoref{chap:STP}), the first term  in \autoref{eq:tra_adv_ubs}
@@ -382,20 +383,16 @@
 
 Note that it is straightforward to rewrite \autoref{eq:tra_adv_ubs} as follows:
-\[
+\begin{gather}
   \label{eq:traadv_ubs2}
-  \tau_u^{ubs} = \tau_u^{cen4} + \frac{1}{12} \left\{
-    \begin{aligned}
-      & + \tau"_i			& \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
-      &  - \tau"_{i+1}		& \quad \text{if }\ u_{i+1/2}       <       0
-    \end{aligned}
-  \right.
-\]
-or equivalently 
-\[
+  \tau_u^{ubs} = \tau_u^{cen4} + \frac{1}{12}
+    \begin{cases}
+      + \tau"_i       & \text{if} \ u_{i + 1/2} \geqslant 0 \\
+      - \tau"_{i + 1} & \text{if} \ u_{i + 1/2} <         0
+    \end{cases}
+  \intertext{or equivalently}
   % \label{eq:traadv_ubs2b}
-  u_{i+1/2} \ \tau_u^{ubs}
-  =u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\delta_i\left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2}
-  - \frac{1}{2} |u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
-\]
+  u_{i + 1/2} \ \tau_u^{ubs} = u_{i + 1/2} \, \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \Big]}^{\,i + 1/2}
+                             - \frac{1}{2} |u|_{i + 1/2} \, \frac{1}{6} \, \delta_{i + 1/2} [\tau"_i] \nonumber
+\end{gather}
 
 \autoref{eq:traadv_ubs2} has several advantages.
@@ -403,8 +400,7 @@
 an upstream-biased diffusion term is added.
 Secondly, this emphasises that the $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has to
-be evaluated at the \emph{now} time step using \autoref{eq:tra_adv_ubs}.
+be evaluated at the \textit{now} time step using \autoref{eq:tra_adv_ubs}.
 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy coefficient which
-is simply proportional to the velocity:
-$A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$.
+is simply proportional to the velocity: $A_u^{lm} = \frac{1}{12} \, {e_{1u}}^3 \, |u|$.
 Note the current version of NEMO uses the computationally more efficient formulation \autoref{eq:tra_adv_ubs}.
 
@@ -412,9 +408,9 @@
 %        QCK scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}\forcode{ = .true.})}
+\subsection{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}~\forcode{= .true.})}
 \label{subsec:TRA_adv_qck}
 
 The Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms (QUICKEST) scheme
-proposed by \citet{Leonard1979} is used when \np{ln\_traadv\_qck}\forcode{ = .true.}.
+proposed by \citet{Leonard1979} is used when \np{ln\_traadv\_qck}~\forcode{= .true.}.
 QUICKEST implementation can be found in the \mdl{traadv\_qck} module.
 
@@ -422,5 +418,5 @@
 \citep{Leonard1991}.
 It has been implemented in NEMO by G. Reffray (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module.
-The resulting scheme is quite expensive but \emph{positive}.
+The resulting scheme is quite expensive but \textit{positive}.
 It can be used on both active and passive tracers.
 However, the intrinsic diffusion of QCK makes its use risky in the vertical direction where
@@ -431,5 +427,4 @@
 
 %%%gmcomment   :  Cross term are missing in the current implementation....
-
 
 % ================================================================
@@ -458,5 +453,5 @@
 except for the pure vertical component that appears when a rotation tensor is used.
 This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}).
-When \np{ln\_traldf\_msc}\forcode{ = .true.}, a Method of Stabilizing Correction is used in which
+When \np{ln\_traldf\_msc}~\forcode{= .true.}, a Method of Stabilizing Correction is used in which
 the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}.
 
@@ -464,24 +459,24 @@
 %        Type of operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Type of operator (\protect\np{ln\_traldf}\{\_NONE,\_lap,\_blp\}\})]
-			     {Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) } 
+\subsection[Type of operator (\protect\np{ln\_traldf}\{\_NONE,\_lap,\_blp\}\})]{Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) } 
 \label{subsec:TRA_ldf_op}
 
 Three operator options are proposed and, one and only one of them must be selected:
+
 \begin{description}
-\item[\np{ln\_traldf\_NONE}\forcode{ = .true.}:]
+\item[\np{ln\_traldf\_NONE}~\forcode{= .true.}:]
   no operator selected, the lateral diffusive tendency will not be applied to the tracer equation.
   This option can be used when the selected advection scheme is diffusive enough (MUSCL scheme for example).
-\item[\np{ln\_traldf\_lap}\forcode{ = .true.}:]
+\item[\np{ln\_traldf\_lap}~\forcode{= .true.}:]
   a laplacian operator is selected.
-  This harmonic operator takes the following expression:  $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $,
+  This harmonic operator takes the following expression:  $\mathpzc{L}(T) = \nabla \cdot A_{ht} \; \nabla T $,
   where the gradient operates along the selected direction (see \autoref{subsec:TRA_ldf_dir}),
   and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}).
-\item[\np{ln\_traldf\_blp}\forcode{ = .true.}]:
+\item[\np{ln\_traldf\_blp}~\forcode{= .true.}]:
   a bilaplacian operator is selected.
   This biharmonic operator takes the following expression:
-  $\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$
+  $\mathpzc{B} = - \mathpzc{L}(\mathpzc{L}(T)) = - \nabla \cdot b \nabla (\nabla \cdot b \nabla T)$
   where the gradient operats along the selected direction,
-  and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$  (see \autoref{chap:LDF}).
+  and $b^2 = B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see \autoref{chap:LDF}).
   In the code, the bilaplacian operator is obtained by calling the laplacian twice.
 \end{description}
@@ -495,16 +490,14 @@
 whereas the laplacian damping time scales only like $\lambda^{-2}$.
 
-
 % -------------------------------------------------------------------------------------------------------------
 %        Direction of action
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Action direction (\protect\np{ln\_traldf}\{\_lev,\_hor,\_iso,\_triad\})]
-			     {Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) } 
+\subsection[Action direction (\protect\np{ln\_traldf}\{\_lev,\_hor,\_iso,\_triad\})]{Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) } 
 \label{subsec:TRA_ldf_dir}
 
 The choice of a direction of action determines the form of operator used.
 The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane when
-iso-level option is used (\np{ln\_traldf\_lev}\forcode{ = .true.}) or
-when a horizontal (\ie geopotential) operator is demanded in \zstar-coordinate
+iso-level option is used (\np{ln\_traldf\_lev}~\forcode{= .true.}) or
+when a horizontal (\ie geopotential) operator is demanded in \textit{z}-coordinate
 (\np{ln\_traldf\_hor} and \np{ln\_zco} equal \forcode{.true.}).
 The associated code can be found in the \mdl{traldf\_lap\_blp} module.
@@ -521,11 +514,10 @@
 
 The resulting discret form of the three operators (one iso-level and two rotated one) is given in
-the next two sub-sections. 
-
+the next two sub-sections.
 
 % -------------------------------------------------------------------------------------------------------------
 %       iso-level operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Iso-level (bi-)laplacian operator ( \protect\np{ln\_traldf\_iso}) }
+\subsection{Iso-level (bi -)laplacian operator ( \protect\np{ln\_traldf\_iso}) }
 \label{subsec:TRA_ldf_lev}
 
@@ -533,71 +525,59 @@
 \begin{equation}
   \label{eq:tra_ldf_lap}
-  D_t^{lT} =\frac{1}{b_t} \left( \;
-    \delta_{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta_{i+1/2} [T] \right]
-    + \delta_{j}\left[ A_v^{lT} \;  \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta_{j+1/2} [T] \right]  \;\right)
+  D_t^{lT} = \frac{1}{b_t} \Bigg(   \delta_{i} \lt[ A_u^{lT} \; \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [T] \rt]
+                                  + \delta_{j} \lt[ A_v^{lT} \; \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [T] \rt] \Bigg)
 \end{equation}
-where  $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells and
+where $b_t = e_{1t} \, e_{2t} \, e_{3t}$  is the volume of $T$-cells and
 where zero diffusive fluxes is assumed across solid boundaries,
 first (and third in bilaplacian case) horizontal tracer derivative are masked.
 It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module.
 The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap} in order to
-compute the iso-level bilaplacian operator. 
-
-It is a \emph{horizontal} operator (\ie acting along geopotential surfaces) in
+compute the iso-level bilaplacian operator.
+
+It is a \textit{horizontal} operator (\ie acting along geopotential surfaces) in
 the $z$-coordinate with or without partial steps, but is simply an iso-level operator in the $s$-coordinate.
-It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}\forcode{ = .true.},
-we have \np{ln\_traldf\_lev}\forcode{ = .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}\forcode{ = .true.}.
+It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}~\forcode{= .true.},
+we have \np{ln\_traldf\_lev}~\forcode{= .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~\forcode{= .true.}.
 In both cases, it significantly contributes to diapycnal mixing.
 It is therefore never recommended, even when using it in the bilaplacian case.
 
-Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}),
+Note that in the partial step $z$-coordinate (\np{ln\_zps}~\forcode{= .true.}),
 tracers in horizontally adjacent cells are located at different depths in the vicinity of the bottom.
 In this case, horizontal derivatives in (\autoref{eq:tra_ldf_lap}) at the bottom level require a specific treatment.
 They are calculated in the \mdl{zpshde} module, described in \autoref{sec:TRA_zpshde}.
 
-
 % -------------------------------------------------------------------------------------------------------------
 %         Rotated laplacian operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Standard and triad (bi-)laplacian operator}
+\subsection{Standard and triad (bi -)laplacian operator}
 \label{subsec:TRA_ldf_iso_triad}
 
-%&&    Standard rotated (bi-)laplacian operator
+%&&    Standard rotated (bi -)laplacian operator
 %&& ----------------------------------------------
-\subsubsection{Standard rotated (bi-)laplacian operator (\protect\mdl{traldf\_iso})}
+\subsubsection{Standard rotated (bi -)laplacian operator (\protect\mdl{traldf\_iso})}
 \label{subsec:TRA_ldf_iso}
-
 The general form of the second order lateral tracer subgrid scale physics (\autoref{eq:PE_zdf})
-takes the following semi-discrete space form in $z$- and $s$-coordinates:
+takes the following semi -discrete space form in $z$- and $s$-coordinates:
 \begin{equation}
   \label{eq:tra_ldf_iso}
   \begin{split}
-    D_T^{lT} = \frac{1}{b_t}   & \left\{   \,\;\delta_i \left[   A_u^{lT}  \left(
-          \frac{e_{2u}\,e_{3u}}{e_{1u}} \,\delta_{i+1/2}[T]
-          - e_{2u}\;r_{1u} \,\overline{\overline{ \delta_{k+1/2}[T] }}^{\,i+1/2,k}
-        \right)   \right]   \right.    \\
-    &             +\delta_j \left[ A_v^{lT} \left(
-        \frac{e_{1v}\,e_{3v}}{e_{2v}}  \,\delta_{j+1/2} [T]
-        - e_{1v}\,r_{2v} \,\overline{\overline{ \delta_{k+1/2} [T] }}^{\,j+1/2,k}
-      \right)   \right]                 \\
-    & +\delta_k \left[ A_w^{lT} \left(
-        -\;e_{2w}\,r_{1w} \,\overline{\overline{ \delta_{i+1/2} [T] }}^{\,i,k+1/2}
-      \right.   \right.                 \\
-    & \qquad \qquad \quad
-    - e_{1w}\,r_{2w} \,\overline{\overline{ \delta_{j+1/2} [T] }}^{\,j,k+1/2}     \\
-    & \left. {\left. {   \qquad \qquad \ \ \ \left. {
-                +\;\frac{e_{1w}\,e_{2w}}{e_{3w}} \,\left( r_{1w}^2 + r_{2w}^2 \right)
-                \,\delta_{k+1/2} [T] } \right) } \right] \quad } \right\}
+    D_T^{lT} = \frac{1}{b_t} \Bigg[ \quad &\delta_i A_u^{lT} \lt( \frac{e_{2u} e_{3u}}{e_{1u}}                      \, \delta_{i + 1/2} [T]
+                                                                  - e_{2u} r_{1u} \, \overline{\overline{\delta_{k + 1/2} [T]}}^{\,i + 1/2,k} \rt) \Bigg. \\
+                                    +     &\delta_j A_v^{lT} \lt( \frac{e_{1v} e_{3v}}{e_{2v}}                       \, \delta_{j + 1/2} [T]
+                                                                  - e_{1v} r_{2v} \, \overline{\overline{\delta_{k + 1/2} [T]}}^{\,j + 1/2,k} \rt)        \\
+                                   +     &\delta_k A_w^{lT} \lt( \frac{e_{1w} e_{2w}}{e_{3w}} (r_{1w}^2 + r_{2w}^2) \, \delta_{k + 1/2} [T] \rt.           \\
+                                   & \qquad \quad \Bigg. \lt.     - e_{2w} r_{1w} \, \overline{\overline{\delta_{i + 1/2} [T]}}^{\,i,k + 1/2}
+                                                                  - e_{1w} r_{2w} \, \overline{\overline{\delta_{j + 1/2} [T]}}^{\,j,k + 1/2} \rt) \Bigg]
   \end{split}
 \end{equation}
-where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells,
+where $b_t = e_{1t} \, e_{2t} \, e_{3t}$  is the volume of $T$-cells,
 $r_1$ and $r_2$ are the slopes between the surface of computation ($z$- or $s$-surfaces) and
 the surface along which the diffusion operator acts (\ie horizontal or iso-neutral surfaces).
-It is thus used when, in addition to \np{ln\_traldf\_lap}\forcode{ = .true.},
-we have \np{ln\_traldf\_iso}\forcode{ = .true.},
-or both \np{ln\_traldf\_hor}\forcode{ = .true.} and \np{ln\_zco}\forcode{ = .true.}.
+It is thus used when, in addition to \np{ln\_traldf\_lap}~\forcode{= .true.},
+we have \np{ln\_traldf\_iso}~\forcode{= .true.},
+or both \np{ln\_traldf\_hor}~\forcode{= .true.} and \np{ln\_zco}~\forcode{= .true.}.
 The way these slopes are evaluated is given in \autoref{sec:LDF_slp}.
 At the surface, bottom and lateral boundaries, the turbulent fluxes of heat and salt are set to zero using
-the mask technique (see \autoref{sec:LBC_coast}). 
+the mask technique (see \autoref{sec:LBC_coast}).
 
 The operator in \autoref{eq:tra_ldf_iso} involves both lateral and vertical derivatives.
@@ -606,24 +586,23 @@
 For computer efficiency reasons, this term is not computed in the \mdl{traldf\_iso} module,
 but in the \mdl{trazdf} module where, if iso-neutral mixing is used,
-the vertical mixing coefficient is simply increased by
-$\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$. 
+the vertical mixing coefficient is simply increased by $\frac{e_{1w} e_{2w}}{e_{3w}}(r_{1w}^2 + r_{2w}^2)$.
 
 This formulation conserves the tracer but does not ensure the decrease of the tracer variance.
 Nevertheless the treatment performed on the slopes (see \autoref{chap:LDF}) allows the model to run safely without
-any additional background horizontal diffusion \citep{Guilyardi_al_CD01}. 
-
-Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}),
+any additional background horizontal diffusion \citep{Guilyardi_al_CD01}.
+
+Note that in the partial step $z$-coordinate (\np{ln\_zps}~\forcode{= .true.}),
 the horizontal derivatives at the bottom level in \autoref{eq:tra_ldf_iso} require a specific treatment.
 They are calculated in module zpshde, described in \autoref{sec:TRA_zpshde}.
 
-%&&     Triad rotated (bi-)laplacian operator
+%&&     Triad rotated (bi -)laplacian operator
 %&&  -------------------------------------------
-\subsubsection{Triad rotated (bi-)laplacian operator (\protect\np{ln\_traldf\_triad})}
+\subsubsection{Triad rotated (bi -)laplacian operator (\protect\np{ln\_traldf\_triad})}
 \label{subsec:TRA_ldf_triad}
 
-If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}\forcode{ = .true.}; see \autoref{apdx:triad})
+If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}~\forcode{= .true.}; see \autoref{apdx:triad})
 
 An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
-is also available in \NEMO (\np{ln\_traldf\_grif}\forcode{ = .true.}).
+is also available in \NEMO (\np{ln\_traldf\_grif}~\forcode{= .true.}).
 A complete description of the algorithm is given in \autoref{apdx:triad}.
 
@@ -635,5 +614,5 @@
 It requires an additional assumption on boundary conditions:
 first and third derivative terms normal to the coast,
-normal to the bottom and normal to the surface are set to zero. 
+normal to the bottom and normal to the surface are set to zero.
 
 %&&    Option for the rotated operators
@@ -642,14 +621,11 @@
 \label{subsec:TRA_ldf_options}
 
-\np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators)
-
-\np{rn\_slpmax} = slope limit (both operators)
-
-\np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only)
-
-\np{rn\_sw\_triad} =1 switching triad;
-                   =0 all 4 triads used (triad only) 
-
-\np{ln\_botmix\_triad} = lateral mixing on bottom (triad only)
+\begin{itemize}
+\item \np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators)
+\item \np{rn\_slpmax} = slope limit (both operators)
+\item \np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only)
+\item \np{rn\_sw\_triad} $= 1$ switching triad; $= 0$ all 4 triads used (triad only) 
+\item \np{ln\_botmix\_triad} = lateral mixing on bottom (triad only)
+\end{itemize}
 
 % ================================================================
@@ -666,33 +642,32 @@
 The formulation of the vertical subgrid scale tracer physics is the same for all the vertical coordinates,
 and is based on a laplacian operator.
-The vertical diffusion operator given by (\autoref{eq:PE_zdf}) takes the following semi-discrete space form:
-\[
+The vertical diffusion operator given by (\autoref{eq:PE_zdf}) takes the following semi -discrete space form:
+\begin{gather*}
   % \label{eq:tra_zdf}
-  \begin{split}
-    D^{vT}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vT}_w}{e_{3w}}  \delta_{k+1/2}[T] \;\right]    \\
-    D^{vS}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vS}_w}{e_{3w}}  \delta_{k+1/2}[S] \;\right]
-  \end{split}
-\]
+    D^{vT}_T = \frac{1}{e_{3t}} \, \delta_k \lt[ \, \frac{A^{vT}_w}{e_{3w}} \delta_{k + 1/2}[T] \, \rt] \\
+    D^{vS}_T = \frac{1}{e_{3t}} \; \delta_k \lt[ \, \frac{A^{vS}_w}{e_{3w}} \delta_{k + 1/2}[S] \, \rt]
+\end{gather*}
 where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity coefficients on temperature and salinity,
 respectively.
-Generally, $A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is parameterised (\ie \key{zdfddm} is defined).
+Generally, $A_w^{vT} = A_w^{vS}$ except when double diffusive mixing is parameterised
+(\ie \key{zdfddm} is defined).
 The way these coefficients are evaluated is given in \autoref{chap:ZDF} (ZDF).
 Furthermore, when iso-neutral mixing is used, both mixing coefficients are increased by
-$\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ to account for
-the vertical second derivative of \autoref{eq:tra_ldf_iso}.
+$\frac{e_{1w} e_{2w}}{e_{3w} }({r_{1w}^2 + r_{2w}^2})$ to account for the vertical second derivative of
+\autoref{eq:tra_ldf_iso}.
 
 At the surface and bottom boundaries, the turbulent fluxes of heat and salt must be specified.
 At the surface they are prescribed from the surface forcing and added in a dedicated routine
 (see \autoref{subsec:TRA_sbc}), whilst at the bottom they are set to zero for heat and salt unless
-a geothermal flux forcing is prescribed as a bottom boundary condition (see \autoref{subsec:TRA_bbc}). 
+a geothermal flux forcing is prescribed as a bottom boundary condition (see \autoref{subsec:TRA_bbc}).
 
 The large eddy coefficient found in the mixed layer together with high vertical resolution implies that
-in the case of explicit time stepping (\np{ln\_zdfexp}\forcode{ = .true.})
+in the case of explicit time stepping (\np{ln\_zdfexp}~\forcode{= .true.})
 there would be too restrictive a constraint on the time step.
 Therefore, the default implicit time stepping is preferred for the vertical diffusion since
 it overcomes the stability constraint.
-A forward time differencing scheme (\np{ln\_zdfexp}\forcode{ = .true.}) using
+A forward time differencing scheme (\np{ln\_zdfexp}~\forcode{= .true.}) using
 a time splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative.
-Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 
+Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.
 
 % ================================================================
@@ -712,5 +687,5 @@
 This has been found to enhance readability of the code.
 The two formulations are completely equivalent;
-the forcing terms in trasbc are the surface fluxes divided by the thickness of the top model layer. 
+the forcing terms in trasbc are the surface fluxes divided by the thickness of the top model layer.
 
 Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components
@@ -724,32 +699,36 @@
 The surface module (\mdl{sbcmod}, see \autoref{chap:SBC}) provides the following forcing fields (used on tracers):
 
-$\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface
-(\ie the difference between the total surface heat flux and the fraction of the short wave flux that 
-penetrates into the water column, see \autoref{subsec:TRA_qsr})
-plus the heat content associated with of the mass exchange with the atmosphere and lands.
-
-$\bullet$ $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...)
-
-$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) and
-possibly with the sea-ice and ice-shelves.
-
-$\bullet$ \textit{rnf}, the mass flux associated with runoff
-(see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
-
-$\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt,
-(see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied).
+\begin{itemize}
+\item
+  $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface
+  (\ie the difference between the total surface heat flux and the fraction of the short wave flux that
+  penetrates into the water column, see \autoref{subsec:TRA_qsr})
+  plus the heat content associated with of the mass exchange with the atmosphere and lands.
+\item
+  $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...)
+\item
+  \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) and
+  possibly with the sea-ice and ice-shelves.
+\item
+  \textit{rnf}, the mass flux associated with runoff
+  (see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
+\item
+  \textit{fwfisf}, the mass flux associated with ice shelf melt,
+  (see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied).
+\end{itemize}
 
 The surface boundary condition on temperature and salinity is applied as follows:
 \begin{equation}
   \label{eq:tra_sbc}
-  \begin{aligned}
-    &F^T = \frac{ 1 }{\rho_o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns}       }^t  & \\
-    & F^S =\frac{ 1 }{\rho_o  \,      \left. e_{3t} \right|_{k=1} }  &\overline{ \textit{sfx} }^t   & \\
-  \end{aligned}
-\end{equation} 
-where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps ($t-\rdt/2$ and $t+\rdt/2$).
+  \begin{alignedat}{2}
+    F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} &\overline{Q_{ns}      }^t \\
+    F^S &=               &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} &\overline{\textit{sfx}}^t
+  \end{alignedat}
+\end{equation}
+where $\overline x^t$ means that $x$ is averaged over two consecutive time steps
+($t - \rdt / 2$ and $t + \rdt / 2$).
 Such time averaging prevents the divergence of odd and even time step (see \autoref{chap:STP}).
 
-In the linear free surface case (\np{ln\_linssh}\forcode{ = .true.}), an additional term has to be added on
+In the linear free surface case (\np{ln\_linssh}~\forcode{= .true.}), an additional term has to be added on
 both temperature and salinity.
 On temperature, this term remove the heat content associated with mass exchange that has been added to $Q_{ns}$.
@@ -759,11 +738,10 @@
 \begin{equation}
   \label{eq:tra_sbc_lin}
-  \begin{aligned}
-    &F^T = \frac{ 1 }{\rho_o \;C_p \,\left. e_{3t} \right|_{k=1} }
-    &\overline{ \left( Q_{ns} - \textit{emp}\;C_p\,\left. T \right|_{k=1} \right) }^t  & \\
-    %
-    & F^S =\frac{ 1 }{\rho_o \,\left. e_{3t} \right|_{k=1} }
-    &\overline{ \left( \;\textit{sfx} - \textit{emp} \;\left. S \right|_{k=1}  \right) }^t   & \\
-  \end{aligned}
+  \begin{alignedat}{2}
+    F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}}
+          &\overline{(Q_{ns}       - C_p \, \textit{emp} \lt. T \rt|_{k = 1})}^t \\
+    F^S &=               &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}}
+          &\overline{(\textit{sfx} -        \textit{emp} \lt. S \rt|_{k = 1})}^t
+  \end{alignedat}
 \end{equation} 
 Note that an exact conservation of heat and salt content is only achieved with non-linear free surface.
@@ -783,7 +761,7 @@
 
 Options are defined through the \ngn{namtra\_qsr} namelist variables.
-When the penetrative solar radiation option is used (\np{ln\_flxqsr}\forcode{ = .true.}),
+When the penetrative solar radiation option is used (\np{ln\_flxqsr}~\forcode{= .true.}),
 the solar radiation penetrates the top few tens of meters of the ocean.
-If it is not used (\np{ln\_flxqsr}\forcode{ = .false.}) all the heat flux is absorbed in the first ocean level.
+If it is not used (\np{ln\_flxqsr}~\forcode{= .false.}) all the heat flux is absorbed in the first ocean level.
 Thus, in the former case a term is added to the time evolution equation of temperature \autoref{eq:PE_tra_T} and
 the surface boundary condition is modified to take into account only the non-penetrative part of the surface 
@@ -791,15 +769,15 @@
 \begin{equation}
   \label{eq:PE_qsr}
-  \begin{split}
-    \frac{\partial T}{\partial t} &= {\ldots} + \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k}	\\
-    Q_{ns} &= Q_\text{Total} - Q_{sr}
-  \end{split}
+  \begin{gathered}
+    \pd[T]{t} = \ldots + \frac{1}{\rho_o \, C_p \, e_3} \; \pd[I]{k} \\
+    Q_{ns} = Q_\text{Total} - Q_{sr}
+  \end{gathered}
 \end{equation}
 where $Q_{sr}$ is the penetrative part of the surface heat flux (\ie the shortwave radiation) and
-$I$ is the downward irradiance ($\left. I \right|_{z=\eta}=Q_{sr}$).
+$I$ is the downward irradiance ($\lt. I \rt|_{z = \eta} = Q_{sr}$).
 The additional term in \autoref{eq:PE_qsr} is discretized as follows:
 \begin{equation}
   \label{eq:tra_qsr}
-  \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k} \equiv \frac{1}{\rho_o\, C_p\, e_{3t}} \delta_k \left[ I_w \right]
+  \frac{1}{\rho_o \, C_p \, e_3} \, \pd[I]{k} \equiv \frac{1}{\rho_o \, C_p \, e_{3t}} \delta_k [I_w]
 \end{equation}
 
@@ -810,18 +788,18 @@
 (specified through namelist parameter \np{rn\_abs}).
 It is assumed to penetrate the ocean with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$,
-of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the \ngn{namtra\_qsr} namelist).
+of a few tens of centimetres (typically $\xi_0 = 0.35~m$ set as \np{rn\_si0} in the \ngn{namtra\_qsr} namelist).
 For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy propagates to
 larger depths where it contributes to local heating.
 The way this second part of the solar energy penetrates into the ocean depends on which formulation is chosen.
-In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}\forcode{ = .true.})
+In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}~\forcode{= .true.})
 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths,
 leading to the following expression \citep{Paulson1977}:
 \[
   % \label{eq:traqsr_iradiance}
-  I(z) = Q_{sr} \left[Re^{-z / \xi_0} + \left( 1-R\right) e^{-z / \xi_1} \right]
+  I(z) = Q_{sr} \lt[ Re^{- z / \xi_0} + (1 - R) e^{- z / \xi_1} \rt]
 \]
 where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths.
 It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter.
-The set of default values ($\xi_0$, $\xi_1$, $R$) corresponds to a Type I water in Jerlov's (1968) classification
+The set of default values ($\xi_0, \xi_1, R$) corresponds to a Type I water in Jerlov's (1968) classification
 (oligotrophic waters).
 
@@ -840,25 +818,27 @@
 reproduces quite closely the light penetration profiles predicted by the full spectal model,
 but with much greater computational efficiency.
-The 2-bands formulation does not reproduce the full model very well. 
-
-The RGB formulation is used when \np{ln\_qsr\_rgb}\forcode{ = .true.}.
+The 2-bands formulation does not reproduce the full model very well.
+
+The RGB formulation is used when \np{ln\_qsr\_rgb}~\forcode{= .true.}.
 The RGB attenuation coefficients (\ie the inverses of the extinction length scales) are tabulated over
 61 nonuniform chlorophyll classes ranging from 0.01 to 10 g.Chl/L
 (see the routine \rou{trc\_oce\_rgb} in \mdl{trc\_oce} module).
 Four types of chlorophyll can be chosen in the RGB formulation:
-\begin{description} 
-\item[\np{nn\_chdta}\forcode{ = 0}]
+
+\begin{description}
+\item[\np{nn\_chdta}~\forcode{= 0}]
   a constant 0.05 g.Chl/L value everywhere ; 
-\item[\np{nn\_chdta}\forcode{ = 1}]
+\item[\np{nn\_chdta}~\forcode{= 1}]
   an observed time varying chlorophyll deduced from satellite surface ocean color measurement spread uniformly in
   the vertical direction;
-\item[\np{nn\_chdta}\forcode{ = 2}]
+\item[\np{nn\_chdta}~\forcode{= 2}]
   same as previous case except that a vertical profile of chlorophyl is used.
   Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value;
-\item[\np{ln\_qsr\_bio}\forcode{ = .true.}]
+\item[\np{ln\_qsr\_bio}~\forcode{= .true.}]
   simulated time varying chlorophyll by TOP biogeochemical model.
   In this case, the RGB formulation is used to calculate both the phytoplankton light limitation in
-  PISCES or LOBSTER and the oceanic heating rate. 
+  PISCES or LOBSTER and the oceanic heating rate.
 \end{description} 
+
 The trend in \autoref{eq:tra_qsr} associated with the penetration of the solar radiation is added to
 the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}.
@@ -871,10 +851,10 @@
 Finally, note that when the ocean is shallow ($<$ 200~m), part of the solar radiation can reach the ocean floor.
 In this case, we have chosen that all remaining radiation is absorbed in the last ocean level
-(\ie $I$ is masked). 
+(\ie $I$ is masked).
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[width=1.0\textwidth]{Fig_TRA_Irradiance}
+    \includegraphics[]{Fig_TRA_Irradiance}
     \caption{
       \protect\label{fig:traqsr_irradiance}
@@ -903,5 +883,5 @@
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[width=1.0\textwidth]{Fig_TRA_geoth}
+    \includegraphics[]{Fig_TRA_geoth}
     \caption{
       \protect\label{fig:geothermal}
@@ -917,8 +897,8 @@
 This is the default option in \NEMO, and it is implemented using the masking technique.
 However, there is a non-zero heat flux across the seafloor that is associated with solid earth cooling.
-This flux is weak compared to surface fluxes (a mean global value of $\sim0.1\;W/m^2$ \citep{Stein_Stein_Nat92}),
+This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{Stein_Stein_Nat92}),
 but it warms systematically the ocean and acts on the densest water masses.
 Taking this flux into account in a global ocean model increases the deepest overturning cell
-(\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups  \citep{Emile-Geay_Madec_OS09}. 
+(\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OS09}.
 
 Options are defined through the  \ngn{namtra\_bbc} namelist variables.
@@ -939,5 +919,5 @@
 %--------------------------------------------------------------------------------------------------------------
 
-Options are defined through the  \ngn{nambbl} namelist variables.
+Options are defined through the \ngn{nambbl} namelist variables.
 In a $z$-coordinate configuration, the bottom topography is represented by a series of discrete steps.
 This is not adequate to represent gravity driven downslope flows.
@@ -951,5 +931,5 @@
 sometimes over a thickness much larger than the thickness of the observed gravity plume.
 A similar problem occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly downstream of
-a sill \citep{Willebrand_al_PO01}, and the thickness of the plume is not resolved. 
+a sill \citep{Willebrand_al_PO01}, and the thickness of the plume is not resolved.
 
 The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{Beckmann_Doscher1997},
@@ -964,5 +944,5 @@
 %        Diffusive BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}\forcode{ = 1})}
+\subsection{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}~\forcode{= 1})}
 \label{subsec:TRA_bbl_diff}
 
@@ -971,19 +951,19 @@
 \[
   % \label{eq:tra_bbl_diff}
-  {\rm {\bf F}}_\sigma=A_l^\sigma \; \nabla_\sigma T
+  \vect F_\sigma = A_l^\sigma \, \nabla_\sigma T
 \]
-with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells,
-and  $A_l^\sigma$ the lateral diffusivity in the BBL.
+with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, and
+$A_l^\sigma$ the lateral diffusivity in the BBL.
 Following \citet{Beckmann_Doscher1997}, the latter is prescribed with a spatial dependence,
 \ie in the conditional form
 \begin{equation}
   \label{eq:tra_bbl_coef}
-  A_l^\sigma (i,j,t)=\left\{ {
-      \begin{array}{l}
-        A_{bbl}  \quad \quad   \mbox{if}  \quad   \nabla_\sigma \rho  \cdot  \nabla H<0 \\ \\
-        0\quad \quad \;\,\mbox{otherwise} \\
-      \end{array}}
-  \right.
-\end{equation} 
+  A_l^\sigma (i,j,t) =
+      \begin{cases}
+        A_{bbl} & \text{if~} \nabla_\sigma \rho \cdot \nabla H < 0 \\
+        \\
+        0      & \text{otherwise} \\
+      \end{cases}
+\end{equation}
 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist parameter \np{rn\_ahtbbl} and
 usually set to a value much larger than the one used for lateral mixing in the open ocean.
@@ -995,26 +975,28 @@
 \[
   % \label{eq:tra_bbl_Drho}
-  \nabla_\sigma \rho / \rho = \alpha \,\nabla_\sigma T + \beta   \,\nabla_\sigma S
+  \nabla_\sigma \rho / \rho = \alpha \, \nabla_\sigma T + \beta \, \nabla_\sigma S
 \]
-where $\rho$, $\alpha$ and $\beta$ are functions of $\overline{T}^\sigma$,
-$\overline{S}^\sigma$ and $\overline{H}^\sigma$, the along bottom mean temperature, salinity and depth, respectively.
+where $\rho$, $\alpha$ and $\beta$ are functions of $\overline T^\sigma$, $\overline S^\sigma$ and
+$\overline H^\sigma$, the along bottom mean temperature, salinity and depth, respectively.
 
 % -------------------------------------------------------------------------------------------------------------
 %        Advective BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Advective bottom boundary layer  (\protect\np{nn\_bbl\_adv}\forcode{ = 1..2})}
+\subsection{Advective bottom boundary layer  (\protect\np{nn\_bbl\_adv}~\forcode{= 1..2})}
 \label{subsec:TRA_bbl_adv}
 
-%\sgacomment{"downsloping flow" has been replaced by "downslope flow" in the following
-%if this is not what is meant then "downwards sloping flow" is also a possibility"}
+%\sgacomment{
+%  "downsloping flow" has been replaced by "downslope flow" in the following
+%  if this is not what is meant then "downwards sloping flow" is also a possibility"
+%}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[width=0.7\textwidth]{Fig_BBL_adv}
+    \includegraphics[]{Fig_BBL_adv}
     \caption{
       \protect\label{fig:bbl}
       Advective/diffusive Bottom Boundary Layer.
-      The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$.
+      The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i + 1}_{kdnw}$.
       Red arrows indicate the additional overturning circulation due to the advective BBL.
       The transport of the downslope flow is defined either as the transport of the bottom ocean cell (black arrow),
@@ -1026,32 +1008,31 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-
 %!!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
 %!!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
-%!!        \ie transport proportional to the along-slope density gradient
+%!!        i.e. transport proportional to the along-slope density gradient
 
 %%%gmcomment   :  this section has to be really written
 
-When applying an advective BBL (\np{nn\_bbl\_adv}\forcode{ = 1..2}), an overturning circulation is added which
+When applying an advective BBL (\np{nn\_bbl\_adv}~\forcode{= 1..2}), an overturning circulation is added which
 connects two adjacent bottom grid-points only if dense water overlies less dense water on the slope.
-The density difference causes dense water to move down the slope. 
-
-\np{nn\_bbl\_adv}\forcode{ = 1}:
+The density difference causes dense water to move down the slope.
+
+\np{nn\_bbl\_adv}~\forcode{= 1}:
 the downslope velocity is chosen to be the Eulerian ocean velocity just above the topographic step
 (see black arrow in \autoref{fig:bbl}) \citep{Beckmann_Doscher1997}.
 It is a \textit{conditional advection}, that is, advection is allowed only
-if dense water overlies less dense water on the slope (\ie $\nabla_\sigma \rho  \cdot  \nabla H<0$) and
-if the velocity is directed towards greater depth (\ie $\vect{U}  \cdot  \nabla H>0$).
-
-\np{nn\_bbl\_adv}\forcode{ = 2}:
+if dense water overlies less dense water on the slope (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$) and
+if the velocity is directed towards greater depth (\ie $\vect U \cdot \nabla H > 0$).
+
+\np{nn\_bbl\_adv}~\forcode{= 2}:
 the downslope velocity is chosen to be proportional to $\Delta \rho$,
 the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}.
 The advection is allowed only  if dense water overlies less dense water on the slope
-(\ie $\nabla_\sigma \rho  \cdot  \nabla H<0$).
+(\ie $\nabla_\sigma \rho \cdot \nabla H < 0$).
 For example, the resulting transport of the downslope flow, here in the $i$-direction (\autoref{fig:bbl}),
 is simply given by the following expression:
 \[
   % \label{eq:bbl_Utr}
-  u^{tr}_{bbl} = \gamma \, g \frac{\Delta \rho}{\rho_o}  e_{1u} \; min \left( {e_{3u}}_{kup},{e_{3u}}_{kdwn} \right)
+  u^{tr}_{bbl} = \gamma g \frac{\Delta \rho}{\rho_o} e_{1u} \, min ({e_{3u}}_{kup},{e_{3u}}_{kdwn})
 \]
 where $\gamma$, expressed in seconds, is the coefficient of proportionality provided as \np{rn\_gambbl},
@@ -1062,5 +1043,5 @@
 The possible values for $\gamma$ range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}.
 
-Scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$ using the upwind scheme.
+Scalar properties are advected by this additional transport $(u^{tr}_{bbl},v^{tr}_{bbl})$ using the upwind scheme.
 Such a diffusive advective scheme has been chosen to mimic the entrainment between the downslope plume and
 the surrounding water at intermediate depths.
@@ -1071,19 +1052,21 @@
 the downslope flow \autoref{eq:bbl_dw}, the horizontal \autoref{eq:bbl_hor} and
 the upward \autoref{eq:bbl_up} return flows as follows: 
-\begin{align}
+\begin{alignat}{3}
+  \label{eq:bbl_dw}
   \partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw}
-                            +  \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}}  \left( T^{sh}_{kup} - T^{do}_{kdw} \right)  \label{eq:bbl_dw} \\
-                            %
+                                &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}} &&\lt( T^{sh}_{kup} - T^{do}_{kdw} \rt) \\
+  \label{eq:bbl_hor}
   \partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup}
-                            + \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}}   \left( T^{do}_{kup} - T^{sh}_{kup} \right)   \label{eq:bbl_hor} \\
-                            %
+                                &&+ \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}} &&\lt( T^{do}_{kup} - T^{sh}_{kup} \rt) \\
+  %
   \intertext{and for $k =kdw-1,\;..., \; kup$ :}
   %
+  \label{eq:bbl_up}
   \partial_t T^{do}_{k} &\equiv \partial_t S^{do}_{k}
-                          + \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   \left( T^{do}_{k+1} - T^{sh}_{k} \right)   \label{eq:bbl_up}
-\end{align}
-where $b_t$ is the $T$-cell volume. 
-
-Note that the BBL transport, $( u^{tr}_{bbl}, v^{tr}_{bbl} )$, is available in the model outputs.
+                                &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   &&\lt( T^{do}_{k +1} - T^{sh}_{k}   \rt)
+\end{alignat}
+where $b_t$ is the $T$-cell volume.
+
+Note that the BBL transport, $(u^{tr}_{bbl},v^{tr}_{bbl})$, is available in the model outputs.
 It has to be used to compute the effective velocity as well as the effective overturning circulation.
 
@@ -1101,8 +1084,8 @@
 \begin{equation}
   \label{eq:tra_dmp}
-  \begin{split}
-    \frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right)  \\
-    \frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right)
-  \end{split}
+  \begin{gathered}
+    \pd[T]{t} = \cdots - \gamma (T - T_o) \\
+    \pd[S]{t} = \cdots - \gamma (S - S_o)
+  \end{gathered}
 \end{equation} 
 where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ are given temperature and salinity fields
@@ -1111,5 +1094,5 @@
 The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true.
 It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_tradmp} are set to true in
-\textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are correctly set
+\ngn{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are correctly set
 (\ie that $T_o$ and $S_o$ are provided in input files and read using \mdl{fldread},
 see \autoref{subsec:SBC_fldread}).
@@ -1128,5 +1111,5 @@
 The second case corresponds to the use of the robust diagnostic method \citep{Sarmiento1982}.
 It allows us to find the velocity field consistent with the model dynamics whilst
-having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$). 
+having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$).
 
 The robust diagnostic method is very efficient in preventing temperature drift in intermediate waters but
@@ -1140,52 +1123,6 @@
 \citep{Madec_al_JPO96}.
 
-\subsection{Generating \ifile{resto} using DMP\_TOOLS}
-
-DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$.
-Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled and
-run on the same machine as the NEMO model.
-A \ifile{mesh\_mask} file for the model configuration is required as an input.
-This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1.
-The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work.
-The \ngn{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for
-the restoration coefficient.
-
-%--------------------------------------------nam_dmp_create-------------------------------------------------
-%\namtools{namelist_dmp}
-%-------------------------------------------------------------------------------------------------------
-
-\np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and
-should be the same as specified in \ngn{namcfg}.
-The variable \np{lzoom} is used to specify that the damping is being used as in case \textit{a} above to
-provide boundary conditions to a zoom configuration.
-In the case of the arctic or antarctic zoom configurations this includes some specific treatment.
-Otherwise damping is applied to the 6 grid points along the ocean boundaries.
-The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in
-the \ngn{nam\_zoom\_dmp} name list.
-
-The remaining switch namelist variables determine the spatial variation of the restoration coefficient in
-non-zoom configurations.
-\np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain.
-\np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea for
-the ORCA4, ORCA2 and ORCA05 configurations.
-If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as
-a function of the model number.
-This option is included to allow backwards compatability of the ORCA2 reference configurations with
-previous model versions.
-\np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines.
-This option only has an effect if \np{ln\_full\_field} is true.
-\np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer.
-Finally \np{ln\_custom} specifies that the custom module will be called.
-This module is contained in the file \mdl{custom} and can be edited by users.
-For example damping could be applied in a specific region.
-
-The restoration coefficient can be set to zero in equatorial regions by
-specifying a positive value of \np{nn\_hdmp}. 
-Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to 
-the full values of a 10\deg latitud band. 
-This is often used because of the short adjustment time scale in the equatorial region
-\citep{Reverdin1991, Fujio1991, Marti_PhD92}.
-The time scale associated with the damping depends on the depth as a hyperbolic tangent,
-with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.  
+For generating \ifile{resto}, see the documentation for the DMP tool provided with the source code under
+\path{./tools/DMP_TOOLS}.
 
 % ================================================================
@@ -1199,15 +1136,14 @@
 %--------------------------------------------------------------------------------------------------------------
 
-Options are defined through the  \ngn{namdom} namelist variables.
+Options are defined through the \ngn{namdom} namelist variables.
 The general framework for tracer time stepping is a modified leap-frog scheme \citep{Leclair_Madec_OM09},
 \ie a three level centred time scheme associated with a Asselin time filter (cf. \autoref{sec:STP_mLF}):
 \begin{equation}
   \label{eq:tra_nxt}
-  \begin{aligned}
-    (e_{3t}T)^{t+\rdt} &= (e_{3t}T)_f^{t-\rdt} &+ 2 \, \rdt  \,e_{3t}^t\ \text{RHS}^t &	\\ \\
-    (e_{3t}T)_f^t  \;\ \quad &= (e_{3t}T)^t \;\quad
-    &+\gamma \,\left[ {(e_{3t}T)_f^{t-\rdt} -2(e_{3t}T)^t+(e_{3t}T)^{t+\rdt}} \right] &  \\
-    & &- \gamma\,\rdt \, \left[ Q^{t+\rdt/2} -  Q^{t-\rdt/2} \right]  &
-  \end{aligned}
+  \begin{alignedat}{3}
+    &(e_{3t}T)^{t + \rdt} &&= (e_{3t}T)_f^{t - \rdt} &&+ 2 \, \rdt \,e_{3t}^t \ \text{RHS}^t \\
+    &(e_{3t}T)_f^t        &&= (e_{3t}T)^t            &&+ \, \gamma \, \lt[ (e_{3t}T)_f^{t - \rdt} - 2(e_{3t}T)^t + (e_{3t}T)^{t + \rdt} \rt] \\
+    &                     &&                         &&- \, \gamma \, \rdt \, \lt[ Q^{t + \rdt/2} - Q^{t - \rdt/2} \rt]  
+  \end{alignedat}
 \end{equation} 
 where RHS is the right hand side of the temperature equation, the subscript $f$ denotes filtered values,
@@ -1215,15 +1151,15 @@
 (\ie fluxes plus content in mass exchanges).
 $\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter).
-Its default value is \np{rn\_atfp}\forcode{ = 10.e-3}.
+Its default value is \np{rn\_atfp}~\forcode{= 10.e-3}.
 Note that the forcing correction term in the filter is not applied in linear free surface
-(\jp{lk\_vvl}\forcode{ = .false.}) (see \autoref{subsec:TRA_sbc}.
+(\jp{lk\_vvl}~\forcode{= .false.}) (see \autoref{subsec:TRA_sbc}).
 Not also that in constant volume case, the time stepping is performed on $T$, not on its content, $e_{3t}T$.
 
-When the vertical mixing is solved implicitly,
-the update of the \textit{next} tracer fields is done in module \mdl{trazdf}.
+When the vertical mixing is solved implicitly, the update of the \textit{next} tracer fields is done in
+\mdl{trazdf} module.
 In this case only the swapping of arrays and the Asselin filtering is done in the \mdl{tranxt} module.
 
 In order to prepare for the computation of the \textit{next} time step, a swap of tracer arrays is performed:
-$T^{t-\rdt} = T^t$ and $T^t = T_f$. 
+$T^{t - \rdt} = T^t$ and $T^t = T_f$.
 
 % ================================================================
@@ -1240,5 +1176,5 @@
 %        Equation of State
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Equation of seawater (\protect\np{nn\_eos}\forcode{ = -1..1})}
+\subsection{Equation of seawater (\protect\np{nn\_eos}~\forcode{= -1..1})}
 \label{subsec:TRA_eos}
 
@@ -1264,5 +1200,5 @@
 To that purposed, a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available.
 
-In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density.
+In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density.
 Called \textit{rau0} in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$.
 This is a sensible choice for the reference density used in a Boussinesq ocean climate model, as,
@@ -1270,8 +1206,9 @@
 density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}.
 
-Options are defined through the  \ngn{nameos} namelist variables, and in particular \np{nn\_eos} which
+Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} which
 controls the EOS used (\forcode{= -1} for TEOS10 ; \forcode{= 0} for EOS-80 ; \forcode{= 1} for S-EOS).
+
 \begin{description}
-\item[\np{nn\_eos}\forcode{ = -1}]
+\item[\np{nn\_eos}~\forcode{= -1}]
   the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used.
   The accuracy of this approximation is comparable to the TEOS-10 rational function approximation,
@@ -1282,9 +1219,8 @@
   the TEOS-10 rational function approximation for hydrographic data analysis \citep{TEOS10}.
   A key point is that conservative state variables are used:
-  Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \deg{C}, notation: $\Theta$).
+  Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \degC, notation: $\Theta$).
   The pressure in decibars is approximated by the depth in meters.
   With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant.
-  It is set to $C_p=3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}.
-
+  It is set to $C_p = 3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}.
   Choosing polyTEOS10-bsq implies that the state variables used by the model are $\Theta$ and $S_A$.
   In particular, the initial state deined by the user have to be given as \textit{Conservative} Temperature and
@@ -1293,6 +1229,5 @@
   either computing the air-sea and ice-sea fluxes (forced mode) or
   sending the SST field to the atmosphere (coupled mode).
-
-\item[\np{nn\_eos}\forcode{ = 0}]
+\item[\np{nn\_eos}~\forcode{= 0}]
   the polyEOS80-bsq equation of seawater is used.
   It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized to
@@ -1305,7 +1240,6 @@
   pressure \citep{UNESCO1983}.
   Nevertheless, a severe assumption is made in order to have a heat content ($C_p T_p$) which
-  is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value. 
- 
-\item[\np{nn\_eos}\forcode{ = 1}]
+  is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value.
+\item[\np{nn\_eos}~\forcode{= 1}]
   a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen,
   the coefficients of which has been optimized to fit the behavior of TEOS10
@@ -1317,34 +1251,44 @@
   as well as between \textit{absolute} and \textit{practical} salinity.
   S-EOS takes the following expression:
-  \[
+  \begin{gather*}
     % \label{eq:tra_S-EOS}
-    \begin{split}
-      d_a(T,S,z)  =  ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a  \\
-      & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a  \\
-      & - \nu \; T_a \; S_a \;  ) \; / \; \rho_o                     \\
-      with \ \  T_a = T-10  \; ;  & \;  S_a = S-35  \; ;\;  \rho_o = 1026~Kg/m^3
-    \end{split}
-  \]
+    \begin{alignedat}{2}
+    &d_a(T,S,z) = \frac{1}{\rho_o} \big[ &- a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * &T_a \big. \\
+    &                                    &+ b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * &S_a       \\  
+    &                              \big. &- \nu \;                           T_a                  &S_a \big] \\
+    \end{alignedat}
+    \\
+    \text{with~} T_a = T - 10 \, ; \, S_a = S - 35 \, ; \, \rho_o = 1026~Kg/m^3
+  \end{gather*}
   where the computer name of the coefficients as well as their standard value are given in \autoref{tab:SEOS}.
-  In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing the associated coefficients.
-  Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS.
-  setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS.
+  In fact, when choosing S-EOS, various approximation of EOS can be specified simply by
+  changing the associated coefficients.
+  Setting to zero the two thermobaric coefficients $(\mu_1,\mu_2)$ remove thermobaric effect from S-EOS.
+  setting to zero the three cabbeling coefficients $(\lambda_1,\lambda_2,\nu)$ remove cabbeling effect from
+  S-EOS.
   Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S.
 \end{description}
-
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{table}[!tb]
   \begin{center}
-    \begin{tabular}{|p{26pt}|p{72pt}|p{56pt}|p{136pt}|}
+    \begin{tabular}{|l|l|l|l|}
       \hline
-      coeff.	& computer name   & S-EOS		&  description		                  \\ \hline
-      $a_0$       & \np{rn\_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff. 	\\ \hline
-      $b_0$	      & \np{rn\_b0} 	   & 7.6554 $10^{-1}$ &  linear haline  expansion coeff. 	\\ \hline
-      $\lambda_1$	& \np{rn\_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$ 	      \\ \hline
-      $\lambda_2$	& \np{rn\_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$	 	   \\ \hline
-      $\nu$       & \np{rn\_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$ 	   \\ \hline
-      $\mu_1$     & \np{rn\_mu1} 	& 1.4970 $10^{-4}$ &  thermobaric coeff. in T    	   \\ \hline
-      $\mu_2$     & \np{rn\_mu2} 	& 1.1090 $10^{-5}$ &  thermobaric coeff. in S   	      \\ \hline
+      coeff.	   & computer name   & S-EOS		      & description		                 \\
+      \hline
+      $a_0$       & \np{rn\_a0}     & $1.6550~10^{-1}$ & linear thermal expansion coeff. \\
+      \hline
+      $b_0$	      & \np{rn\_b0} 	   & $7.6554~10^{-1}$ & linear haline  expansion coeff. \\
+      \hline
+      $\lambda_1$	& \np{rn\_lambda1}& $5.9520~10^{-2}$ & cabbeling coeff. in $T^2$ 	     \\
+      \hline
+      $\lambda_2$	& \np{rn\_lambda2}& $5.4914~10^{-4}$ & cabbeling coeff. in $S^2$	 	  \\
+      \hline
+      $\nu$       & \np{rn\_nu}     & $2.4341~10^{-3}$ & cabbeling coeff. in $T \, S$ 	  \\
+      \hline
+      $\mu_1$     & \np{rn\_mu1} 	& $1.4970~10^{-4}$ & thermobaric coeff. in T    	  \\
+      \hline
+      $\mu_2$     & \np{rn\_mu2} 	& $1.1090~10^{-5}$ & thermobaric coeff. in S   	     \\
+      \hline
     \end{tabular}
     \caption{
@@ -1352,16 +1296,15 @@
       Standard value of S-EOS coefficients.
     }
-  \end{center}
+\end{center}
 \end{table}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-
 % -------------------------------------------------------------------------------------------------------------
 %        Brunt-V\"{a}is\"{a}l\"{a} Frequency
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos}\forcode{ = 0..2})}
+\subsection{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos}~\forcode{= 0..2})}
 \label{subsec:TRA_bn2}
 
-An accurate computation of the ocean stability (\ie of $N$, the brunt-V\"{a}is\"{a}l\"{a} frequency) is of
+An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a} frequency) is of
 paramount importance as determine the ocean stratification and is used in several ocean parameterisations
 (namely TKE, GLS, Richardson number dependent vertical diffusion, enhanced vertical diffusion,
@@ -1372,10 +1315,10 @@
 \[
   % \label{eq:tra_bn2}
-  N^2 =	\frac{g}{e_{3w}} \left(   \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T]   \right)
+  N^2 = \frac{g}{e_{3w}} \lt( \beta \; \delta_{k + 1/2}[S] - \alpha \; \delta_{k + 1/2}[T] \rt)
 \]
-where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS,
-and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients.
-The coefficients are a polynomial function of temperature, salinity and depth which
-expression depends on the chosen EOS.
+where $(T,S) = (\Theta,S_A)$ for TEOS10, $(\theta,S_p)$ for TEOS-80, or $(T,S)$ for S-EOS, and,
+$\alpha$ and $\beta$ are the thermal and haline expansion coefficients.
+The coefficients are a polynomial function of temperature, salinity and depth which expression depends on
+the chosen EOS.
 They are computed through \textit{eos\_rab}, a \fortran function that can be found in \mdl{eosbn2}.
 
@@ -1390,15 +1333,15 @@
   \label{eq:tra_eos_fzp}
   \begin{split}
-    T_f (S,p) = \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} -  2.154996 \;10^{-4} \,S  \right) \ S    \\
-    - 7.53\,10^{-3} \ \ p
-  \end{split}
+    &T_f (S,p) = \lt( a + b \, \sqrt{S} + c \, S \rt) \, S + d \, p \\
+    &\text{where~} a = -0.0575, \, b = 1.710523~10^{-3}, \, c = -2.154996~10^{-4} \\ 
+    &\text{and~} d = -7.53~10^{-3}
+    \end{split}
 \end{equation}
 
 \autoref{eq:tra_eos_fzp} is only used to compute the potential freezing point of sea water
-(\ie referenced to the surface $p=0$),
+(\ie referenced to the surface $p = 0$),
 thus the pressure dependent terms in \autoref{eq:tra_eos_fzp} (last term) have been dropped.
 The freezing point is computed through \textit{eos\_fzp},
-a \fortran function that can be found in \mdl{eosbn2}.  
-
+a \fortran function that can be found in \mdl{eosbn2}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1411,5 +1354,4 @@
 %
 
-
 % ================================================================
 % Horizontal Derivative in zps-coordinate 
@@ -1421,9 +1363,9 @@
 I've changed "derivative" to "difference" and "mean" to "average"}
 
-With partial cells (\np{ln\_zps}\forcode{ = .true.}) at bottom and top (\np{ln\_isfcav}\forcode{ = .true.}),
+With partial cells (\np{ln\_zps}~\forcode{= .true.}) at bottom and top (\np{ln\_isfcav}~\forcode{= .true.}),
 in general, tracers in horizontally adjacent cells live at different depths.
 Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) and
 the hydrostatic pressure gradient calculations (\mdl{dynhpg} module).
-The partial cell properties at the top (\np{ln\_isfcav}\forcode{ = .true.}) are computed in the same way as
+The partial cell properties at the top (\np{ln\_isfcav}~\forcode{= .true.}) are computed in the same way as
 for the bottom.
 So, only the bottom interpolation is explained below.
@@ -1432,18 +1374,18 @@
 a linear interpolation in the vertical is used to approximate the deeper tracer as if
 it actually lived at the depth of the shallower tracer point (\autoref{fig:Partial_step_scheme}).
-For example, for temperature in the $i$-direction the needed interpolated temperature, $\widetilde{T}$, is:
+For example, for temperature in the $i$-direction the needed interpolated temperature, $\widetilde T$, is:
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!p]
   \begin{center}
-    \includegraphics[width=0.9\textwidth]{Fig_partial_step_scheme}
+    \includegraphics[]{Fig_partial_step_scheme}
     \caption{
       \protect\label{fig:Partial_step_scheme}
       Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate
-      (\protect\np{ln\_zps}\forcode{ = .true.}) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$.
-      A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$,
+      (\protect\np{ln\_zps}~\forcode{= .true.}) in the case $(e3w_k^{i + 1} - e3w_k^i) > 0$.
+      A linear interpolation is used to estimate $\widetilde T_k^{i + 1}$,
       the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points.
-      The horizontal difference is then given by: $\delta_{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and
-      the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$.
+      The horizontal difference is then given by: $\delta_{i + 1/2} T_k = \widetilde T_k^{\, i + 1} -T_k^{\, i}$ and
+      the average by: $\overline T_k^{\, i + 1/2} = (\widetilde T_k^{\, i + 1/2} - T_k^{\, i}) / 2$.
     }
   \end{center}
@@ -1451,42 +1393,43 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \[
-  \widetilde{T}= \left\{
-    \begin{aligned}
-      &T^{\,i+1}      -\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right)}{ e_{3w}^{i+1} }\;\delta_k T^{i+1}
-      && \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$   } 	\\ \\
-      &T^{\,i} \ \ \ \,+\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right) }{e_{3w}^i       }\;\delta_k T^{i+1}
-      && \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   }
-    \end{aligned}
-  \right.
+  \widetilde T = \lt\{
+    \begin{alignedat}{2}
+      &T^{\, i + 1} &-\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt) }{ e_{3w}^{i + 1} } \; \delta_k T^{i + 1}
+      & \quad \text{if $e_{3w}^{i + 1} \geq e_{3w}^i$} \\ \\
+      &T^{\, i}     &+\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt )}{e_{3w}^i       } \; \delta_k T^{i + 1}
+      & \quad \text{if $e_{3w}^{i + 1} <    e_{3w}^i$}
+    \end{alignedat}
+  \rt.
 \]
 and the resulting forms for the horizontal difference and the horizontal average value of $T$ at a $U$-point are: 
 \begin{equation}
   \label{eq:zps_hde}
-  \begin{aligned}
-    \delta_{i+1/2} T=
+  \begin{split}
+    \delta_{i + 1/2} T       &=
     \begin{cases}
-      \ \ \ \widetilde {T}\quad\ -T^i	 	& \ \ \quad\quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\ \\
-      \ \ \ T^{\,i+1}-\widetilde{T}		& \ \ \quad\quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   }
+                                \widetilde T - T^i          & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\
+                                \\
+                                T^{\, i + 1} - \widetilde T & \text{if~} e_{3w}^{i + 1} <    e_{3w}^i
     \end{cases}
-    \\ \\
-    \overline {T}^{\,i+1/2} \ =
+    \\
+    \overline T^{\, i + 1/2} &=
     \begin{cases}
-      ( \widetilde {T}\ \ \;\,-T^{\,i})	 / 2	& \;\ \ \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\ \\
-      ( T^{\,i+1}-\widetilde{T} ) / 2		& \;\ \ \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   }
+                                (\widetilde T - T^{\, i}   ) / 2 & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\
+                                \\
+                                (T^{\, i + 1} - \widetilde T) / 2 & \text{if~} e_{3w}^{i + 1} <   e_{3w}^i
     \end{cases}
-  \end{aligned}
+  \end{split}
 \end{equation}
 
 The computation of horizontal derivative of tracers as well as of density is performed once for all at
 each time step in \mdl{zpshde} module and stored in shared arrays to be used when needed.
-It has to be emphasized that the procedure used to compute the interpolated density, $\widetilde{\rho}$,
+It has to be emphasized that the procedure used to compute the interpolated density, $\widetilde \rho$,
 is not the same as that used for $T$ and $S$.
-Instead of forming a linear approximation of density, we compute $\widetilde{\rho }$ from the interpolated values of
+Instead of forming a linear approximation of density, we compute $\widetilde \rho$ from the interpolated values of
 $T$ and $S$, and the pressure at a $u$-point
-(in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos} ): 
+(in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos}): 
 \[
   % \label{eq:zps_hde_rho}
-  \widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u })
-  \quad \text{where }\  z_u = \min \left( {z_T^{i+1} ,z_T^i } \right)
+  \widetilde \rho = \rho (\widetilde T,\widetilde S,z_u) \quad \text{where~} z_u = \min \lt( z_T^{i + 1},z_T^i \rt)
 \]
 
