Changeset 10501

Timestamp:

2019-01-10T17:30:42+01:00 (5 years ago)

Author:

nicolasmartin

Message:

Global work on math environnements for equations (partial commits)

• Enhance chosen environnement regarding alignement and numbering needs (alignat, gather and respective alignedat and gathered)
• Simplify the syntax by using new commands:
  • \pd[...]{...} for \frac{\partial ...}{\partial ...}
  • \lt-\rt in place of \right-\left for automatic sizing of delimiters
  • \vect ... for {\rm {\bf ...}}} and similar ones
• Align continuation lines to improve code readability

Location:

NEMO/trunk/doc/latex/NEMO/subfiles

Files:

• chap_model_basics.tex (modified) (50 diffs)
• chap_time_domain.tex (modified) (19 diffs)

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

@@ -2 +2 @@
 
 \begin{document}
+
 % ================================================================
-% Chapter 1 Model Basics
+% Chapter 1  Model Basics
 % ================================================================
-
 \chapter{Model Basics}
 \label{chap:PE}
-
 \minitoc
 
@@ -26 +25 @@
 \label{subsec:PE_Vector}
 
-
 The ocean is a fluid that can be described to a good approximation by the primitive equations,
 \ie the Navier-Stokes equations along with a nonlinear equation of state which
@@ -32 +30 @@
 plus the following additional assumptions made from scale considerations:
 
-\textit{(1) spherical earth approximation:} the geopotential surfaces are assumed to be spheres so that
-gravity (local vertical) is parallel to the earth's radius
-
-\textit{(2) thin-shell approximation:} the ocean depth is neglected compared to the earth's radius
-
-\textit{(3) turbulent closure hypothesis:} the turbulent fluxes
-(which represent the effect of small scale processes on the large-scale) are expressed in terms of
-large-scale features
-
-\textit{(4) Boussinesq hypothesis:} density variations are neglected except in their contribution to
-the buoyancy force
-
-\textit{(5) Hydrostatic hypothesis:} the vertical momentum equation is reduced to a balance between
-the vertical pressure gradient and the buoyancy force
-(this removes convective processes from the initial Navier-Stokes equations and so
-convective processes must be parameterized instead)
-
-\textit{(6) Incompressibility hypothesis:} the three dimensional divergence of the velocity vector is assumed to
-be zero.
+\begin{enumerate}
+\item
+  \textit{spherical earth approximation}: the geopotential surfaces are assumed to be spheres so that
+  gravity (local vertical) is parallel to the earth's radius
+\item
+  \textit{thin-shell approximation}: the ocean depth is neglected compared to the earth's radius
+\item
+  \textit{turbulent closure hypothesis}: the turbulent fluxes
+  (which represent the effect of small scale processes on the large-scale)
+  are expressed in terms of large-scale features
+\item
+  \textit{Boussinesq hypothesis}: density variations are neglected except in their contribution to
+  the buoyancy force
+  \begin{equation}
+    \label{eq:PE_eos}
+    \rho = \rho \ (T,S,p)
+  \end{equation}
+\item
+  \textit{Hydrostatic hypothesis}: the vertical momentum equation is reduced to a balance between
+  the vertical pressure gradient and the buoyancy force
+  (this removes convective processes from the initial Navier-Stokes equations and so
+  convective processes must be parameterized instead)
+  \begin{equation}
+    \label{eq:PE_hydrostatic}
+    \pd[p]{z} = - \rho \ g
+  \end{equation}
+\item
+  \textit{Incompressibility hypothesis}: the three dimensional divergence of the velocity vector $\vect U$
+  is assumed to be zero.
+  \begin{equation}
+    \label{eq:PE_continuity}
+    \nabla \cdot \vect U = 0
+  \end{equation}
+\end{enumerate}
 
 Because the gravitational force is so dominant in the equations of large-scale motions,
-it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to
-the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to
-\textbf{k}, \ie tangent to the geopotential surfaces.
-Let us define the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$ 
-(the subscript $h$ denotes the local horizontal vector, \ie over the (\textbf{i},\textbf{j}) plane), 
-$T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density.
-The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) vector system
-provides the following six equations
-(namely the momentum balance, the hydrostatic equilibrium, the incompressibility equation,
-the heat and salt conservation equations and an equation of state):
+it is useful to choose an orthogonal set of unit vectors $(i,j,k)$ linked to the earth such that
+$k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$,
+\ie tangent to the geopotential surfaces.
+Let us define the following variables: $\vect U$ the vector velocity, $\vect U = \vect U_h + w \, \vect k$
+(the subscript $h$ denotes the local horizontal vector, \ie over the $(i,j)$ plane), 
+$T$ the potential temperature, $S$ the salinity, $\rho$ the \textit{in situ} density.
+The vector invariant form of the primitive equations in the $(i,j,k)$ vector system provides
+the following equations:
 \begin{subequations}
   \label{eq:PE}
-  \begin{equation}
+  \begin{gather}
+    \intertext{$-$ the momentum balance}
     \label{eq:PE_dyn}
-    \frac{\partial {\rm {\bf U}}_h }{\partial t}=
-    -\left[    {\left( {\nabla \times {\rm {\bf U}}} \right)\times {\rm {\bf U}}
-        +\frac{1}{2}\nabla \left( {{\rm {\bf U}}^2} \right)}    \right]_h
-    -f\;{\rm {\bf k}}\times {\rm {\bf U}}_h
-    -\frac{1}{\rho_o }\nabla _h p + {\rm {\bf D}}^{\rm {\bf U}} + {\rm {\bf F}}^{\rm {\bf U}}
-  \end{equation}
-  \begin{equation}
-    \label{eq:PE_hydrostatic}
-    \frac{\partial p }{\partial z} = - \rho \ g
-  \end{equation}
-  \begin{equation}
-    \label{eq:PE_continuity}
-    \nabla \cdot {\bf U}=  0
-  \end{equation}
-  \begin{equation}
+    \pd[\vect U_h]{t} = - \lt[ (\nabla \times \vect U) \times \vect U + \frac{1}{2} \nabla \lt( \vect U^2 \rt) \rt]_h
+                        - f \; k \times \vect U_h - \frac{1}{\rho_o} \nabla_h p
+                        + \vect D^{\vect U} + \vect F^{\vect U} \\
+    \intertext{$-$ the heat and salt conservation equations}
     \label{eq:PE_tra_T}
-    \frac{\partial T}{\partial t} = - \nabla \cdot  \left( T \ \rm{\bf U} \right) + D^T + F^T
-  \end{equation}
-  \begin{equation}
+    \pd[T]{t} = - \nabla \cdot (T \ \vect U) + D^T + F^T \\
     \label{eq:PE_tra_S}
-    \frac{\partial S}{\partial t} = - \nabla \cdot  \left( S \ \rm{\bf U} \right) + D^S + F^S
-  \end{equation}
-  \begin{equation}
-    \label{eq:PE_eos}
-    \rho = \rho \left( T,S,p \right)
-  \end{equation}
+    \pd[S]{t} = - \nabla \cdot (S \ \vect U) + D^S + F^S
+  \end{gather}
 \end{subequations}
-where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions, $t$ is the time,
-$z$ is the vertical coordinate, $\rho $ is the \textit{in situ} density given by the equation of state
+where $\nabla$ is the generalised derivative vector operator in $(i,j,k)$ directions, $t$ is the time,
+$z$ is the vertical coordinate, $\rho$ is the \textit{in situ} density given by the equation of state
 (\autoref{eq:PE_eos}), $\rho_o$ is a reference density, $p$ the pressure,
-$f=2 \bf \Omega \cdot \bf k$ is the Coriolis acceleration
-(where $\bf \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration.
-${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum,
-temperature and salinity, and ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$ surface forcing terms.
+$f = 2 \vect \Omega \cdot k$ is the Coriolis acceleration
+(where $\vect \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration.
+$\vect D^{\vect U}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum,
+temperature and salinity, and $\vect F^{\vect U}$, $F^T$ and $F^S$ surface forcing terms.
 Their nature and formulation are discussed in \autoref{sec:PE_zdf_ldf} and \autoref{subsec:PE_boundary_condition}.
-
-
 
 % -------------------------------------------------------------------------------------------------------------
@@ -113 +106 @@
 An ocean is bounded by complex coastlines, bottom topography at its base and
 an air-sea or ice-sea interface at its top.
-These boundaries can be defined by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$,
+These boundaries can be defined by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,k,t)$,
 where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface.
-Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$, chosen as a mean sea surface
+Both $H$ and $\eta$ are usually referenced to a given surface, $z = 0$, chosen as a mean sea surface
 (\autoref{fig:ocean_bc}).
 Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with
@@ -127 +120 @@
 \begin{figure}[!ht]
   \begin{center}
-    \includegraphics[width=0.90\textwidth]{Fig_I_ocean_bc}
-    \caption{   \protect\label{fig:ocean_bc}
-      The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,t)$,
+    \includegraphics[]{Fig_I_ocean_bc}
+    \caption{
+      \protect\label{fig:ocean_bc}
+      The ocean is bounded by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,t)$,
       where $H$ is the depth of the sea floor and $\eta$ the height of the sea surface.
-      Both $H$ and $\eta$ are referenced to $z=0$.
+      Both $H$ and $\eta$ are referenced to $z = 0$.
     }
   \end{center}
 \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
 
 \begin{description}
@@ -161 +154 @@
   \begin{equation}
     \label{eq:PE_w_bbc}
-    w = -{\rm {\bf U}}_h \cdot  \nabla _h \left( H \right)
+    w = - \vect U_h \cdot \nabla_h (H)
   \end{equation}
   In addition, the ocean exchanges momentum with the earth through frictional processes.
@@ -167 +160 @@
   It must be parameterized in terms of turbulent fluxes using bottom and/or lateral boundary conditions.
   Its specification depends on the nature of the physical parameterisation used for
-  ${\rm {\bf D}}^{\rm {\bf U}}$ in \autoref{eq:PE_dyn}.
+  $\vect D^{\vect U}$ in \autoref{eq:PE_dyn}.
   It is discussed in \autoref{eq:PE_zdf}.% and Chap. III.6 to 9.
 \item[Atmosphere - ocean interface:]
@@ -174 +167 @@
   \[
     % \label{eq:PE_w_sbc}
-    w = \frac{\partial \eta }{\partial t}
-    + \left. {{\rm {\bf U}}_h } \right|_{z=\eta } \cdot  \nabla _h \left( \eta \right)
-    + P-E
+    w = \pd[\eta]{t} + \lt. \vect U_h \rt|_{z = \eta} \cdot \nabla_h (\eta) + P - E
   \]
   The dynamic boundary condition, neglecting the surface tension (which removes capillary waves from the system)
-  leads to the continuity of pressure across the interface $z=\eta$.
+  leads to the continuity of pressure across the interface $z = \eta$.
   The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat.
 \item[Sea ice - ocean interface:]
   the ocean and sea ice exchange heat, salt, fresh water and momentum.
   The sea surface temperature is constrained to be at the freezing point at the interface.
-  Sea ice salinity is very low ($\sim4-6 \,psu$) compared to those of the ocean ($\sim34 \,psu$).
+  Sea ice salinity is very low ($\sim4-6 \, psu$) compared to those of the ocean ($\sim34 \, psu$).
   The cycle of freezing/melting is associated with fresh water and salt fluxes that cannot be neglected.
 \end{description}
-
-
-%\newpage
 
 % ================================================================
 % The Horizontal Pressure Gradient
 % ================================================================
-\section{Horizontal pressure gradient }
+\section{Horizontal pressure gradient}
 \label{sec:PE_hor_pg}
 
@@ -204 +192 @@
 
 The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at
-a reference geopotential surface ($z=0$) and a hydrostatic pressure $p_h$ such that:
-$p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$.
+a reference geopotential surface ($z = 0$) and a hydrostatic pressure $p_h$ such that:
+$p(i,j,k,t) = p_s(i,j,t) + p_h(i,j,k,t)$.
 The latter is computed by integrating (\autoref{eq:PE_hydrostatic}),
 assuming that pressure in decibars can be approximated by depth in meters in (\autoref{eq:PE_eos}).
@@ -211 +199 @@
 \[
   % \label{eq:PE_pressure}
-  p_h \left( {i,j,z,t} \right)
-  = \int_{\varsigma =z}^{\varsigma =0} {g\;\rho \left( {T,S,\varsigma} \right)\;d\varsigma }
+  p_h (i,j,z,t) = \int_{\varsigma = z}^{\varsigma = 0} g \; \rho (T,S,\varsigma) \; d \varsigma
 \]
 Two strategies can be considered for the surface pressure term:
@@ -224 +211 @@
 The phase speed of such waves is high (some hundreds of metres per second) so that
 the time step would have to be very short if they were present in the model.
-The latter strategy filters out these waves since the rigid lid approximation implies $\eta=0$,
-\ie the sea surface is the surface $z=0$.
+The latter strategy filters out these waves since the rigid lid approximation implies $\eta = 0$,
+\ie the sea surface is the surface $z = 0$.
 This well known approximation increases the surface wave speed to infinity and
 modifies certain other longwave dynamics (\eg barotropic Rossby or planetary waves).
 The rigid-lid hypothesis is an obsolescent feature in modern OGCMs.
-It has been available until the release 3.1 of  \NEMO, and it has been removed in release 3.2 and followings.
+It has been available until the release 3.1 of \NEMO, and it has been removed in release 3.2 and followings.
 Only the free surface formulation is now described in the this document (see the next sub-section).
 
@@ -244 +231 @@
 \begin{equation}
   \label{eq:PE_ssh}
-  \frac{\partial \eta }{\partial t}=-D+P-E
-  \quad \text{where} \
-  D=\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right]
+  \pd[\eta]{t} = - D + P - E \quad \text{where} \quad D = \nabla \cdot \lt[ (H + \eta) \; \overline{U}_h \, \rt]
 \end{equation}
 and using (\autoref{eq:PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$.
@@ -256 +241 @@
 
 Two choices can be made regarding the implementation of the free surface in the model,
-depending on the physical processes of interest. 
+depending on the physical processes of interest.
 
 $\bullet$ If one is interested in EGWs, in particular the tides and their interaction with
 the baroclinic structure of the ocean (internal waves) possibly in shallow seas,
 then a non linear free surface is the most appropriate.
-This means that no approximation is made in (\autoref{eq:PE_ssh}) and that
+This means that no approximation is made in \autoref{eq:PE_ssh} and that
 the variation of the ocean volume is fully taken into account.
 Note that in order to study the fast time scales associated with EGWs it is necessary to
@@ -272 +257 @@
 not altering the slow barotropic Rossby waves.
 If further, an approximative conservation of heat and salt contents is sufficient for the problem solved,
-then it is sufficient to solve a linearized version of (\autoref{eq:PE_ssh}),
+then it is sufficient to solve a linearized version of \autoref{eq:PE_ssh},
 which still allows to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}.
 Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost.
@@ -286 +271 @@
 (see \autoref{subsec:DYN_spg_ts}).
 
-%\newpage
-
 % ================================================================
 % Curvilinear z-coordinate System
@@ -293 +276 @@
 \section{Curvilinear \textit{z-}coordinate system}
 \label{sec:PE_zco}
-
 
 % -------------------------------------------------------------------------------------------------------------
@@ -318 +300 @@
 The general case is detailed by \citet{Eiseman1980} in their survey of the conservation laws of fluid dynamics.
 
-Let (\textit{i},\textit{j},\textit{k}) be a set of orthogonal curvilinear coordinates on
+Let $(i,j,k)$ be a set of orthogonal curvilinear coordinates on
 the sphere associated with the positively oriented orthogonal set of unit vectors
-(\textbf{i},\textbf{j},\textbf{k}) linked to the earth such that
-\textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k},
+$(i,j,k)$ linked to the earth such that
+$k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$,
 \ie along geopotential surfaces (\autoref{fig:referential}).
 Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined by
 the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and
