Index: branches/2017/dev_merge_2017/DOC/tex_main/NEMO_manual.sty
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_main/NEMO_manual.sty	(revision 9413)
+++ branches/2017/dev_merge_2017/DOC/tex_main/NEMO_manual.sty	(revision 9414)
@@ -21,8 +21,7 @@
 \usepackage{subfiles}	%% subdocs
 
-
 %% Extensions in bundle package
 
-\usepackage{amssymb, graphicx, longtable, makeidx, xspace}
+\usepackage{amssymb, graphicx, makeidx, tabularx, xspace}
 
 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/annex_A.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/annex_A.tex	(revision 9413)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/annex_A.tex	(revision 9414)
@@ -16,5 +16,5 @@
 % ================================================================
 \section{Chain rule for $s-$coordinates}
-\label{sec:A_continuity}
+\label{sec:A_chain}
 
 In order to establish the set of Primitive Equation in curvilinear $s$-coordinates
@@ -382,5 +382,5 @@
 Applying similar manipulation to the second component and replacing 
 $\sigma _1$ and $\sigma _2$ by their expression \autoref{apdx:A_s_slope}, it comes:
-\begin{equation} \label{apdx:A_grad_p}
+\begin{equation} \label{apdx:A_grad_p_1}
 \begin{split}
  -\frac{1}{\rho _o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z
@@ -394,5 +394,5 @@
 \end{equation}
 
-An additional term appears in (\autoref{apdx:A_grad_p}) which accounts for the 
+An additional term appears in (\autoref{apdx:A_grad_p_1}) which accounts for the 
 tilt of $s-$surfaces with respect to geopotential $z-$surfaces.
 
@@ -416,7 +416,7 @@
 \end{equation*}
 
-Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p} and using the definition of 
+Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and using the definition of 
 the density anomaly it comes the expression in two parts:
-\begin{equation} \label{apdx:A_grad_p}
+\begin{equation} \label{apdx:A_grad_p_2}
 \begin{split}
  -\frac{1}{\rho _o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z
@@ -430,5 +430,5 @@
 \end{equation}
 This formulation of the pressure gradient is characterised by the appearance of a term depending on the 
-the sea surface height only (last term on the right hand side of expression \autoref{apdx:A_grad_p}).
+the sea surface height only (last term on the right hand side of expression \autoref{apdx:A_grad_p_2}).
 This term will be loosely termed \textit{surface pressure gradient}
 whereas the first term will be termed the 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/annex_C.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/annex_C.tex	(revision 9413)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/annex_C.tex	(revision 9414)
@@ -122,9 +122,9 @@
 
 advection term (vector invariant form):
-\begin{equation} \label{eq:E_tot_vect_vor}
+\begin{equation} \label{eq:E_tot_vect_vor_1}
 \int\limits_D  \zeta \; \left( \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0   \\
 \end{equation}
 %
-\begin{equation} \label{eq:E_tot_vect_adv}
+\begin{equation} \label{eq:E_tot_vect_adv_1}
    \int\limits_D  \textbf{U}_h \cdot \nabla_h \left( \frac{{\textbf{U}_h}^2}{2} \right)     dv 
 + \int\limits_D  \textbf{U}_h \cdot \nabla_z \textbf{U}_h  \;dv    
@@ -151,5 +151,5 @@
 
 pressure gradient:
-\begin{equation} \label{eq:E_tot_pg}
+\begin{equation} \label{eq:E_tot_pg_1}
    - \int\limits_D  \left. \nabla p \right|_z \cdot \textbf{U}_h \;dv 
 = - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
@@ -172,13 +172,13 @@
 Vector invariant form:
 \begin{subequations} \label{eq:E_tot_vect}
-\begin{equation} \label{eq:E_tot_vect_vor}
+\begin{equation} \label{eq:E_tot_vect_vor_2}
 \int\limits_D   \textbf{U}_h \cdot \text{VOR} \;dv   = 0   \\
 \end{equation}
-\begin{equation} \label{eq:E_tot_vect_adv}
+\begin{equation} \label{eq:E_tot_vect_adv_2}
    \int\limits_D  \textbf{U}_h \cdot \text{KEG}  \;dv 
 + \int\limits_D  \textbf{U}_h \cdot \text{ZAD}  \;dv    
 -  \int\limits_D { \frac{{\textbf{U}_h}^2}{2} \frac{1}{e_3} \partial_t e_3 \;dv }   = 0   \\
 \end{equation}
-\begin{equation} \label{eq:E_tot_pg}
+\begin{equation} \label{eq:E_tot_pg_2}
    - \int\limits_D  \textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv 
 = - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
@@ -189,12 +189,12 @@
 Flux form:
 \begin{subequations} \label{eq:E_tot_flux}
-\begin{equation} \label{eq:E_tot_flux_metric}
+\begin{equation} \label{eq:E_tot_flux_metric_2}
 \int\limits_D  \textbf{U}_h \cdot \text {COR} \;  dv   = 0   \\
 \end{equation}
-\begin{equation} \label{eq:E_tot_flux_adv}
+\begin{equation} \label{eq:E_tot_flux_adv_2}
    \int\limits_D \textbf{U}_h \cdot \text{ADV}   \;dv 
 +   \frac{1}{2} \int\limits_D {  {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t e_3  \;dv } =\;0  \\
 \end{equation}
-\begin{equation} \label{eq:E_tot_pg}
+\begin{equation} \label{eq:E_tot_pg_3}
    - \int\limits_D  \textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv 
 = - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
@@ -207,6 +207,6 @@
 
 
-\autoref{eq:E_tot_pg} is the balance between the conversion KE to PE and PE to KE. 
-Indeed the left hand side of \autoref{eq:E_tot_pg} can be transformed as follows:
+\autoref{eq:E_tot_pg_3} is the balance between the conversion KE to PE and PE to KE. 
+Indeed the left hand side of \autoref{eq:E_tot_pg_3} can be transformed as follows:
 \begin{flalign*}
 \partial_t  \left( \int\limits_D { \rho \, g \, z  \;dv} \right) 
@@ -223,5 +223,5 @@
 the vertical velocity referenced to the fixe $z$-coordinate system (see \autoref{apdx:A_w_s}). 
  
-The left hand side of \autoref{eq:E_tot_pg} can be transformed as follows:
+The left hand side of \autoref{eq:E_tot_pg_3} can be transformed as follows:
 \begin{flalign*}
 - \int\limits_D  \left. \nabla p \right|_z & \cdot \textbf{U}_h \;dv  
@@ -325,5 +325,5 @@
 % ================================================================
 \section{Discrete total energy conservation: vector invariant form}
-\label{sec:C.1}
+\label{sec:C.2}
 
 % -------------------------------------------------------------------------------------------------------------
@@ -331,5 +331,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Total energy conservation}
-\label{subsec:C_KE+PE}
+\label{subsec:C_KE+PE_vect}
 
 The discrete form of the total energy conservation, \autoref{eq:Tot_Energy}, is given by:
@@ -401,8 +401,8 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsubsection{Vorticity term with EEN scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})}
-\label{subsec:C_vorEEN} 
+\label{subsec:C_vorEEN_vect} 
 
 With the EEN scheme, the vorticity terms are represented as: 
-\begin{equation} \label{eq:dynvor_een}
+\begin{equation} \tag{\ref{eq:dynvor_een}}
 \left\{ {    \begin{aligned}
  +q\,e_3 \, v 	&\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}} 
@@ -415,5 +415,5 @@
 $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$,
 and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by: 
-\begin{equation} \label{eq:Q_triads}
+\begin{equation} \tag{\ref{eq:Q_triads}}
 _i^j \mathbb{Q}^{i_p}_{j_p}
 = \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right)
@@ -580,5 +580,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Pressure gradient term}
-\label{subsec:C.1.4}
+\label{subsec:C.2.6}
 
 \gmcomment{
@@ -733,5 +733,5 @@
 % ================================================================
 \section{Discrete total energy conservation: flux form}
-\label{sec:C.1}
+\label{sec:C.3}
 
 % -------------------------------------------------------------------------------------------------------------
@@ -739,5 +739,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Total energy conservation}
-\label{subsec:C_KE+PE}
+\label{subsec:C_KE+PE_flux}
 
 The discrete form of the total energy conservation, \autoref{eq:Tot_Energy}, is given by:
@@ -763,5 +763,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Coriolis and advection terms: flux form}
-\label{subsec:C.1.3}
+\label{subsec:C.3.2}
 
 % -------------------------------------------------------------------------------------------------------------
@@ -769,5 +769,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsubsection{Coriolis plus ``metric'' term}
-\label{subsec:C.1.3.1} 
+\label{subsec:C.3.3} 
 
 In flux from the vorticity term reduces to a Coriolis term in which the Coriolis 
@@ -789,5 +789,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsubsection{Flux form advection}
-\label{subsec:C.1.3.2} 
+\label{subsec:C.3.4} 
 
 The flux form operator of the momentum advection is evaluated using a 
@@ -877,5 +877,5 @@
 % ================================================================
 \section{Discrete enstrophy conservation}
-\label{sec:C.1}
+\label{sec:C.4}
 
 
@@ -887,5 +887,5 @@
 
 In the ENS scheme, the vorticity term is descretized as follows:
-\begin{equation} \label{eq:dynvor_ens}
+\begin{equation} \tag{\ref{eq:dynvor_ens}}
 \left\{   \begin{aligned}
 +\frac{1}{e_{1u}} & \overline{q}^{\,i}  & {\overline{ \overline{\left( e_{1v}\,e_{3v}\;  v \right) } } }^{\,i, j+1/2}    \\
@@ -947,5 +947,5 @@
 
 With the EEN scheme, the vorticity terms are represented as: 
-\begin{equation} \label{eq:dynvor_een}
+\begin{equation} \tag{\ref{eq:dynvor_een}}
 \left\{ {    \begin{aligned}
  +q\,e_3 \, v 	&\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}} 
@@ -958,5 +958,5 @@
 $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$,
 and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by: 
-\begin{equation} \label{eq:Q_triads}
+\begin{equation} \tag{\ref{eq:Q_triads}}
 _i^j \mathbb{Q}^{i_p}_{j_p}
 = \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right)
@@ -1017,5 +1017,5 @@
 % ================================================================
 \section{Conservation properties on tracers}
-\label{sec:C.2}
+\label{sec:C.5}
 
 
@@ -1033,5 +1033,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Advection term}
-\label{subsec:C.2.1}
+\label{subsec:C.5.1}
 
 conservation of a tracer, $T$:
@@ -1123,5 +1123,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Conservation of potential vorticity}
-\label{subsec:C.3.1}
+\label{subsec:C.6.1}
 
 The lateral momentum diffusion term conserves the potential vorticity :
@@ -1157,5 +1157,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Dissipation of horizontal kinetic energy}
-\label{subsec:C.3.2}
+\label{subsec:C.6.2}
 
 The lateral momentum diffusion term dissipates the horizontal kinetic energy:
@@ -1209,5 +1209,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Dissipation of enstrophy}
-\label{subsec:C.3.3}
+\label{subsec:C.6.3}
 
 The lateral momentum diffusion term dissipates the enstrophy when the eddy 
@@ -1234,5 +1234,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Conservation of horizontal divergence}
-\label{subsec:C.3.4}
+\label{subsec:C.6.4}
 
 When the horizontal divergence of the horizontal diffusion of momentum 
@@ -1263,5 +1263,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Dissipation of horizontal divergence variance}
-\label{subsec:C.3.5}
+\label{subsec:C.6.5}
 
 \begin{flalign*}
@@ -1289,5 +1289,5 @@
 % ================================================================
 \section{Conservation properties on vertical momentum physics}
-\label{sec:C_4}
+\label{sec:C.7}
 
 As for the lateral momentum physics, the continuous form of the vertical diffusion 
@@ -1461,5 +1461,5 @@
 % ================================================================
 \section{Conservation properties on tracer physics}
-\label{sec:C.5}
+\label{sec:C.8}
 
 The numerical schemes used for tracer subgridscale physics are written such 
@@ -1473,5 +1473,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Conservation of tracers}
-\label{subsec:C.5.1}
+\label{subsec:C.8.1}
 
 constraint of conservation of tracers:
@@ -1507,5 +1507,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Dissipation of tracer variance}
-\label{subsec:C.5.2}
+\label{subsec:C.8.2}
 
 constraint on the dissipation of tracer variance:
Index: branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex	(revision 9413)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex	(revision 9414)
@@ -61,6 +61,6 @@
 The iso-neutral second order tracer diffusive operator for small
 angles between iso-neutral surfaces and geopotentials is given by
-\autoref{eq:PE_iso_tensor}:
-\begin{subequations} \label{eq:PE_iso_tensor}
+\autoref{eq:iso_tensor_1}:
+\begin{subequations} \label{eq:iso_tensor_1}
   \begin{equation}
     D^{lT}=-\Div\vect{f}^{lT}\equiv
@@ -73,5 +73,5 @@
   \end{equation}
   \begin{equation}
-    \label{eq:PE_iso_tensor:c}
+    \label{eq:iso_tensor_2}
     \mbox{with}\quad \;\;\Re =
     \begin{pmatrix}
@@ -118,5 +118,5 @@
 
 The off-diagonal terms of the small angle diffusion tensor
-\autoref{eq:PE_iso_tensor}, \autoref{eq:PE_iso_tensor:c} produce skew-fluxes along the
+\autoref{eq:iso_tensor_1}, \autoref{eq:iso_tensor_2} produce skew-fluxes along the
 $i$- and $j$-directions resulting from the vertical tracer gradient:
 \begin{align}
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_ASM.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_ASM.tex	(revision 9413)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_ASM.tex	(revision 9414)
@@ -90,5 +90,5 @@
 % Divergence damping description %%%
 \section{Divergence damping initialisation}
-\label{sec:ASM_details}
+\label{sec:ASM_div_dmp}
 
 The velocity increments may be initialized by the iterative application of 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIA.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIA.tex	(revision 9413)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIA.tex	(revision 9414)
@@ -1617,10 +1617,9 @@
 A non-Boussinesq fluid conserves mass. It satisfies the following relations:
 
-\begin{equation} \label{eq:MV_nBq}
-	\begin{split}
+\[	\begin{split}
 		\mathcal{M} &=  \mathcal{V}  \;\bar{\rho} \\
 		\mathcal{V} &=  \mathcal{A}  \;\bar{\eta}
-	\end{split}
-\end{equation}
+	\end{split}  \label{eq:MV_nBq}
+\]
 
 Temporal changes in total mass is obtained from the density conservation equation :
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex	(revision 9413)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex	(revision 9414)
@@ -557,5 +557,5 @@
 Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level functions for 
 (a) T-point depth and (b) the associated scale factor as computed 
-from \autoref{eq:DOM_zgr_ana} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate.}
+from \autoref{eq:DOM_zgr_ana_1} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -579,5 +579,5 @@
 near the ocean surface. The following function is proposed as a standard for a 
 $z$-coordinate (with either full or partial steps): 
-\begin{equation} \label{eq:DOM_zgr_ana}
+\begin{equation} \label{eq:DOM_zgr_ana_1}
 \begin{split}
  z_0 (k) 	&= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\ 
@@ -592,5 +592,5 @@
 If the ice shelf cavities are opened (\np{ln\_isfcav}\forcode{ = .true.}), the definition of $z_0$ is the same. 
 However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to:
-\begin{equation} \label{eq:DOM_zgr_ana}
+\begin{equation} \label{eq:DOM_zgr_ana_2}
 \begin{split}
  e_3^T(k) &= z_W (k+1) - z_W (k)   \\
@@ -616,5 +616,5 @@
 With the choice of the stretching $h_{cr} =3$ and the number of levels 
 \jp{jpk}=$31$, the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in 
-\autoref{eq:DOM_zgr_ana} have been determined such that \autoref{eq:DOM_zgr_coef} is 
+\autoref{eq:DOM_zgr_ana_2} have been determined such that \autoref{eq:DOM_zgr_coef} is 
 satisfied, through an optimisation procedure using a bisection method. For the first 
 standard ORCA2 vertical grid this led to the following values: $h_{sur} =4762.96$, 
@@ -677,5 +677,5 @@
 \caption{ \protect\label{tab:orca_zgr}   
 Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed 
-from \autoref{eq:DOM_zgr_ana} using the coefficients given in \autoref{eq:DOM_zgr_coef}}
+from \autoref{eq:DOM_zgr_ana_2} using the coefficients given in \autoref{eq:DOM_zgr_coef}}
 \end{table}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -761,5 +761,5 @@
 \begin{equation}
   s = -\frac{k}{n-1} \quad \text{ and } \quad 0 \leq k \leq n-1