-the distance from the centre of the earth $a+z(k)$ where $a$ is the earth's radius and
+the distance from the centre of the earth $a + z(k)$ where $a$ is the earth's radius and
 $z$ the altitude above a reference sea level (\autoref{fig:referential}).
 The local deformation of the curvilinear coordinate system is given by $e_1$, $e_2$ and $e_3$,
@@ -332 +314 @@
   \label{eq:scale_factors}
   \begin{aligned}
-    e_1 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda}{\partial i}\cos \varphi } \right)^2
-        +\left( {\frac{\partial \varphi }{\partial i}} \right)^2} \right]^{1/2} \\
-    e_2 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda }{\partial j}\cos \varphi } \right)^2+
-        \left( {\frac{\partial \varphi }{\partial j}} \right)^2} \right]^{1/2} \\
-    e_3 &=\left( {\frac{\partial z}{\partial k}} \right) \\
+    e_1 &= (a + z) \lt[ \lt( \pd[\lambda]{i} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{i} \rt)^2 \rt]^{1/2} \\
+    e_2 &= (a + z) \lt[ \lt( \pd[\lambda]{j} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{j} \rt)^2 \rt]^{1/2} \\
+    e_3 &= \lt( \pd[z]{k} \rt)
   \end{aligned}
 \end{equation}
@@ -343 +323 @@
 \begin{figure}[!tb]
   \begin{center}
-    \includegraphics[width=0.60\textwidth]{Fig_I_earth_referential}
-    \caption{  \protect\label{fig:referential}
+    \includegraphics[]{Fig_I_earth_referential}
+    \caption{
+      \protect\label{fig:referential}
       the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear
-      coordinate system (\textbf{i},\textbf{j},\textbf{k}). }
+      coordinate system $(i,j,k)$.
+    }
   \end{center}
 \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-Since the ocean depth is far smaller than the earth's radius, $a+z$, can be replaced by $a$ in
+Since the ocean depth is far smaller than the earth's radius, $a + z$, can be replaced by $a$ in
 (\autoref{eq:scale_factors}) (thin-shell approximation).
 The resulting horizontal scale factors $e_1$, $e_2$  are independent of $k$ while
-the vertical scale factor is a single function of $k$ as \textbf{k} is parallel to \textbf{z}.
+the vertical scale factor is a single function of $k$ as $k$ is parallel to $z$.
 The scalar and vector operators that appear in the primitive equations
 (\autoref{eq:PE_dyn} to \autoref{eq:PE_eos}) can be written in the tensorial form,
@@ -360 +342 @@
 \begin{subequations}
   % \label{eq:PE_discrete_operators}
-  \begin{equation}
+  \begin{gather}
     \label{eq:PE_grad}
-    \nabla q=\frac{1}{e_1 }\frac{\partial q}{\partial i}\;{\rm {\bf
-        i}}+\frac{1}{e_2 }\frac{\partial q}{\partial j}\;{\rm {\bf j}}+\frac{1}{e_3
-    }\frac{\partial q}{\partial k}\;{\rm {\bf k}}    \\
-  \end{equation}
-  \begin{equation}
+    \nabla q =   \frac{1}{e_1} \pd[q]{i} \; \vect i
+               + \frac{1}{e_2} \pd[q]{j} \; \vect j
+               + \frac{1}{e_3} \pd[q]{k} \; \vect k \\
     \label{eq:PE_div}
-    \nabla \cdot {\rm {\bf A}}
-    = \frac{1}{e_1 \; e_2} \left[
-      \frac{\partial \left(e_2 \; a_1\right)}{\partial i }
-      +\frac{\partial \left(e_1 \; a_2\right)}{\partial j }       \right]
-    + \frac{1}{e_3} \left[ \frac{\partial a_3}{\partial k }   \right]
-  \end{equation}
-  \begin{equation}
+    \nabla \cdot \vect A =   \frac{1}{e_1 \; e_2} \lt[ \pd[(e_2 \; a_1)]{\partial i} + \pd[(e_1 \; a_2)]{j} \rt]
+                           + \frac{1}{e_3} \lt[ \pd[a_3]{k} \rt]
+  \end{gather}
+  \begin{multline}
     \label{eq:PE_curl}
-    \begin{split}
-      \nabla \times \vect{A} =
-      \left[ {\frac{1}{e_2 }\frac{\partial a_3}{\partial j}
-          -\frac{1}{e_3 }\frac{\partial a_2 }{\partial k}} \right] \; \vect{i}
-      &+\left[ {\frac{1}{e_3 }\frac{\partial a_1 }{\partial k}
-          -\frac{1}{e_1 }\frac{\partial a_3 }{\partial i}} \right] \; \vect{j}      \\
-      &+\frac{1}{e_1 e_2 } \left[ {\frac{\partial \left( {e_2 a_2 } \right)}{\partial i}
-          -\frac{\partial \left( {e_1 a_1 } \right)}{\partial j}} \right] \; \vect{k}
-    \end{split}
-  \end{equation}
-  \begin{equation}
+      \nabla \times \vect{A} =   \lt[ \frac{1}{e_2} \pd[a_3]{j} - \frac{1}{e_3} \pd[a_2]{k}   \rt] \vect i \\
+                               + \lt[ \frac{1}{e_3} \pd[a_1]{k} - \frac{1}{e_1} \pd[a_3]{i}   \rt] \vect j \\
+                               + \frac{1}{e_1 e_2} \lt[ \pd[(e_2 a_2)]{i} - \pd[(e_1 a_1)]{j} \rt] \vect k
+  \end{multline}
+  \begin{gather}
     \label{eq:PE_lap}
-    \Delta q = \nabla \cdot \left(  \nabla q \right)
-  \end{equation}
-  \begin{equation}
+    \Delta q = \nabla \cdot (\nabla q) \\
     \label{eq:PE_lap_vector}
-    \Delta {\rm {\bf A}} =
-    \nabla \left( \nabla \cdot {\rm {\bf A}} \right)
-    - \nabla \times \left(  \nabla \times {\rm {\bf A}} \right)
-  \end{equation}
+    \Delta \vect A = \nabla (\nabla \cdot \vect A) - \nabla \times (\nabla \times \vect A)
+  \end{gather}
 \end{subequations}
-where $q$ is a scalar quantity and ${\rm {\bf A}}=(a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinate system.
+where $q$ is a scalar quantity and $\vect A = (a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinates system.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -408 +375 @@
 it is necessary to compute the horizontal component of the non-linear and viscous terms of the equation using
 \autoref{eq:PE_grad}) to \autoref{eq:PE_lap_vector}.
-Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate system and
+Let us set $\vect U = (u,v,w) = \vect U_h + w \; \vect k $, the velocity in the $(i,j,k)$ coordinates system and
 define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, by:
-\begin{equation}
+\begin{gather}
   \label{eq:PE_curl_Uh}
-  \zeta =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,v}
-        \right)}{\partial i}-\frac{\partial \left( {e_1 \,u} \right)}{\partial j}}
-  \right]
-\end{equation}
-\begin{equation}
+  \zeta = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, v)]{i} - \pd[(e_1 \, u)]{j} \rt] \\
   \label{eq:PE_div_Uh}
-  \chi =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,u}
-        \right)}{\partial i}+\frac{\partial \left( {e_1 \,v} \right)}{\partial j}}
-  \right]
-\end{equation}
+  \chi  = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, u)]{i} + \pd[(e_1 \, v)]{j} \rt]
+\end{gather}
 
 Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that
 $e_3$  is a function of the single variable $k$,
-the nonlinear term of \autoref{eq:PE_dyn} can be transformed as follows:
-\begin{flalign*}
-  &\left[ {\left( {  \nabla \times {\rm {\bf U}}    } \right) \times {\rm {\bf U}}
-      +\frac{1}{2}   \nabla \left( {{\rm {\bf U}}^2} \right)}   \right]_h        &
-\end{flalign*}
-\begin{flalign*}
-  &\qquad=\left( {{
-        \begin{array}{*{20}c}
-          {\left[    {   \frac{1}{e_3} \frac{\partial u  }{\partial k}
-          -\frac{1}{e_1}   \frac{\partial w  }{\partial i} } \right] w - \zeta \; v }     \\
-          {\zeta \; u - \left[ {  \frac{1}{e_2} \frac{\partial w}{\partial j}
-          -\frac{1}{e_3} \frac{\partial v}{\partial k} } \right] \ w}    \\
-        \end{array}
-      }} \right)
-  +\frac{1}{2} \left( {{
-        \begin{array}{*{20}c}
-          { \frac{1}{e_1}  \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial i}}  \hfill    \\
-          { \frac{1}{e_2}  \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial j}}  \hfill    \\
-        \end{array}
-      }} \right)        &
-\end{flalign*}
-\begin{flalign*}
-  & \qquad =\left( {{
-        \begin{array}{*{20}c}
-          {-\zeta \; v} \hfill \\
-          { \zeta \; u} \hfill \\
-        \end{array}
-      }} \right)
-  +\frac{1}{2}\left( {{
-        \begin{array}{*{20}c}
-          {\frac{1}{e_1 }\frac{\partial \left( {u^2+v^2} \right)}{\partial i}} \hfill  \\
-          {\frac{1}{e_2 }\frac{\partial \left( {u^2+v^2} \right)}{\partial j}} \hfill  \\
-        \end{array}
-      }} \right)
-  +\frac{1}{e_3 }\left( {{
-        \begin{array}{*{20}c}
-          { w \; \frac{\partial u}{\partial k}}    \\
-          { w \; \frac{\partial v}{\partial k}}    \\
-        \end{array}
-      }} \right)
-  -\left( {{
-        \begin{array}{*{20}c}
-          {\frac{w}{e_1}\frac{\partial w}{\partial i} -\frac{1}{2e_1}\frac{\partial w^2}{\partial i}} \hfill \\
-          {\frac{w}{e_2}\frac{\partial w}{\partial j} -\frac{1}{2e_2}\frac{\partial w^2}{\partial j}} \hfill \\
-        \end{array}
-      }} \right)        &
-\end{flalign*}
-
+$NLT$ the nonlinear term of \autoref{eq:PE_dyn} can be transformed as follows:
+\begin{alignat*}{2}
+  &NLT &=   &\lt[ (\nabla \times {\vect U}) \times {\vect U} + \frac{1}{2} \nabla \lt( {\vect U}^2 \rt) \rt]_h \\
+  &    &=   &\lt(
+    \begin{array}{*{20}c}
+                \lt[ \frac{1}{e_3} \pd[u]{k} - \frac{1}{e_1} \pd[w]{i} \rt] w - \zeta \; v   \\
+                \zeta \; u - \lt[ \frac{1}{e_2} \pd[w]{j} - \frac{1}{e_3} \pd[v]{k} \rt] \ w
+    \end{array}
+                                                                                             \rt)
+          + \frac{1}{2} \lt(
+    \begin{array}{*{20}c}
+                             \frac{1}{e_1} \pd[(u^2 + v^2 + w^2)]{i} \\
+                             \frac{1}{e_2} \pd[(u^2 + v^2 + w^2)]{j}
+    \end{array}
+                                                                     \rt) \\
+  &    &=   &\lt(
+    \begin{array}{*{20}c}
+                  -\zeta \; v \\
+                   \zeta \; u
+    \end{array}
+                              \rt)
+          + \frac{1}{2} \lt(
+    \begin{array}{*{20}c}
+                             \frac{1}{e_1} \pd[(u^2 + v^2)]{i} \\
+                             \frac{1}{e_2} \pd[(u^2 + v^2)]{j}
+    \end{array}
+                                                               \rt) \\
+  &    &  &+ \frac{1}{e_3} \lt(
+    \begin{array}{*{20}c}
+                                w \; \pd[u]{k} \\
+                                w \; \pd[v]{k}
+    \end{array}
+                                               \rt)
+           - \lt(
+    \begin{array}{*{20}c}
+                  \frac{w}{e_1} \pd[w]{i} - \frac{1}{2 e_1} \pd[w^2]{i} \\
+                  \frac{w}{e_2} \pd[w]{j} - \frac{1}{2 e_2} \pd[w^2]{j}
+    \end{array}
+                                                                        \rt)
+\end{alignat*}
 The last term of the right hand side is obviously zero, and thus the nonlinear term of
 \autoref{eq:PE_dyn} is written in the $(i,j,k)$ coordinate system:
 \begin{equation}
   \label{eq:PE_vector_form}
-  \left[ {\left( {   \nabla \times {\rm {\bf U}}    } \right) \times {\rm {\bf U}}
-      +\frac{1}{2}   \nabla \left( {{\rm {\bf U}}^2} \right)}   \right]_h
-  =\zeta
-  \;{\rm {\bf k}}\times {\rm {\bf U}}_h +\frac{1}{2}\nabla _h \left( {{\rm
-        {\bf U}}_h^2 } \right)+\frac{1}{e_3 }w\frac{\partial {\rm {\bf U}}_h
-  }{\partial k}
+  NLT =   \zeta \; \vect k \times \vect U_h + \frac{1}{2} \nabla_h \lt( \vect U_h^2 \rt)
+        + \frac{1}{e_3} w \pd[\vect U_h]{k}
 \end{equation}
 
@@ -489 +438 @@
 \ie to write it as the divergence of fluxes.
 For example, the first component of \autoref{eq:PE_vector_form} (the $i$-component) is transformed as follows:
-\begin{flalign*}
-  &{
-    \begin{array}{*{20}l}
-      \left[ {\left( {\nabla \times \vect{U}} \right)\times \vect{U}
-      +\frac{1}{2}\nabla \left( {\vect{U}}^2 \right)} \right]_i   % \\
-  % \\
-      = - \zeta \;v
-      + \frac{1}{2\;e_1 } \frac{\partial \left( {u^2+v^2} \right)}{\partial i}
-      + \frac{1}{e_3}w \ \frac{\partial u}{\partial k}         \\ \\
-      \qquad =\frac{1}{e_1 \; e_2} \left(    -v\frac{\partial \left( {e_2 \,v} \right)}{\partial i}
-      +v\frac{\partial \left( {e_1 \,u} \right)}{\partial j}    \right)
-      +\frac{1}{e_1 e_2 }\left(  +e_2 \; u\frac{\partial u}{\partial i}
-      +e_2 \; v\frac{\partial v}{\partial i}              \right)
-      +\frac{1}{e_3}       \left(   w\;\frac{\partial u}{\partial k}       \right)   \\
+\begin{alignat*}{2}
+  &NLT_i &= &- \zeta \; v + \frac{1}{2 \; e_1} \pd[ (u^2 + v^2) ]{i} + \frac{1}{e_3} w \ \pd[u]{k} \\
+  &      &=  &\frac{1}{e_1 \; e_2} \lt( -v \pd[(e_2 \, v)]{i} + v \pd[(e_1 \, u)]{j} \rt)
+            + \frac{1}{e_1 e_2} \lt( e_2 \; u \pd[u]{i} + e_2 \; v \pd[v]{i} \rt) \\
+  &      & &+ \frac{1}{e_3} \lt( w \; \pd[u]{k} \rt) \\
+  &      &=  &\frac{1}{e_1 \; e_2} \lt[ - \lt( v^2 \pd[e_2]{i} + e_2 \, v \pd[v]{i} \rt)
+                                     + \lt( \pd[ \lt( e_1 \, u \, v \rt)]{j} -         e_1 \, u \pd[v]{j} \rt) \rt. \\
+  &      &                       &\lt. + \lt( \pd[ \lt( e_2 \, u \, u \rt)]{i} - u \pd[ \lt( e_2 u \rt)]{i} \rt)
+                                     + e_2 v \pd[v]{i}                                                         \rt] \\
+  &      & &+ \frac{1}{e_3} \lt( \pd[(w \, u)]{k} - u \pd[w]{k} \rt) \\
+  &      &=  &\frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, u \, v)]{j} \rt)
+            + \frac{1}{e_3} \pd[(w \, u)]{k} \\
+  &      & &+ \frac{1}{e_1 e_2} \lt[ - u \lt( \pd[(e_1 v)]{j} - v \, \pd[e_1]{j} \rt)
+                                  - u \pd[(e_2 u)]{i}                              \rt]
+            - \frac{1}{e_3} \pd[w]{k} u \\
+  &      & &+ \frac{1}{e_1 e_2} \lt( - v^2 \pd[e_2]{i} \rt) \\
+  &      &= &\nabla \cdot (\vect U \, u) - (\nabla \cdot \vect U) \ u
+            + \frac{1}{e_1 e_2} \lt( -v^2 \pd[e_2]{i} + u v \, \pd[e_1]{j} \rt) \\
+  \intertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it comes:}
+  &      &= &\, \nabla \cdot (\vect U \, u) + \frac{1}{e_1 e_2} \lt( v \; \pd[e_2]{i} - u \; \pd[e_1]{j} \rt) (-v)
+\end{alignat*}
+
+The flux form of the momentum advection term is therefore given by:
+\begin{equation}
+  \label{eq:PE_flux_form}
+  NLT =   \nabla \cdot \lt(
+    \begin{array}{*{20}c}
+                            \vect U \, u \\
+                            \vect U \, v
     \end{array}
-  }         &
-\end{flalign*}
-\begin{flalign*}
-  &{
-    \begin{array}{*{20}l}
-      \qquad =\frac{1}{e_1 \; e_2}  \left\{
-      -\left(         v^2  \frac{\partial e_2                                }{\partial i}
-      +e_2 \,v    \frac{\partial v                                   }{\partial i}     \right)
-      +\left(           \frac{\partial \left( {e_1 \,u\,v}  \right)}{\partial j}
-      -e_1 \,u    \frac{\partial v                                   }{\partial j}  \right)  \right. \\
-      \left.   \qquad \qquad \quad
-      +\left(           \frac{\partial \left( {e_2 u\,u}     \right)}{\partial i}
-      -u       \frac{\partial \left( {e_2 u}         \right)}{\partial i}  \right)
-      +e_2 v            \frac{\partial v                                    }{\partial i}
-      \right\}
-      +\frac{1}{e_3} \left(
-      \frac{\partial \left( {w\,u} \right)         }{\partial k}
-      -u       \frac{\partial w                    }{\partial k}  \right) \\
-    \end{array}
-  }      &
-\end{flalign*}
-\begin{flalign*}
-  &
-  {
-    \begin{array}{*{20}l}
-      \qquad =\frac{1}{e_1 \; e_2}  \left(
-      \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i}
-      +        \frac{\partial \left( {e_1 \,u\,v} \right)}{\partial j}  \right)
-      +\frac{1}{e_3 }      \frac{\partial \left( {w\,u       } \right)}{\partial k} \\
-      \qquad \qquad \quad
-      +\frac{1}{e_1 e_2 }     \left(
-      -u \left(   \frac{\partial \left( {e_1 v   } \right)}{\partial j}
-      -v\,\frac{\partial e_1 }{\partial j}             \right)
-      -u       \frac{\partial \left( {e_2 u   } \right)}{\partial i}
-      \right)
-      -\frac{1}{e_3 }      \frac{\partial w}{\partial k} u
-      +\frac{1}{e_1 e_2 }\left(  -v^2\frac{\partial e_2   }{\partial i}     \right)
-    \end{array}
-  }      &
-\end{flalign*}
-\begin{flalign*}
-  &{
-    \begin{array}{*{20}l}
-      \qquad = \nabla \cdot \left( {{\rm {\bf U}}\,u} \right)
-      -   \left( \nabla \cdot {\rm {\bf U}} \right) \ u
-      +\frac{1}{e_1 e_2 }\left(
-      -v^2     \frac{\partial e_2 }{\partial i}
-      +uv   \,    \frac{\partial e_1 }{\partial j}    \right) \\
-    \end{array}
-  }      &
-\end{flalign*}
-as $\nabla \cdot {\rm {\bf U}}\;=0$ (incompressibility) it comes:
-\begin{flalign*}
-  &{
-    \begin{array}{*{20}l}
-      \qquad = \nabla \cdot \left(  {{\rm {\bf U}}\,u}      \right)
-      +  \frac{1}{e_1 e_2 }   \left( v \; \frac{\partial e_2}{\partial i}
-      -u \; \frac{\partial e_1}{\partial j}  \right)  \left( -v \right)
-    \end{array}
-  }      &
-\end{flalign*}
-
-The flux form of the momentum advection term is therefore given by:
-\begin{multline}
-  \label{eq:PE_flux_form}
-  \left[
-    \left(  {\nabla \times {\rm {\bf U}}}    \right) \times {\rm {\bf U}}
-    +\frac{1}{2}  \nabla \left(  {{\rm {\bf U}}^2}    \right)
-  \right]_h \\
-  = \nabla \cdot  \left( {{
-        \begin{array}{*{20}c}
-          {\rm {\bf U}} \, u  \hfill \\
-          {\rm {\bf U}} \, v  \hfill \\
-        \end{array}
-      }}
-  \right)
-  +\frac{1}{e_1 e_2 }      \left(
-    v\frac{\partial e_2}{\partial i}
-    -u\frac{\partial e_1}{\partial j}
-  \right) {\rm {\bf k}} \times {\rm {\bf U}}_h
-\end{multline}
+                                         \rt)
+        + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \vect k \times \vect U_h
+\end{equation}
 
 The flux form has two terms,
 the first one is expressed as the divergence of momentum fluxes (hence the flux form name given to this formulation)
 and the second one is due to the curvilinear nature of the coordinate system used.