-  \label{eq:s}
+  \label{eq:DOM_s}
 \end{equation}
 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_DYN.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_DYN.tex	(revision 9413)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_DYN.tex	(revision 9414)
@@ -414,5 +414,5 @@
 Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) 
 schemes can be used to compute the product of the Coriolis parameter and the 
-vorticity. However, the energy-conserving scheme  (\autoref{eq:dynvor_een}) has 
+vorticity. However, the energy-conserving scheme (\autoref{eq:dynvor_een}) has 
 exclusively been used to date. This term is evaluated using a leapfrog scheme, 
 $i.e.$ the velocity is centred in time (\textit{now} velocity).
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_model_basics.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_model_basics.tex	(revision 9413)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_model_basics.tex	(revision 9414)
@@ -648,8 +648,8 @@
 
 In fact one is totally free to choose any space and time vertical coordinate by introducing an arbitrary vertical coordinate :
-\begin{equation} \label{eq:s}
+\begin{equation} \label{eq:PE_s}
 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$, when $i$, $j$ and $t$ are held fixed. \autoref{eq:s} is a transformation from the $(i,j,k,t)$ coordinate system with independent variables into the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through \autoref{eq:s}.
+with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$, when $i$, $j$ and $t$ are held fixed. \autoref{eq:PE_s} is a transformation from the $(i,j,k,t)$ coordinate system with independent variables into the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through \autoref{eq:PE_s}.
 This so-called \textit{generalised vertical coordinate} \citep{Kasahara_MWR74} is in fact an Arbitrary Lagrangian--Eulerian (ALE) coordinate. Indeed, choosing an expression for $s$ is an arbitrary choice that determines which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and which part will be used to move them (Lagrangian part).
 The coordinate is also sometime referenced as an adaptive coordinate \citep{Hofmeister_al_OM09}, since the coordinate system is adapted in the course of the simulation. Its most often used implementation is via an ALE algorithm, in which a pure lagrangian step is followed by regridding and remapping steps, the later step implicitly embedding the vertical advection \citep{Hirt_al_JCP74, Chassignet_al_JPO03, White_al_JCP09}. Here we follow the \citep{Kasahara_MWR74} strategy : a regridding step (an update of the vertical coordinate) followed by an eulerian step with an explicit computation of vertical advection relative to the moving s-surfaces.
@@ -715,5 +715,5 @@
  \vspace{0.5cm}
 $\bullet$ Vector invariant form of the momentum equation :
-\begin{multline} \label{eq:PE_sco_u}
+\begin{multline} \label{eq:PE_sco_u_vector}
 \frac{\partial  u   }{\partial t}=
 	+   \left( {\zeta +f} \right)\,v                                    
@@ -724,5 +724,5 @@
 	+   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad
 \end{multline}
-\begin{multline} \label{eq:PE_sco_v}
+\begin{multline} \label{eq:PE_sco_v_vector}
 \frac{\partial v }{\partial t}=
 	-   \left( {\zeta +f} \right)\,u   
@@ -735,6 +735,6 @@
 
  \vspace{0.5cm}
-$\bullet$ Vector invariant form of the momentum equation :
-\begin{multline} \label{eq:PE_sco_u}
+$\bullet$ 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 }
@@ -749,5 +749,5 @@
 	+   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad
 \end{multline}
-\begin{multline} \label{eq:PE_sco_v}
+\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}
@@ -1138,5 +1138,5 @@
 rotation between geopotential and $s$-surfaces, while it is only an approximation 
 for the rotation between isoneutral and $z$- or $s$-surfaces. Indeed, in the latter 
-case, two assumptions are made to simplify  \autoref{eq:PE_iso_tensor} \citep{Cox1987}. 
+case, two assumptions are made to simplify \autoref{eq:PE_iso_tensor} \citep{Cox1987}. 
 First, the horizontal contribution of the dianeutral mixing is neglected since the ratio 
 between iso and dia-neutral diffusive coefficients is known to be several orders of 