-The latter is called the \emph{metric} term and can be viewed as a modification of the Coriolis parameter: 
+The latter is called the \textit{metric} term and can be viewed as a modification of the Coriolis parameter: 
 \[
   % \label{eq:PE_cor+metric}
-  f \to f + \frac{1}{e_1\;e_2}  \left(  v \frac{\partial e_2}{\partial i}
-    -u \frac{\partial e_1}{\partial j}  \right)
+  f \to f + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt)
 \]
 
 Note that in the case of geographical coordinate,
-\ie when $(i,j) \to (\lambda ,\varphi )$ and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$,
-we recover the commonly used modification of the Coriolis parameter $f \to f+(u/a) \tan \varphi$.
+\ie when $(i,j) \to (\lambda,\varphi)$ and $(e_1,e_2) \to (a \, \cos \varphi,a)$,
+we recover the commonly used modification of the Coriolis parameter $f \to f + (u / a) \tan \varphi$.
 
 To sum up, the curvilinear $z$-coordinate equations solved by the ocean model can be written in
 the following tensorial formalism:
 
-\vspace{+10pt}
-$\bullet$ \textbf{Vector invariant form of the momentum equations} :
-
-\begin{subequations}
-  \label{eq:PE_dyn_vect}
-  \[
+\begin{itemize}
+\item
+  \textbf{Vector invariant form of the momentum equations}:
+  \begin{equation}
+    \label{eq:PE_dyn_vect}
+    \begin{split}
     % \label{eq:PE_dyn_vect_u}
-    \begin{split}
-      \frac{\partial u}{\partial t}
-      = +   \left( {\zeta +f} \right)\,v
-      -   \frac{1}{2\,e_1}           \frac{\partial}{\partial i} \left(  u^2+v^2   \right)
-      -   \frac{1}{e_3    }  w     \frac{\partial u}{\partial k}      &      \\
-      -   \frac{1}{e_1    }            \frac{\partial}{\partial i} \left( \frac{p_s+p_h }{\rho_o}    \right)
-      &+   D_u^{\vect{U}}  +   F_u^{\vect{U}}      \\ \\
-      \frac{\partial v}{\partial t} =
-      -   \left( {\zeta +f} \right)\,u
-      -   \frac{1}{2\,e_2 }        \frac{\partial }{\partial j}\left(  u^2+v^2  \right)
-      -   \frac{1}{e_3     }   w  \frac{\partial v}{\partial k}     &      \\
-      -   \frac{1}{e_2     }        \frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o}  \right)
-      &+  D_v^{\vect{U}}  +   F_v^{\vect{U}}
+      \pd[u]{t} = &+ (\zeta + f) \, v - \frac{1}{2 e_1} \pd[]{i} (u^2 + v^2)
+                   - \frac{1}{e_3} w \pd[u]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\
+                  &+ D_u^{\vect U} + F_u^{\vect U} \\
+      \pd[v]{t} = &- (\zeta + f) \, u - \frac{1}{2 e_2} \pd[]{j} (u^2 + v^2)
+                   - \frac{1}{e_3} w \pd[v]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 
+                  &+ D_v^{\vect U} + F_v^{\vect U}
     \end{split}
-  \]
-\end{subequations}
-
-
-\vspace{+10pt}
-$\bullet$ \textbf{flux form of the momentum equations} :
-\begin{subequations}
+  \end{equation}
+\item
+  \textbf{flux form of the momentum equations}:
   % \label{eq:PE_dyn_flux}
   \begin{multline*}
     % \label{eq:PE_dyn_flux_u}
-    \frac{\partial u}{\partial t}=
-    +   \left( { f + \frac{1}{e_1 \; e_2}
-        \left(     v \frac{\partial e_2}{\partial i}
-          -u \frac{\partial e_1}{\partial j}    \right)}    \right) \, v    \\
-    - \frac{1}{e_1 \; e_2}    \left(
-      \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i}
-      +        \frac{\partial \left( {e_1 \,v\,u} \right)}{\partial j}  \right)
-    - \frac{1}{e_3 }\frac{\partial \left( {         w\,u} \right)}{\partial k}    \\
-    -   \frac{1}{e_1 }\frac{\partial}{\partial i}\left( \frac{p_s+p_h }{\rho_o}   \right)
-    +   D_u^{\vect{U}} +   F_u^{\vect{U}}
+    \pd[u]{t} = + \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v \\
+                - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) \\
+                - \frac{1}{e_3} \pd[(w \, u)]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt)
+                + D_u^{\vect U} + F_u^{\vect U}
   \end{multline*}
   \begin{multline*}
     % \label{eq:PE_dyn_flux_v}
-    \frac{\partial v}{\partial t}=
-    -   \left( { f + \frac{1}{e_1 \; e_2}
-        \left(     v \frac{\partial e_2}{\partial i}
-          -u \frac{\partial e_1}{\partial j}    \right)}    \right) \, u   \\
-    \frac{1}{e_1 \; e_2}   \left(
-      \frac{\partial \left( {e_2 \,u\,v} \right)}{\partial i}
-      +        \frac{\partial \left( {e_1 \,v\,v} \right)}{\partial j}  \right)
-    - \frac{1}{e_3 } \frac{\partial \left( {        w\,v} \right)}{\partial k}    \\
-    -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o}    \right)
-    +  D_v^{\vect{U}} +  F_v^{\vect{U}}
+    \pd[v]{t} = - \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u \\
+                + \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) \\
+                - \frac{1}{e_3} \pd[(w \, v)]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt)
+                + D_v^{\vect U} + F_v^{\vect U}
   \end{multline*}
-\end{subequations}
-where $\zeta$, the relative vorticity, is given by \autoref{eq:PE_curl_Uh} and
-$p_s $, the surface pressure, is given by:
-\[
+  where $\zeta$, the relative vorticity, is given by \autoref{eq:PE_curl_Uh} and $p_s$, the surface pressure,
+  is given by:
+  \[
   % \label{eq:PE_spg}
-  p_s =  \rho \,g \,\eta
-\]
-with $\eta$ is solution of \autoref{eq:PE_ssh}.
-
-The vertical velocity and the hydrostatic pressure are diagnosed from the following equations:
-\[
+    p_s = \rho \,g \, \eta
+  \]
+  with $\eta$ is solution of \autoref{eq:PE_ssh}.
+
+  The vertical velocity and the hydrostatic pressure are diagnosed from the following equations:
+  \[
   % \label{eq:w_diag}
-  \frac{\partial w}{\partial k}=-\chi \;e_3
-\]
-\[
+    \pd[w]{k} = - \chi \; e_3 \qquad 
   % \label{eq:hp_diag}
-  \frac{\partial p_h }{\partial k}=-\rho \;g\;e_3
-\]
-where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:PE_div_Uh}.
-
-\vspace{+10pt}
-$\bullet$ \textit{tracer equations} :
-\[
-  % \label{eq:S}
-  \frac{\partial T}{\partial t} =
-  -\frac{1}{e_1 e_2 }\left[ {    \frac{\partial \left( {e_2 T\,u} \right)}{\partial i}
-      +\frac{\partial \left( {e_1 T\,v} \right)}{\partial j}} \right]
-  -\frac{1}{e_3 }\frac{\partial \left( {T\,w} \right)}{\partial k} + D^T + F^T
-\]
-\[
-  % \label{eq:T}
-  \frac{\partial S}{\partial t} =
-  -\frac{1}{e_1 e_2 }\left[     {\frac{\partial \left( {e_2 S\,u} \right)}{\partial i}
-      +\frac{\partial \left( {e_1 S\,v} \right)}{\partial j}} \right]
-  -\frac{1}{e_3 }\frac{\partial \left( {S\,w} \right)}{\partial k} + D^S + F^S
-\]
-\[
-  % \label{eq:rho}
-  \rho =\rho \left( {T,S,z(k)} \right)
-\]
-
-The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale parameterisation used.
+    \pd[p_h]{k} = - \rho \; g \; e_3
+  \]
+  where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:PE_div_Uh}.
+\item \textit{tracer equations}:
+  \[
+    %\label{eq:S}
+    \pd[T]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt]
+                - \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\
+    %\label{eq:T}
+    \pd[S]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt]
+                - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S
+    %\label{eq:rho}
+    \rho = \rho \big( T,S,z(k) \big)
+  \]
+\end{itemize}
+
+The expression of $\vect D^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale parameterisation used.
 It will be defined in \autoref{eq:PE_zdf}.
-The nature and formulation of ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms,
+The nature and formulation of $\vect F^{\vect U}$, $F^T$ and $F^S$, the surface forcing terms,
 are discussed in \autoref{chap:SBC}.
 
-
-\newpage 
+\newpage
 
 % ================================================================
@@ -717 +566 @@
 Second the ocean floor depends on the geographical position,
 varying from more than 6,000 meters in abyssal trenches to zero at the coast.
-Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 
+Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing.
 Therefore, in order to represent the ocean with respect to
 the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height
@@ -736 +585 @@
 \begin{equation}
   \label{eq:PE_s}
-  s=s(i,j,k,t)
+  s = s(i,j,k,t)
 \end{equation}
 with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$,
@@ -759 +608 @@
 
 %\gmcomment{
-
 %A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient...
-
 the generalized vertical coordinates used in ocean modelling are not orthogonal,
 which contrasts with many other applications in mathematical physics.
@@ -772 +619 @@
 The key motivation for maintaining the same horizontal velocity component is that
 the hydrostatic and geostrophic balances are dominant in the large-scale ocean.
-Use of an alternative quasi-horizontal velocity, for example one oriented parallel to the generalized surface,
+Use of an alternative quasi -horizontal velocity, for example one oriented parallel to the generalized surface,
 would lead to unacceptable numerical errors.
-Correspondingly, the vertical direction is anti-parallel to the gravitational force in
+Correspondingly, the vertical direction is anti -parallel to the gravitational force in
 all of the coordinate systems.
-We do not choose the alternative of a quasi-vertical direction oriented normal to
-the surface of a constant generalized vertical coordinate. 
+We do not choose the alternative of a quasi -vertical direction oriented normal to
+the surface of a constant generalized vertical coordinate.
 
 It is the method used to measure transport across the generalized vertical coordinate surfaces which differs between
@@ -786 +633 @@
 volume or mass conservation.
 In other models, such as isopycnal layered models, this transport is prescribed based on assumptions about
-the physical processes producing a flux across the layer interfaces. 
-
+the physical processes producing a flux across the layer interfaces.
 
 In this section we first establish the PE in the generalised vertical $s$-coordinate,
-then we discuss the particular cases available in \NEMO, namely $z$, \zstar, $s$, and \ztilde.  
+then we discuss the particular cases available in \NEMO, namely $z$, \zstar, $s$, and \ztilde.
 %}
 
@@ -796 +642 @@
 % The s-coordinate Formulation
 % -------------------------------------------------------------------------------------------------------------
-\subsection{\textit{S-}coordinate formulation}
-
-Starting from the set of equations established in \autoref{sec:PE_zco} for the special case $k=z$ and thus $e_3=1$,
-we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$,
-which includes $z$-, \zstar- and $\sigma-$coordinates as special cases
-($s=z$, $s=\zstar$, and $s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.).
+\subsection{\textit{S}-coordinate formulation}
+
+Starting from the set of equations established in \autoref{sec:PE_zco} for the special case $k = z$ and
+thus $e_3 = 1$, we introduce an arbitrary vertical coordinate $s = s(i,j,k,t)$,
+which includes $z$-, \zstar- and $\sigma$-coordinates as special cases
+($s = z$, $s = \zstar$, and $s = \sigma = z / H$ or $ = z / \lt( H + \eta \rt)$, resp.).
 A formal derivation of the transformed equations is given in \autoref{apdx:A}.
-Let us define the vertical scale factor by $e_3=\partial_s z$  ($e_3$ is now a function of $(i,j,k,t)$ ),
-and the slopes in the (\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by:
+Let us define the vertical scale factor by $e_3 = \partial_s z$  ($e_3$ is now a function of $(i,j,k,t)$ ),
+and the slopes in the $(i,j)$ directions between $s$- and $z$-surfaces by:
 \begin{equation}
   \label{eq:PE_sco_slope}
-  \sigma_1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s
-  \quad \text{, and } \quad
-  \sigma_2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s
+  \sigma_1 = \frac{1}{e_1} \; \lt. \pd[z]{i} \rt|_s \quad \text{and} \quad
+  \sigma_2 = \frac{1}{e_2} \; \lt. \pd[z]{j} \rt|_s
 \end{equation}
-We also introduce  $\omega $, a dia-surface velocity component, defined as the velocity 
+We also introduce $\omega$, a dia-surface velocity component, defined as the velocity 
 relative to the moving $s$-surfaces and normal to them:
 \[
   % \label{eq:PE_sco_w}
-  \omega  = w - e_3 \, \frac{\partial s}{\partial t} - \sigma_1 \,u - \sigma_2 \,v    \\
+  \omega = w - e_3 \, \pd[s]{t} - \sigma_1 \, u - \sigma_2 \, v
 \]
 
-The equations solved by the ocean model \autoref{eq:PE} in $s-$coordinate can be written as follows
+The equations solved by the ocean model \autoref{eq:PE} in $s$-coordinate can be written as follows
 (see \autoref{sec:A_momentum}):
 
- \vspace{0.5cm}
-$\bullet$ Vector invariant form of the momentum equation :
-\begin{multline*}
+\begin{itemize}
+\item \textbf{Vector invariant form of the momentum equation}:
+  \begin{multline*}
   % \label{eq:PE_sco_u_vector}
-  \frac{\partial  u   }{\partial t}=
-  +   \left( {\zeta +f} \right)\,v
-  -   \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left(  u^2+v^2   \right)
-  -   \frac{1}{e_3} \omega \frac{\partial u}{\partial k}       \\
-  -   \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho_o}    \right)
-  +  g\frac{\rho }{\rho_o}\sigma_1
-  +   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad
-\end{multline*}
-\begin{multline*}
+    \pd[u]{t} = + (\zeta + f) \, v - \frac{1}{2 \, e_1} \pd[]{i} (u^2 + v^2) - \frac{1}{e_3} \omega \pd[u]{k} \\
+                - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + g \frac{\rho}{\rho_o} \sigma_1
+                + D_u^{\vect U} + F_u^{\vect U}
+  \end{multline*}
+  \begin{multline*}
   % \label{eq:PE_sco_v_vector}
-  \frac{\partial v }{\partial t}=
-  -   \left( {\zeta +f} \right)\,u
-  -   \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left(  u^2+v^2  \right)
-  -   \frac{1}{e_3 } \omega \frac{\partial v}{\partial k}         \\
-  -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o}  \right)
-  +  g\frac{\rho }{\rho_o }\sigma_2
-  +  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad
-\end{multline*}
-
- \vspace{0.5cm}
-$\bullet$ Flux form of the momentum equation :
-\begin{multline*}
+    \pd[v]{t} = - (\zeta + f) \, u - \frac{1}{2 \, e_2} \pd[]{j}(u^2 + v^2) - \frac{1}{e_3} \omega \pd[v]{k} \\
+                - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + g \frac{\rho}{\rho_o} \sigma_2
+                + D_v^{\vect U} + F_v^{\vect U}
+  \end{multline*}
+\item \textbf{Flux form of the momentum equation}:
+  \begin{multline*}
   % \label{eq:PE_sco_u_flux}
-  \frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t}=
-  +   \left( { f + \frac{1}{e_1 \; e_2 }
-      \left(    v \frac{\partial e_2}{\partial i}
-        -u \frac{\partial e_1}{\partial j}   \right)}    \right) \, v    \\
-  - \frac{1}{e_1 \; e_2 \; e_3 }    \left(
-    \frac{\partial \left( {e_2 \, e_3 \, u\,u} \right)}{\partial i}
-    +       \frac{\partial \left( {e_1 \, e_3 \, v\,u} \right)}{\partial j}   \right)
-  - \frac{1}{e_3 }\frac{\partial \left( { \omega\,u} \right)}{\partial k}    \\
-  - \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho_o}    \right)
-  +  g\frac{\rho }{\rho_o}\sigma_1
-  +   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad
-\end{multline*}
-\begin{multline*}
+    \frac{1}{e_3} \pd[(e_3 \, u)]{t} = + \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v \\
+                                       - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[(e_2 \, e_3 \, u \, u)]{i} + \pd[(e_1 \, e_3 \, v \, u)]{j} \rt) \\
+                                       - \frac{1}{e_3} \pd[(\omega \, u)]{k}
+                                       - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt)
+                                       + g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U}
+  \end{multline*}
+  \begin{multline*}
   % \label{eq:PE_sco_v_flux}
-  \frac{1}{e_3} \frac{\partial \left(  e_3\,v  \right) }{\partial t}=
-  -   \left( { f + \frac{1}{e_1 \; e_2}
-      \left(    v \frac{\partial e_2}{\partial i}
-        -u \frac{\partial e_1}{\partial j}   \right)}    \right) \, u   \\
-  - \frac{1}{e_1 \; e_2 \; e_3 }    \left(
-    \frac{\partial \left( {e_2 \; e_3  \,u\,v} \right)}{\partial i}
-    +       \frac{\partial \left( {e_1 \; e_3  \,v\,v} \right)}{\partial j}   \right)
-  - \frac{1}{e_3 } \frac{\partial \left( { \omega\,v} \right)}{\partial k}    \\
-  -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o}  \right)
-  +  g\frac{\rho }{\rho_o }\sigma_2
-  +  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad
-\end{multline*}
-
-where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient,
-and the hydrostatic pressure have the same expressions as in $z$-coordinates although
-they do not represent exactly the same quantities.
-$\omega$ is provided by the continuity equation (see \autoref{apdx:A}):
-\[
+    \frac{1}{e_3} \pd[(e_3 \, v)]{t} = - \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u \\
+                                       - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[( e_2 \; e_3 \, u \, v)]{i} + \pd[(e_1 \; e_3 \, v \, v)]{j} \rt) \\
+                                       - \frac{1}{e_3} \pd[(\omega \, v)]{k}
+                                       - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt)
+                                       + g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U}
+  \end{multline*}
+  where the relative vorticity, $\zeta$, the surface pressure gradient,
+  and the hydrostatic pressure have the same expressions as in $z$-coordinates although
+  they do not represent exactly the same quantities.
+  $\omega$ is provided by the continuity equation (see \autoref{apdx:A}):
+  \[
   % \label{eq:PE_sco_continuity}
-  \frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0
-  \qquad \text{with }\;\;
-  \chi =\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3 \,u}
-        \right)}{\partial i}+\frac{\partial \left( {e_1 e_3 \,v} \right)}{\partial
-        j}} \right]
-\]
-
- \vspace{0.5cm}
-$\bullet$ tracer equations:
-\begin{multline*}
+    \pd[e_3]{t} + e_3 \; \chi + \pd[\omega]{s} = 0 \quad \text{with} \quad
+    \chi = \frac{1}{e_1 e_2 e_3} \lt( \pd[(e_2 e_3 \, u)]{i} + \pd[(e_1 e_3 \, v)]{j} \rt)
+  \]
+\item \textit{tracer equations}:
+  \begin{multline*}
   % \label{eq:PE_sco_t}
-  \frac{1}{e_3} \frac{\partial \left(  e_3\,T  \right) }{\partial t}=
-  -\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,T} \right)}{\partial i}
-      +\frac{\partial \left( {e_1 e_3\,v\,T} \right)}{\partial j}} \right]   \\
-  -\frac{1}{e_3 }\frac{\partial \left( {T\,\omega } \right)}{\partial k}   + D^T + F^S   \qquad
-\end{multline*}
-
-\begin{multline*}
+    \frac{1}{e_3} \pd[(e_3 \, T)]{t} = - \frac{1}{e_1 e_2 e_3} \lt(   \pd[(e_2 e_3 \, u \, T)]{i}
+                                                                    + \pd[(e_1 e_3 \, v \, T)]{j} \rt) \\
+                                       - \frac{1}{e_3} \pd[(T \, \omega)]{k} + D^T + F^S
+  \end{multline*}
+  \begin{multline}
   % \label{eq:PE_sco_s}
-  \frac{1}{e_3} \frac{\partial \left(  e_3\,S  \right) }{\partial t}=
-  -\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,S} \right)}{\partial i}
-      +\frac{\partial \left( {e_1 e_3\,v\,S} \right)}{\partial j}} \right]    \\
-  -\frac{1}{e_3 }\frac{\partial \left( {S\,\omega } \right)}{\partial k}     + D^S + F^S   \qquad
-\end{multline*}
-
+    \frac{1}{e_3} \pd[(e_3 \, S)]{t} = - \frac{1}{e_1 e_2 e_3} \lt(   \pd[(e_2 e_3 \, u \, S)]{i}
+                                                                    + \pd[(e_1 e_3 \, v \, S)]{j} \rt) \\
+                                       - \frac{1}{e_3} \pd[(S \, \omega)]{k} + D^S + F^S
+  \end{multline}
+\end{itemize}
 The equation of state has the same expression as in $z$-coordinate,
 and similar expressions are used for mixing and forcing terms.
 
 \gmcomment{
-\colorbox{yellow}{ to be updated $= = >$}
-Add a few works on z and zps and s and underlies the differences between all of them
-\colorbox{yellow}{ $< = =$ end update}  }
-
-
+  \colorbox{yellow}{ to be updated $= = >$}
+  Add a few works on z and zps and s and underlies the differences between all of them
+  \colorbox{yellow}{$< = =$ end update}
+}
 
 % -------------------------------------------------------------------------------------------------------------
 % Curvilinear \zstar-coordinate System
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Curvilinear \zstar--coordinate system}
+\subsection{Curvilinear \zstar-coordinate system}
 \label{subsec:PE_zco_star}
 
@@ -925 +738 @@
 \begin{figure}[!b]
   \begin{center}
-    \includegraphics[width=1.0\textwidth]{Fig_z_zstar}
-    \caption{  \protect\label{fig:z_zstar}
+    \includegraphics[]{Fig_z_zstar}
+    \caption{
+      \protect\label{fig:z_zstar}
       (a) $z$-coordinate in linear free-surface case ;
-      (b) $z-$coordinate in non-linear free surface case ;
-      (c) re-scaled height coordinate (become popular as the \zstar-coordinate
-      \citep{Adcroft_Campin_OM04} ).
+      (b) $z$-coordinate in non-linear free surface case ;
+      (c) re-scaled height coordinate
+      (become popular as the \zstar-coordinate \citep{Adcroft_Campin_OM04}).
     }
   \end{center}
@@ -936 +750 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-
 In that case, the free surface equation is nonlinear, and the variations of volume are fully taken into account.
-These coordinates systems is presented in a report \citep{Levier2007} available on the \NEMO web site. 
-
-%\gmcomment{
+These coordinates systems is presented in a report \citep{Levier2007} available on the \NEMO web site.
+
 The \zstar coordinate approach is an unapproximated, non-linear free surface implementation which allows one to
 deal with large amplitude free-surface variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}.
@@ -947 +759 @@
 as in the $z$-coordinate formulation, but is equally distributed over the full water column.
 Thus vertical levels naturally follow sea-surface variations, with a linear attenuation with depth,
-as illustrated by figure fig.1c.
-Note that with a flat bottom, such as in fig.1c, the bottom-following $z$ coordinate and \zstar are equivalent.
-The definition and modified oceanic equations for the rescaled vertical coordinate  \zstar,
+as illustrated by \autoref{fig:z_zstar}.
+Note that with a flat bottom, such as in \autoref{fig:z_zstar}, the bottom-following $z$ coordinate and \zstar are equivalent.
+The definition and modified oceanic equations for the rescaled vertical coordinate \zstar,
 including the treatment of fresh-water flux at the surface, are detailed in Adcroft and Campin (2004).
 The major points are summarized here.
-The position ( \zstar) and vertical discretization (\zstar) are expressed as:
+The position (\zstar) and vertical discretization (\zstar) are expressed as:
 \[
   % \label{eq:z-star}
-  H +  \zstar = (H + z) / r \quad \text{and} \ \delta \zstar = \delta z / r \quad \text{with} \ r = \frac{H+\eta} {H}
-\] 
+  H + \zstar = (H + z)  / r \quad \text{and}  \quad \delta \zstar
+              = \delta z / r \quad \text{with} \quad r
+              = \frac{H + \eta}{H}
+\]
 Since the vertical displacement of the free surface is incorporated in the vertical coordinate \zstar,
 the upper and lower boundaries are at fixed  \zstar position,
-$\zstar = 0$ and  $\zstar = -H$ respectively.
+$\zstar = 0$ and $\zstar = -H$ respectively.
 Also the divergence of the flow field is no longer zero as shown by the continuity equation:
-\[ 
-  \frac{\partial r}{\partial t} = \nabla_{\zstar} \cdot \left( r \; \rm{\bf U}_h \right)
-  \left( r \; w\textit{*} \right) = 0 
-\] 
-%}
-
+\[
+  \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) (r \; w *) = 0
+\]
 
 % from MOM4p1 documentation
-
 To overcome problems with vanishing surface and/or bottom cells, we consider the zstar coordinate 
 \[
   % \label{eq:PE_}
-  z^\star = H \left( \frac{z-\eta}{H+\eta} \right)
+  \zstar = H \lt( \frac{z - \eta}{H + \eta} \rt)
 \]
 
@@ -981 +791 @@
 and it has been recently promoted by Adcroft and Campin (2004) for global climate modelling.
 
-The surfaces of constant $z^\star$ are quasi-horizontal.
-Indeed, the $z^\star$ coordinate reduces to $z$ when $\eta$ is zero.
+The surfaces of constant \zstar are quasi -horizontal.
+Indeed, the \zstar coordinate reduces to $z$ when $\eta$ is zero.
 In general, when noting the large differences between
 undulations of the bottom topography versus undulations in the surface height,
-it is clear that surfaces constant $z^\star$ are very similar to the depth surfaces.
+it is clear that surfaces constant \zstar are very similar to the depth surfaces.
 These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to
 terrain following sigma models discussed in \autoref{subsec:PE_sco}.
-Additionally, since $z^\star$ when $\eta = 0$,
+Additionally, since \zstar when $\eta = 0$,
 no flow is spontaneously generated in an unforced ocean starting from rest, regardless the bottom topography.
 This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in the presence of
 nontrivial topographic variations can generate nontrivial spontaneous flow from a resting state,
 depending on the sophistication of the pressure gradient solver.
-The quasi-horizontal nature of the coordinate surfaces also facilitates the implementation of
-neutral physics parameterizations in $z^\star$ models using the same techniques as in $z$-models
+The quasi -horizontal nature of the coordinate surfaces also facilitates the implementation of
+neutral physics parameterizations in \zstar models using the same techniques as in $z$-models
 (see Chapters 13-16 of \cite{Griffies_Bk04}) for a discussion of neutral physics in $z$-models,
 as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO).
 
-The range over which $z^\star$ varies is time independent $-H \leq z^\star \leq 0$.
+The range over which \zstar varies is time independent $-H \leq \zstar \leq 0$.
 Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > ?H$.
-This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$. 
-
-Because $z^\star$ has a time independent range, all grid cells have static increments ds,
-and the sum of the ver tical increments yields the time independent ocean depth. %k ds = H. 
-The $z^\star$ coordinate is therefore invisible to undulations of the free surface,
+This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$.
+
+Because \zstar has a time independent range, all grid cells have static increments ds,
+and the sum of the ver tical increments yields the time independent ocean depth. %k ds = H.
+The \zstar coordinate is therefore invisible to undulations of the free surface,
 since it moves along with the free surface.
-This proper ty means that no spurious vertical transport is induced across surfaces of constant $z^\star$ by
+This proper ty means that no spurious vertical transport is induced across surfaces of constant \zstar by
 the motion of external gravity waves.
 Such spurious transpor t can be a problem in z-models, especially those with tidal forcing.
-Quite generally, the time independent range for the $z^\star$ coordinate is a very convenient property that
+Quite generally, the time independent range for the \zstar coordinate is a very convenient property that
 allows for a nearly arbitrary ver tical resolution even in the presence of large amplitude fluctuations of
 the surface height, again so long as $\eta > -H$.
-
 %end MOM doc %%%
 
-
-
-\newpage
+\newpage 
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1035 +842 @@
 For example, the topographic $\beta$-effect is usually larger than the planetary one along continental slopes.
 Topographic Rossby waves can be excited and can interact with the mean current.
-In the $z-$coordinate system presented in the previous section (\autoref{sec:PE_zco}),
-$z-$surfaces are geopotential surfaces.
+In the $z$-coordinate system presented in the previous section (\autoref{sec:PE_zco}),
+$z$-surfaces are geopotential surfaces.
 The bottom topography is discretised by steps.
 This often leads to a misrepresentation of a gradually sloping bottom and to
@@ -1043 +850 @@
 One solution to strongly reduce this error is to use a partial step representation of bottom topography instead of
 a full step one \cite{Pacanowski_Gnanadesikan_MWR98}.
-Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate).
+Another solution is to introduce a terrain-following coordinate system (hereafter $s$-coordinate).
 
 The $s$-coordinate avoids the discretisation error in the depth field since the layers of
@@ -1050 +857 @@
 which would be ignored in typical $z$-model applications with the largest grid spacing at greatest depths,
 can easily be represented (with relatively low vertical resolution).
-A terrain-following model (hereafter $s-$model) also facilitates the modelling of the boundary layer flows over
+A terrain-following model (hereafter $s$-model) also facilitates the modelling of the boundary layer flows over
 a large depth range, which in the framework of the $z$-model would require high vertical resolution over
 the whole depth range.
@@ -1063 +870 @@
 \begin{equation}
   \label{eq:PE_p_sco}
-  \left. {\nabla p} \right|_z =\left. {\nabla p} \right|_s -\frac{\partial
-    p}{\partial s}\left. {\nabla z} \right|_s 
+  \nabla p |_z = \nabla p |_s - \pd[p]{s} \nabla z |_s
 \end{equation}
 
 The second term in \autoref{eq:PE_p_sco} depends on the tilt of the coordinate surface and
 introduces a truncation error that is not present in a $z$-model.
-In the special case of a $\sigma-$coordinate (\ie depth-normalised coordinate system $\sigma = z/H$),
+In the special case of a $\sigma$-coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$),
 \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude of this truncation error.
 It depends on topographic slope, stratification, horizontal and vertical resolution, the equation of state,
@@ -1085 +891 @@
 (see \autoref{sec:DOM_zgr}.
 
-For numerical reasons a minimum of diffusion is required along the coordinate surfaces of any finite difference model.
+For numerical reasons a minimum of diffusion is required along the coordinate surfaces of
+any finite difference model.
 It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces.
 This is the case for a $z$-model as well as for a $s$-model.
-However, density varies more strongly on $s-$surfaces than on horizontal surfaces in regions of
+However, density varies more strongly on $s$-surfaces than on horizontal surfaces in regions of
 large topographic slopes, implying larger diapycnal diffusion in a $s$-model than in a $z$-model.
 Whereas such a diapycnal diffusion in a $z$-model tends to weaken horizontal density (pressure) gradients and thus
@@ -1101 +908 @@
 (see \autoref{subsec:PE_ldf}).
 Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large,
-strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}). 
-
-The $s-$coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two aspects from
+strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}).
+
+The $s$-coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two aspects from
 similar models:
 it allows a representation of bottom topography with mixed full or partial step-like/terrain following topography;
 It also offers a completely general transformation, $s=s(i,j,z)$ for the vertical coordinate.
 
-
-\newpage
-
 % -------------------------------------------------------------------------------------------------------------
 % Curvilinear z-tilde coordinate System
 % -------------------------------------------------------------------------------------------------------------
-\subsection{\texorpdfstring{Curvilinear \ztilde--coordinate}{}}
+\subsection{\texorpdfstring{Curvilinear \ztilde-coordinate}{}}
 \label{subsec:PE_zco_tilde}
 
-The \ztilde-coordinate has been developed by \citet{Leclair_Madec_OM11}.
+The \ztilde -coordinate has been developed by \citet{Leclair_Madec_OM11}.
 It is available in \NEMO since the version 3.4.
 Nevertheless, it is currently not robust enough to be used in all possible configurations.
 Its use is therefore not recommended.
-
 
 \newpage 
@@ -1134 +937 @@
 a few kilometres in the horizontal, a few meters in the vertical and a few minutes.
 They are usually solved at larger scales: the specified grid spacing and time step of the numerical model.
-The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations) must be represented entirely in terms of large-scale patterns to close the equations.
+The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations)
+must be represented entirely in terms of large-scale patterns to close the equations.
 These effects appear in the equations as the divergence of turbulent fluxes
 (\ie fluxes associated with the mean correlation of small scale perturbations).
@@ -1144 +948 @@
 
 The control exerted by gravity on the flow induces a strong anisotropy between the lateral and vertical motions.
-Therefore subgrid-scale physics \textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$  in
+Therefore subgrid-scale physics \textbf{D}$^{\vect U}$, $D^{S}$ and $D^{T}$  in
 \autoref{eq:PE_dyn}, \autoref{eq:PE_tra_T} and \autoref{eq:PE_tra_S} are divided into
-a lateral part  \textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and
-a vertical part  \textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$.
+a lateral part \textbf{D}$^{l \vect U}$, $D^{l S}$ and $D^{l T}$ and
+a vertical part \textbf{D}$^{v \vect U}$, $D^{v S}$ and $D^{v T}$.
 The formulation of these terms and their underlying physics are briefly discussed in the next two subsections.
 
@@ -1160 +964 @@
 Turbulent motions are thus never explicitly solved, even partially, but always parameterized.
 The vertical turbulent fluxes are assumed to depend linearly on the gradients of large-scale quantities
-(for example, the turbulent heat flux is given by $\overline{T'w'}=-A^{vT} \partial_z \overline T$,
-where $A^{vT}$ is an eddy coefficient).
+(for example, the turbulent heat flux is given by $\overline{T' w'} = -A^{v T} \partial_z \overline T$,
+where $A^{v T}$ is an eddy coefficient).
 This formulation is analogous to that of molecular diffusion and dissipation.
 This is quite clearly a necessary compromise: considering only the molecular viscosity acting on
@@ -1169 +973 @@
 \begin{equation}
   \label{eq:PE_zdf}
-  \begin{split}
-    {\vect{D}}^{v \vect{U}} &=\frac{\partial }{\partial z}\left( {A^{vm}\frac{\partial {\vect{U}}_h }{\partial z}} \right) \ , \\
-    D^{vT}                         &= \frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial T}{\partial z}} \right) \ ,
-    \quad
-    D^{vS}=\frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial S}{\partial z}} \right)
-  \end{split}
+  \begin{gathered}
+    \vect D^{v \vect U} = \pd[]{z} \lt( A^{vm} \pd[\vect U_h]{z} \rt) \ , \\
+          D^{vT}       = \pd[]{z} \lt( A^{vT} \pd[T]{z}         \rt) \quad \text{and} \quad
+          D^{vS}       = \pd[]{z} \lt( A^{vT} \pd[S]{z}         \rt)
+  \end{gathered}
 \end{equation}
 where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, respectively.
@@ -1232 +1035 @@
 one uses a advective scheme which is diffusive enough to maintain the model stability.
 It must be emphasised that then, all the sub-grid scale physics is included in the formulation of
-the advection scheme. 
+the advection scheme.
 
 All these parameterisations of subgrid scale physics have advantages and drawbacks.
@@ -1246 +1049 @@
 \begin{equation}
   \label{eq:PE_iso_tensor}
-  D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
-  \mbox{with}\quad \;\;\Re =\left( {{
-        \begin{array}{*{20}c}
-          1 \hfill & 0 \hfill & {-r_1 } \hfill \\
-          0 \hfill & 1 \hfill & {-r_2 } \hfill \\
-          {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\
-        \end{array}
-      }} \right)
+  D^{lT} = \nabla \vect . \lt( A^{lT} \; \Re \; \nabla T \rt) \quad \text{with} \quad
+  \Re =
+    \begin{pmatrix}
+      1    & 0    & -r_1          \\
+      0    & 1    & -r_2          \\
+      -r_1 & -r_2 & r_1^2 + r_2^2 \\
+    \end{pmatrix}
 \end{equation}
-where $r_1 \;\mbox{and}\;r_2 $ are the slopes between the surface along which the diffusive operator acts and
-the model level ($e. g.$ $z$- or $s$-surfaces).
+where $r_1$ and $r_2$ are the slopes between the surface along which the diffusive operator acts and
+the model level (\eg $z$- or $s$-surfaces).
 Note that the formulation \autoref{eq:PE_iso_tensor} is exact for
 the rotation between geopotential and $s$-surfaces,
@@ -1267 +1069 @@
 
 For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero.
-$\Re $ reduces to the identity in the horizontal direction, no rotation is applied. 
+$\Re$ reduces to the identity in the horizontal direction, no rotation is applied.
 
 For \textit{geopotential} diffusion,
@@ -1278 +1080 @@
 \begin{equation}
   \label{eq:PE_iso_slopes}
-  r_1 =\frac{e_3 }{e_1 }   \left( \pd[\rho]{i} \right) \left( \pd[\rho]{k} \right)^{-1} \, \quad
-  r_2 =\frac{e_3 }{e_2 }   \left( \pd[\rho]{j} \right) \left( \pd[\rho]{k} \right)^{-1} \,
+  r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad
+  r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1}
 \end{equation}
 while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$.
 
 \subsubsection{Eddy induced velocity}
+
 When the \textit{eddy induced velocity} parametrisation (eiv) \citep{Gent1990} is used,
 an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers:
 \[
   % \label{eq:PE_iso+eiv}
-  D^{lT}=\nabla \cdot \left( {A^{lT}\;\Re \;\nabla T} \right)
-  +\nabla \cdot \left( {{\vect{U}}^\ast \,T} \right)
+  D^{lT} = \nabla \cdot \lt( A^{lT} \; \Re \; \nabla T \rt) + \nabla \cdot \lt( \vect U^\ast \, T \rt)
 \]
-where ${\vect{U}}^\ast =\left( {u^\ast ,v^\ast ,w^\ast } \right)$ is a non-divergent,
+where $ \vect U^\ast = \lt( u^\ast,v^\ast,w^\ast \rt)$ is a non-divergent,
 eddy-induced transport velocity. This velocity field is defined by:
 \[
   % \label{eq:PE_eiv}
-  \begin{split}
-    u^\ast  &= +\frac{1}{e_3       }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_1 } \right] \\
-    v^\ast  &= +\frac{1}{e_3       }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_2 } \right] \\
-    w^\ast &=  -\frac{1}{e_1 e_2 }\left[
-      \frac{\partial }{\partial i}\left( {A^{eiv}\;e_2\,\tilde{r}_1 } \right)
-      +\frac{\partial }{\partial j}\left( {A^{eiv}\;e_1\,\tilde{r}_2 } \right)      \right]
-  \end{split}
+  u^\ast &=   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_1 \rt) \\
+  v^\ast &=   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_2 \rt) \\
+  w^\ast &= - \frac{1}{e_1 e_2} \lt[   \pd[]{i} \lt( A^{eiv} \; e_2 \, \tilde{r}_1 \rt)
+                                     + \pd[]{j} \lt( A^{eiv} \; e_1 \, \tilde{r}_2 \rt) \rt]
 \]
 where $A^{eiv}$ is the eddy induced velocity coefficient
 (or equivalently the isoneutral thickness diffusivity coefficient),
-and $\tilde{r}_1$ and $\tilde{r}_2$ are the slopes between isoneutral and \emph{geopotential} surfaces.
+and $\tilde r_1$ and $\tilde r_2$ are the slopes between isoneutral and \textit{geopotential} surfaces.
 Their values are thus independent of the vertical coordinate, but their expression depends on the coordinate: 
 \begin{align}
   \label{eq:PE_slopes_eiv}
   \tilde{r}_n =
-  \begin{cases}
-    r_n            &      \text{in $z$-coordinate}    \\
-    r_n + \sigma_n &      \text{in \zstar and $s$-coordinates}
-  \end{cases}
-                     \quad \text{where } n=1,2
+    \begin{cases}
+      r_n            & \text{in $z$-coordinate}                \\
+      r_n + \sigma_n & \text{in \zstar- and $s$-coordinates}
+    \end{cases}
+  \quad \text{where~} n = 1, 2
 \end{align}
 
 The normal component of the eddy induced velocity is zero at all the boundaries.
-This can be achieved in a model by tapering either the eddy coefficient or
-the slopes to zero in the vicinity of the boundaries.
+This can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in the vicinity of
+the boundaries.
 The latter strategy is used in \NEMO (cf. \autoref{chap:LDF}).
 
@@ -1327 +1126 @@
 \[
   % \label{eq:PE_bilapT}
-  D^{lT}= - \Delta \left( \;\Delta T \right)
-  \qquad \text{where} \;\; \Delta \bullet = \nabla \left( {\sqrt{B^{lT}\,}\;\Re \;\nabla \bullet} \right)
+  D^{lT}= - \Delta \; (\Delta T) \quad \text{where} \quad
+  \Delta \bullet = \nabla \lt( \sqrt{B^{lT}} \; \Re \; \nabla \bullet \rt)
 \]
 It is the Laplacian operator given by \autoref{eq:PE_iso_tensor} applied twice with
-the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 
-
+the harmonic eddy diffusion coefficient set to the square root of the biharmonic one.
 
 \subsubsection{Lateral Laplacian momentum diffusive operator}
@@ -1338 +1136 @@
 The Laplacian momentum diffusive operator along $z$- or $s$-surfaces is found by
 applying \autoref{eq:PE_lap_vector} to the horizontal velocity vector (see \autoref{apdx:B}):
-\[
+\begin{align*}
   % \label{eq:PE_lapU}
-  \begin{split}
-    {\rm {\bf D}}^{l{\rm {\bf U}}}
-    &= \quad \  \nabla _h \left( {A^{lm}\chi } \right)
-    \ - \ \nabla _h \times \left( {A^{lm}\,\zeta \;{\rm {\bf k}}} \right)     \\
-    &=   \left(
-      \begin{aligned}
-        \frac{1}{e_1      } \frac{\partial \left( A^{lm} \chi          \right)}{\partial i}
-        &-\frac{1}{e_2 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial j}  \\
-        \frac{1}{e_2      }\frac{\partial \left( {A^{lm} \chi         } \right)}{\partial j}
-        &+\frac{1}{e_1 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial i}
-      \end{aligned}
-    \right)
-  \end{split}
-\]
+  \vect D^{l \vect U} &=   \nabla_h        \big( A^{lm}    \chi             \big)
+                         - \nabla_h \times \big( A^{lm} \, \zeta \; \vect k \big) \\
+                      &= \lt(   \frac{1}{e_1}     \pd[ \lt( A^{lm}    \chi      \rt) ]{i} \rt.
+                              - \frac{1}{e_2 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{j}
+                                \frac{1}{e_2}     \pd[ \lt( A^{lm}    \chi      \rt) ]{j}
+                         \lt. + \frac{1}{e_1 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{i} \rt)
+\end{align*}
 
 Such a formulation ensures a complete separation between the vorticity and horizontal divergence fields
@@ -1359 +1150 @@
 Unfortunately, it is only available in \textit{iso-level} direction.
 When a rotation is required
-(\ie geopotential diffusion in $s-$coordinates or isoneutral diffusion in both $z$- and $s$-coordinates),
-the $u$ and $v-$fields are considered as independent scalar fields, so that the diffusive operator is given by:
-\[
+(\ie geopotential diffusion in $s$-coordinates or isoneutral diffusion in both $z$- and $s$-coordinates),
+the $u$ and $v$-fields are considered as independent scalar fields, so that the diffusive operator is given by:
+\begin{gather*}
   % \label{eq:PE_lapU_iso}
-  \begin{split}
-    D_u^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla u} \right) \\
-    D_v^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla v} \right)
-  \end{split}
-\]
+    D_u^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla u \rt) \\
+    D_v^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla v \rt)
+\end{gather*}
 where $\Re$ is given by \autoref{eq:PE_iso_tensor}.
 It is the same expression as those used for diffusive operator on tracers.
 It must be emphasised that such a formulation is only exact in a Cartesian coordinate system,
-\ie on a $f-$ or $\beta-$plane, not on the sphere.
+\ie on a $f$- or $\beta$-plane, not on the sphere.
 It is also a very good approximation in vicinity of the Equator in
 a geographical coordinate system \citep{Lengaigne_al_JGR03}.
@@ -1386 +1175 @@
 
 \end{document}
-

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

@@ -4 +4 @@
 
 % ================================================================
-% Chapter 2 Time Domain (step.F90)
+% Chapter 2 ——— Time Domain (step.F90)
 % ================================================================
-\chapter{Time Domain (STP) }
+\chapter{Time Domain (STP)}
 \label{chap:STP}
-
 \minitoc
 
@@ -35 +34 @@
 \begin{equation}
   \label{eq:STP}
-  x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \  \text{RHS}_x^{t-\rdt,\,t,\,t+\rdt}
+  x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t - \rdt, \, t, \, t + \rdt}
 \end{equation} 
 where $x$ stands for $u$, $v$, $T$ or $S$;
@@ -85 +84 @@
 \begin{equation}
   \label{eq:STP_asselin}
-  x_F^t  = x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right]
-\end{equation} 
+  x_F^t = x^t + \gamma \, \lt[ x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt]
+\end{equation}
 where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin coefficient.
 $\gamma$ is initialized as \np{rn\_atfp} (namelist parameter).
-Its default value is \np{rn\_atfp}\forcode{ = 10.e-3} (see \autoref{sec:STP_mLF}),
+Its default value is \np{rn\_atfp}~\forcode{= 10.e-3} (see \autoref{sec:STP_mLF}),
 causing only a weak dissipation of high frequency motions (\citep{Farge1987}).
 The addition of a time filter degrades the accuracy of the calculation from second to first order.
@@ -111 +110 @@
 (when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used :
 \[
-  % \label{eq:STP_euler}
-   x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \ {D_x}^{t-\rdt}
-\] 
+  %\label{eq:STP_euler}
+  x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ D_x^{t - \rdt}
+\]
 
 This is diffusive in time and conditionally stable.
@@ -119 +118 @@
 \begin{equation}
   \label{eq:STP_euler_stability}
-  A^h < \left\{
-    \begin{aligned}
-      &\frac{e^2}{  8 \; \rdt }  &&\quad \text{laplacian diffusion}  \\
-      &\frac{e^4}{64 \; \rdt }   &&\quad \text{bilaplacian diffusion}
-    \end{aligned}
-  \right.
+  A^h <
+  \begin{cases}
+    \frac{e^2}{ 8 \, \rdt} & \text{laplacian diffusion} \\
+    \frac{e^4}{64 \, \rdt} & \text{bilaplacian diffusion}
+  \end{cases}
 \end{equation}
 where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient.
@@ -134 +132 @@
 but usually the numerical stability condition imposes a strong constraint on the time step.
 Two solutions are available in \NEMO to overcome the stability constraint:
-$(a)$ a forward time differencing scheme using a time splitting technique (\np{ln\_zdfexp}\forcode{ = .true.}) or
-$(b)$ a backward (or implicit) time differencing scheme                   (\np{ln\_zdfexp}\forcode{ = .false.}).
-In $(a)$, the master time step $\Delta $t is cut into $N$ fractional time steps so that
+$(a)$ a forward time differencing scheme using a time splitting technique (\np{ln\_zdfexp}~\forcode{= .true.}) or
+$(b)$ a backward (or implicit) time differencing scheme                   (\np{ln\_zdfexp}~\forcode{= .false.}).
+In $(a)$, the master time step $\Delta$t is cut into $N$ fractional time steps so that
 the stability criterion is reduced by a factor of $N$.
 The computation is performed as follows:
-\[
+\begin{alignat*}{2}
   % \label{eq:STP_ts}
-  \begin{split}
-    & x_\ast ^{t-\rdt} = x^{t-\rdt}   \\
-    & x_\ast ^{t-\rdt+L\frac{2\rdt}{N}}=x_\ast ^{t-\rdt+\left( {L-1}
-      \right)\frac{2\rdt}{N}}+\frac{2\rdt}{N}\;\text{DF}^{t-\rdt+\left( {L-1} \right)\frac{2\rdt}{N}}
-    \quad \text{for $L=1$ to $N$}      \\
-    & x^{t+\rdt} = x_\ast^{t+\rdt}
-  \end{split}
-\]
+  &x_\ast^{t - \rdt}                      &= &x^{t - \rdt} \\
+  &x_\ast^{t - \rdt + L \frac{2 \rdt}{N}} &=   &x_\ast ^{t - \rdt + (L - 1) \frac{2 \rdt}{N}}
+                                             + \frac{2 \rdt}{N} \; DF^{t - \rdt + (L - 1) \frac{2 \rdt}{N}}
+  \quad \text{for $L = 1$ to $N$} \\
+  &x^{t + \rdt}                           &= &x_\ast^{t + \rdt}
+\end{alignat*}
 with DF a vertical diffusion term.
 The number of fractional time steps, $N$, is given by setting \np{nn\_zdfexp}, (namelist parameter).
@@ -154 +150 @@
 \begin{equation}
   \label{eq:STP_imp}
-  x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \  \text{RHS}_x^{t+\rdt}
-\end{equation} 
+  x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t + \rdt}
+\end{equation}
 
 %%gm
@@ -166 +162 @@
 \[
   % \label{eq:STP_imp_zdf}
-  \frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\rdt}\equiv \text{RHS}+\frac{1}{e_{3t} }\delta
-  _k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta_{k+1/2} \left[ {T^{t+1}} \right]}
-  \right]
+  \frac{T(k)^{t + 1} - T(k)^{t - 1}}{2 \; \rdt}
+  \equiv
+  \text{RHS} + \frac{1}{e_{3t}} \delta_k \lt[ \frac{A_w^{vT}}{e_{3w} } \delta_{k + 1/2} \lt[ T^{t + 1} \rt] \rt]
 \]
 where RHS is the right hand side of the equation except for the vertical diffusion term.
@@ -174 +170 @@
 \begin{equation}
   \label{eq:STP_imp_mat}
-  -c(k+1)\;T^{t+1}(k+1) + d(k)\;T^{t+1}(k) - \;c(k)\;T^{t+1}(k-1) \equiv b(k)
+  -c(k + 1) \; T^{t + 1}(k + 1) + d(k) \; T^{t + 1}(k) - \; c(k) \; T^{t + 1}(k - 1) \equiv b(k)
 \end{equation}
 where 
-\begin{align*}
+\begin{align*} 
   c(k) &= A_w^{vT} (k) \, / \, e_{3w} (k)     \\
-  d(k) &= e_{3t} (k)       \, / \, (2\rdt) + c_k + c_{k+1}    \\
-  b(k) &= e_{3t} (k) \; \left( T^{t-1}(k) \, / \, (2\rdt) + \text{RHS} \right)    
+  d(k) &= e_{3t}   (k)       \, / \, (2 \rdt) + c_k + c_{k + 1}    \\
+  b(k) &= e_{3t}   (k) \; \lt( T^{t - 1}(k) \, / \, (2 \rdt) + \text{RHS} \rt)
 \end{align*}
 
@@ -189 +185 @@
 (see for example \citet{Richtmyer1967}).
 
-
-
 % -------------------------------------------------------------------------------------------------------------
 %        Surface Pressure gradient
@@ -203 +197 @@
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[width=0.7\textwidth]{Fig_TimeStepping_flowchart}
+    \includegraphics[]{Fig_TimeStepping_flowchart}
     \caption{
       \protect\label{fig:TimeStep_flowchart}
       Sketch of the leapfrog time stepping sequence in \NEMO from \citet{Leclair_Madec_OM09}.
-      The use of a semi-implicit computation of the hydrostatic pressure gradient requires the tracer equation to
+      The use of a semi -implicit computation of the hydrostatic pressure gradient requires the tracer equation to
       be stepped forward prior to the momentum equation.
       The need for knowledge of the vertical scale factor (here denoted as $h$) requires the sea surface height and
@@ -225 +219 @@
 \label{sec:STP_mLF}
 
-Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the LF-RA scheme in order to ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter.
+Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the LF-RA scheme in order to
+ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter.
 The modifications affect both the forcing and filtering treatments in the LF-RA scheme.
 
 In a classical LF-RA environment, the forcing term is centred in time,
-\ie it is time-stepped over a $2\rdt$ period:
-$x^t  = x^t + 2\rdt Q^t $ where $Q$ is the forcing applied to $x$,
+\ie it is time-stepped over a $2 \rdt$ period:
+$x^t = x^t + 2 \rdt Q^t$ where $Q$ is the forcing applied to $x$,
 and the time filter is given by \autoref{eq:STP_asselin} so that $Q$ is redistributed over several time step.
 In the modified LF-RA environment, these two formulations have been replaced by:
-\begin{align}
-  x^{t+\rdt}  &= x^{t-\rdt} + \rdt \left( Q^{t-\rdt/2} + Q^{t+\rdt/2} \right)                   \label{eq:STP_forcing} \\
-  %
-  x_F^t  &= x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right]
-           - \gamma\,\rdt \, \left[ Q^{t+\rdt/2} -  Q^{t-\rdt/2} \right]                          \label{eq:STP_RA}
-\end{align}
+\begin{gather}
+  \label{eq:STP_forcing}
+  x^{t + \rdt} = x^{t - \rdt} + \rdt \lt( Q^{t - \rdt / 2} + Q^{t + \rdt / 2} \rt)  \\
+  \label{eq:STP_RA}
+  x_F^t       = x^t + \gamma \, \lt( x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt)
+                    - \gamma \, \rdt \, \lt( Q^{t + \rdt / 2} - Q^{t - \rdt / 2} \rt)
+\end{gather}
 The change in the forcing formulation given by \autoref{eq:STP_forcing} (see \autoref{fig:MLF_forcing})
 has a significant effect:
-the forcing term no longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}. 
+the forcing term no longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}.
 % forcing seen by the model....
 This property improves the LF-RA scheme in two respects.
@@ -250 +246 @@
 Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme,
 the modified formulation becomes conservative \citep{Leclair_Madec_OM09}.
-Second, the LF-RA becomes a truly quasi-second order scheme.
+Second, the LF-RA becomes a truly quasi -second order scheme.
 Indeed, \autoref{eq:STP_forcing} used in combination with a careful treatment of static instability
 (\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}),
 the two other main sources of time step divergence,
-allows a reduction by two orders of magnitude of the Asselin filter parameter. 
+allows a reduction by two orders of magnitude of the Asselin filter parameter.
 
 Note that the forcing is now provided at the middle of a time step:
-$Q^{t+\rdt/2}$ is the forcing applied over the $[t,t+\rdt]$ time interval.
+$Q^{t + \rdt / 2}$ is the forcing applied over the $[t,t + \rdt]$ time interval.
 This and the change in the time filter, \autoref{eq:STP_RA},
 allows an exact evaluation of the contribution due to the forcing term between any two time steps,
@@ -265 +261 @@
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[width=0.90\textwidth]{Fig_MLF_forcing}
+    \includegraphics[]{Fig_MLF_forcing}
     \caption{
       \protect\label{fig:MLF_forcing}
@@ -271 +267 @@
       (top) ''Traditional'' formulation:
       the forcing is defined at the same time as the variable to which it is applied
-      (integer value of the time step index) and it is applied over a $2\rdt$ period.
+      (integer value of the time step index) and it is applied over a $2 \rdt$ period.
       (bottom)  modified formulation:
       the forcing is defined in the middle of the time (integer and a half value of the time step index) and
-      the mean of two successive forcing values ($n-1/2$, $n+1/2$) is applied over a $2\rdt$ period.
+      the mean of two successive forcing values ($n - 1 / 2$, $n + 1 / 2$) is applied over a $2 \rdt$ period.
     }
   \end{center}
@@ -297 +293 @@
 \]
 This is done simply by keeping the leapfrog environment (\ie the \autoref{eq:STP} three level time stepping) but
-setting all $x^0$ (\textit{before}) and $x^{1}$ (\textit{now}) fields equal at the first time step and
+setting all $x^0$ (\textit{before}) and $x^1$ (\textit{now}) fields equal at the first time step and
 using half the value of $\rdt$.
 
@@ -305 +301 @@
 running the model for $2N$ time steps in one go,
 or by performing two consecutive experiments of $N$ steps with a restart.
-This requires saving two time levels and many auxiliary data in the restart files in machine precision. 
-
-Note that when a semi-implicit scheme is used to evaluate the hydrostatic pressure gradient
+This requires saving two time levels and many auxiliary data in the restart files in machine precision.
+
+Note that when a semi -implicit scheme is used to evaluate the hydrostatic pressure gradient
 (see \autoref{subsec:DYN_hpg_imp}), an extra three-dimensional field has to
 be added to the restart file to ensure an exact restartability.
-This is done optionally via the \np{nn\_dynhpg\_rst} namelist parameter,
+This is done optionally via the  \np{nn\_dynhpg\_rst} namelist parameter,
 so that the size of the restart file can be reduced when restartability is not a key issue
 (operational oceanography or in ensemble simulations for seasonal forecasting).
 
 Note the size of the time step used, $\rdt$, is also saved in the restart file.
-When restarting, if the the time step has been changed, a restart using an Euler time stepping scheme is imposed. 
-Options are defined through the \ngn{namrun} namelist variables.
+When restarting, if the the time step has been changed, a restart using an Euler time stepping scheme is imposed.
+Options are defined through the  \ngn{namrun} namelist variables.
 %%%
 \gmcomment{
@@ -360 +356 @@
 
 \begin{flalign*}
-  &\frac{\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}}{2\rdt}
-  \equiv \text{RHS}+ \delta_k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k+1/2} \left[ {T^{t+1}} \right]}
-  \right]      \\
-  &\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}
-  \equiv {2\rdt} \ \text{RHS}+ {2\rdt} \ \delta_k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k+1/2} \left[ {T^{t+1}} \right]}
-  \right]      \\
-  &\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}
+  &\frac{\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}}{2\rdt}
+  \equiv \text{RHS}+ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]}
+  \rt]      \\
+  &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}
+  \equiv {2\rdt} \ \text{RHS}+ {2\rdt} \ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]}
+  \rt]      \\
+  &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}
   \equiv 2\rdt \ \text{RHS}
-  + 2\rdt \ \left\{ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} [ T_{k+1}^{t+1} - T_k      ^{t+1} ]
-    - \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} [ T_k       ^{t+1} - T_{k-1}^{t+1} ]  \right\}     \\
+  + 2\rdt \ \lt\{ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} [ T_{k +1}^{t+1} - T_k      ^{t+1} ]
+    - \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} [ T_k       ^{t+1} - T_{k -1}^{t+1} ]  \rt\}     \\
   &\\
-  &\left( e_{3t}\,T \right)_k^{t+1}
-  -  {2\rdt} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}                  T_{k+1}^{t+1}
-  + {2\rdt} \ \left\{  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}
-    +  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}     \right\}   T_{k    }^{t+1}
-  -  {2\rdt} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                  T_{k-1}^{t+1}      \\
-  &\equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}    \\
-  % 
+  &\lt( e_{3t}\,T \rt)_k^{t+1}
+  -  {2\rdt} \           \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}                  T_{k +1}^{t+1}
+  + {2\rdt} \ \lt\{  \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}
+    +  \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}     \rt\}   T_{k    }^{t+1}
+  -  {2\rdt} \           \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}                  T_{k -1}^{t+1}      \\
+  &\equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}    \\
+  %
 \end{flalign*}
-
 \begin{flalign*}
   \allowdisplaybreaks
   \intertext{ Tracer case }
   %
-  &  \qquad \qquad  \quad   -  {2\rdt}                  \ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}
-  \qquad \qquad \qquad  \qquad  T_{k+1}^{t+1}   \\
-  &+ {2\rdt} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}
-  +   \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\
-  & \qquad \qquad  \qquad \qquad \qquad \quad \ \ -  {2\rdt} \                          \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                          \quad \ \ T_{k-1}^{t+1}
-  \ \equiv \ \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}  \\
+  &  \qquad \qquad  \quad   -  {2\rdt}                  \ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}
+  \qquad \qquad \qquad  \qquad  T_{k +1}^{t+1}   \\
+  &+ {2\rdt} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}
+  +   \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\
+  & \qquad \qquad  \qquad \qquad \qquad \quad \ \ -  {2\rdt} \                          \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}                          \quad \ \ T_{k -1}^{t+1}
+  \ \equiv \ \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}  \\
   %
 \end{flalign*}
@@ -396 +391 @@
   \intertext{ Tracer content case }
   %
-  & -  {2\rdt} \              & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_{k+1}^{t+1}}  && \  \left( e_{3t}\,T \right)_{k+1}^{t+1}   &\\
-  & + {2\rdt} \ \left[ 1  \right.+ & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_k^{t+1}}
-  + & \frac{(A_w^{vt})_{k -1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_k^{t+1}}  \left.  \right]  & \left( e_{3t}\,T \right)_{k   }^{t+1}  &\\
-  & -  {2\rdt} \               & \frac{(A_w^{vt})_{k -1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_{k-1}^{t+1}}     &\  \left( e_{3t}\,T \right)_{k-1}^{t+1}
-  \equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}  &
+  & -  {2\rdt} \              & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_{k +1}^{t+1}}  && \  \lt( e_{3t}\,T \rt)_{k +1}^{t+1}   &\\
+  & + {2\rdt} \ \lt[ 1  \rt.+ & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_k^{t+1}}
+  + & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_k^{t+1}}  \lt.  \rt]  & \lt( e_{3t}\,T \rt)_{k   }^{t+1}  &\\
+  & -  {2\rdt} \               & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_{k -1}^{t+1}}     &\  \lt( e_{3t}\,T \rt)_{k -1}^{t+1}
+  \equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}  &
 \end{flalign*}
 
 %%
 }
-%%
+
 \biblio