Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex	(revision 9393)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex	(revision 9407)
@@ -5,5 +5,5 @@
 % ================================================================
 \chapter{Space Domain (DOM)}
-\label{DOM}
+\label{chap:DOM}
 \minitoc
 
@@ -20,6 +20,6 @@
 $\ $\newline    % force a new line
 
-Having defined the continuous equations in Chap.~\ref{PE} and chosen a time 
-discretization Chap.~\ref{STP}, we need to choose a discretization on a grid, 
+Having defined the continuous equations in \autoref{chap:PE} and chosen a time 
+discretization \autoref{chap:STP}, we need to choose a discretization on a grid, 
 and numerical algorithms. In the present chapter, we provide a general description 
 of the staggered grid used in \NEMO, and other information relevant to the main 
@@ -32,5 +32,5 @@
 % ================================================================
 \section{Fundamentals of the discretisation}
-\label{DOM_basics}
+\label{sec:DOM_basics}
 
 % -------------------------------------------------------------------------------------------------------------
@@ -38,10 +38,10 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Arrangement of variables}
-\label{DOM_cell}
+\label{subsec:DOM_cell}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!tb]    \begin{center}
 \includegraphics[width=0.90\textwidth]{Fig_cell}
-\caption{ \protect\label{Fig_cell}    
+\caption{ \protect\label{fig:cell}    
 Arrangement of variables. $t$ indicates scalar points where temperature, 
 salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$) 
@@ -56,5 +56,5 @@
 space directions. The arrangement of variables is the same in all directions. 
 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector 
-points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}). 
+points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:cell}). 
 This is the generalisation to three dimensions of the well-known ``C'' grid in 
 Arakawa's classification \citep{Mesinger_Arakawa_Bk76}. The relative and 
@@ -66,8 +66,8 @@
 by the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$. 
 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as 
-indicated on Table \ref{Tab_cell}. In all the following, subscripts $u$, $v$, $w$, 
+indicated on \autoref{tab:cell}. In all the following, subscripts $u$, $v$, $w$, 
 $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale 
 factors are defined. Each scale factor is defined as the local analytical value 
-provided by \eqref{Eq_scale_factors}. As a result, the mesh on which partial 
+provided by \autoref{eq:scale_factors}. As a result, the mesh on which partial 
 derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$, and 
 $\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity. 
@@ -78,7 +78,7 @@
 from their analytical expression. This preserves the symmetry of the discrete set 
 of equations and therefore satisfies many of the continuous properties (see 
-Appendix~\ref{Apdx_C}). A similar, related remark can be made about the domain 
+\autoref{apdx:C}). A similar, related remark can be made about the domain 
 size: when needed, an area, volume, or the total ocean depth must be evaluated 
-as the sum of the relevant scale factors (see \eqref{DOM_bar}) in the next section). 
+as the sum of the relevant scale factors (see \autoref{eq:DOM_bar}) in the next section). 
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -95,9 +95,9 @@
 fw	& $i+1/2$ 	& $j+1/2$	& $k+1/2$ 	\\ \hline
 \end{tabular}
-\caption{ \protect\label{Tab_cell}
+\caption{ \protect\label{tab:cell}
 Location of grid-points as a function of integer or integer and a half value of the column, 
 line or level. This indexing is only used for the writing of the semi-discrete equation. 
 In the code, the indexing uses integer values only and has a reverse direction 
-in the vertical (see \S\ref{DOM_Num_Index})}
+in the vertical (see \autoref{subsec:DOM_Num_Index})}
 \end{center}
 \end{table}
@@ -108,9 +108,9 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Discrete operators}
-\label{DOM_operators}
+\label{subsec:DOM_operators}
 
 Given the values of a variable $q$ at adjacent points, the differencing and 
 averaging operators at the midpoint between them are:
-\begin{subequations} \label{Eq_di_mi}
+\begin{subequations} \label{eq:di_mi}
 \begin{align}
  \delta _i [q]       &=  \  \    q(i+1/2)  - q(i-1/2)		\\
@@ -120,14 +120,14 @@
 
 Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and 
-$k+1/2$. Following \eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a 
+$k+1/2$. Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a 
 variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$- 
 and $w$-points while its Laplacien is defined at $t$-point. These operators have 
 the following discrete forms in the curvilinear $s$-coordinate system:
-\begin{equation} \label{Eq_DOM_grad}
+\begin{equation} \label{eq:DOM_grad}
 \nabla q\equiv 	\frac{1}{e_{1u} } \delta _{i+1/2 } [q] \;\,\mathbf{i}
 		+	\frac{1}{e_{2v} } \delta _{j+1/2 } [q] \;\,\mathbf{j}
 		+	\frac{1}{e_{3w}} \delta _{k+1/2} [q] \;\,\mathbf{k}
 \end{equation}
-\begin{multline} \label{Eq_DOM_lap}
+\begin{multline} \label{eq:DOM_lap}
 \Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
        \;\left(          \delta_i  \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right]
@@ -136,8 +136,8 @@
 \end{multline}
 
-Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ 
+Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ 
 defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, 
 and $f$-points, and its divergence defined at $t$-points:
-\begin{eqnarray}  \label{Eq_DOM_curl}
+\begin{eqnarray}  \label{eq:DOM_curl}
  \nabla \times {\rm{\bf A}}\equiv &
       \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right)  &\ \mathbf{i} \\ 
@@ -145,5 +145,5 @@
  +& \frac{1}{e_{1f} \,e_{2f}    } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2  \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right)  &\ \mathbf{k}
  \end{eqnarray}
-\begin{eqnarray} \label{Eq_DOM_div}
+\begin{eqnarray} \label{eq:DOM_div}
 \nabla \cdot \rm{\bf A} \equiv 
     \frac{1}{e_{1t}\,e_{2t}\,e_{3t}} \left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right]
@@ -153,5 +153,5 @@
 The vertical average over the whole water column denoted by an overbar becomes 
 for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area):
-\begin{equation} \label{DOM_bar}
+\begin{equation} \label{eq:DOM_bar}
 \bar q 	=         \frac{1}{H}    \int_{k^b}^{k^o} {q\;e_{3q} \,dk} 
 		\equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} }
@@ -163,8 +163,8 @@
 
 In continuous form, the following properties are satisfied:
-\begin{equation} \label{Eq_DOM_curl_grad}
+\begin{equation} \label{eq:DOM_curl_grad}
 \nabla \times \nabla q ={\rm {\bf {0}}}
 \end{equation}
-\begin{equation} \label{Eq_DOM_div_curl}
+\begin{equation} \label{eq:DOM_div_curl}
 \nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0
 \end{equation}
@@ -181,8 +181,8 @@
 operators, $i.e.$
 \begin{align} 
-\label{DOM_di_adj}
+\label{eq:DOM_di_adj}
 \sum\limits_i { a_i \;\delta _i \left[ b \right]} 
 	&\equiv -\sum\limits_i {\delta _{i+1/2} \left[ a \right]\;b_{i+1/2} }      \\
-\label{DOM_mi_adj}
+\label{eq:DOM_mi_adj}
 \sum\limits_i { a_i \;\overline b^{\,i}} 
 	& \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} } 
@@ -192,5 +192,5 @@
 $\delta_i^*=\delta_{i+1/2}$ and 
 ${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively. 
-These two properties will be used extensively in the Appendix~\ref{Apdx_C} to 
+These two properties will be used extensively in the \autoref{apdx:C} to 
 demonstrate integral conservative properties of the discrete formulation chosen.
 
@@ -199,10 +199,10 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Numerical indexing}
-\label{DOM_Num_Index}
+\label{subsec:DOM_Num_Index}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!tb]  \begin{center}
 \includegraphics[width=0.90\textwidth]{Fig_index_hor}
-\caption{   \protect\label{Fig_index_hor}    
+\caption{   \protect\label{fig:index_hor}    
 Horizontal integer indexing used in the \textsc{Fortran} code. The dashed area indicates 
 the cell in which variables contained in arrays have the same $i$- and $j$-indices}
@@ -211,5 +211,5 @@
 
 The array representation used in the \textsc{Fortran} code requires an integer 
-indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is 
+indexing while the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is 
 associated with the use of integer values for $t$-points and both integer and 
 integer and a half values for all the other points. Therefore a specific integer 
@@ -222,9 +222,9 @@
 % -----------------------------------
 \subsubsection{Horizontal indexing}
-\label{DOM_Num_Index_hor}
-
-The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}. 
+\label{subsec:DOM_Num_Index_hor}
+
+The indexing in the horizontal plane has been chosen as shown in \autoref{fig:index_hor}. 
 For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point 
-(northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}). 
+(northward $v$-point) have the same index (see the dashed area in \autoref{fig:index_hor}). 
 A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices.
 
@@ -233,18 +233,18 @@
 % -----------------------------------
 \subsubsection{Vertical indexing}
-\label{DOM_Num_Index_vertical}
+\label{subsec:DOM_Num_Index_vertical}
 
 In the vertical, the chosen indexing requires special attention since the 
 $k$-axis is re-orientated downward in the \textsc{Fortran} code compared 
-to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}. 
+to the indexing used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}. 
 The sea surface corresponds to the $w$-level $k=1$ which is the same index 
-as $t$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$) 
+as $t$-level just below (\autoref{fig:index_vert}). The last $w$-level ($k=jpk$) 
 either corresponds to the ocean floor or is inside the bathymetry while the last 
-$t$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that 
+$t$-level is always inside the bathymetry (\autoref{fig:index_vert}). Note that 
 for an increasing $k$ index, a $w$-point and the $t$-point just below have the 
 same $k$ index, in opposition to what is done in the horizontal plane where 
 it is the $t$-point and the nearest velocity points in the direction of the horizontal 
 axis that have the same $i$ or $j$ index (compare the dashed area in 
-Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are 
+\autoref{fig:index_hor} and \autoref{fig:index_vert}). Since the scale factors are 
 chosen to be strictly positive, a \emph{minus sign} appears in the \textsc{Fortran} 
 code \emph{before all the vertical derivatives} of the discrete equations given in 
@@ -254,5 +254,5 @@
 \begin{figure}[!pt]    \begin{center}
 \includegraphics[width=.90\textwidth]{Fig_index_vert}
-\caption{ \protect\label{Fig_index_vert}     
+\caption{ \protect\label{fig:index_vert}     
 Vertical integer indexing used in the \textsc{Fortran } code. Note that 
 the $k$-axis is orientated downward. The dashed area indicates the cell in 
@@ -265,5 +265,5 @@
 % -----------------------------------
 \subsubsection{Domain size}
-\label{DOM_size}
+\label{subsec:DOM_size}
 
 The total size of the computational domain is set by the parameters \np{jpiglo}, 
@@ -273,5 +273,5 @@
 %%%
 Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when the code is 
-run in parallel using domain decomposition (\key{mpp\_mpi} defined, see \S\ref{LBC_mpp}).
+run in parallel using domain decomposition (\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}).
 
 
@@ -282,10 +282,10 @@
 % ================================================================
 \section{Needed fields}
-\label{DOM_fields}
+\label{sec:DOM_fields}
 The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined 
 by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
 The grid-points are located at integer or integer and a half values of as indicated 
-in Table~\ref{Tab_cell}. The associated scale factors are defined using the  
-analytical first derivative of the transformation \eqref{Eq_scale_factors}. 
+in \autoref{tab:cell}. The associated scale factors are defined using the  
+analytical first derivative of the transformation \autoref{eq:scale_factors}. 
 Necessary fields for configuration definition are: \\
 Geographic position :
@@ -316,5 +316,5 @@
 % -------------------------------------------------------------------------------------------------------------
 %\subsection{List of needed fields to build DOMAIN}
-%\label{DOM_fields_list}
+%\label{subsec:DOM_fields_list}
 
 
@@ -323,5 +323,5 @@
 % ================================================================
 \section{Horizontal grid mesh (\protect\mdl{domhgr})}
-\label{DOM_hgr}
+\label{sec:DOM_hgr}
 
 % -------------------------------------------------------------------------------------------------------------
@@ -329,11 +329,11 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Coordinates and scale factors}
-\label{DOM_hgr_coord_e}
+\label{subsec:DOM_hgr_coord_e}
 
 The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined 
 by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
 The grid-points are located at integer or integer and a half values of as indicated 
-in Table~\ref{Tab_cell}. The associated scale factors are defined using the 
-analytical first derivative of the transformation \eqref{Eq_scale_factors}. These 
+in \autoref{tab:cell}. The associated scale factors are defined using the 
+analytical first derivative of the transformation \autoref{eq:scale_factors}. These 
 definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, which 
 provide the horizontal and vertical meshes, respectively. This section deals with 
@@ -343,5 +343,5 @@
 analytical expressions of the longitude $\lambda$ and  latitude $\varphi$ as a 
 function of  $(i,j)$. The horizontal scale factors are calculated using 
-\eqref{Eq_scale_factors}. For example, when the longitude and latitude are 
+\autoref{eq:scale_factors}. For example, when the longitude and latitude are 
 function of a single value ($i$ and $j$, respectively) (geographical configuration 
 of the mesh), the horizontal mesh definition reduces to define the wanted 
@@ -382,9 +382,9 @@
 allowing the user to set arbitrary jumps in thickness between adjacent layers) 
 \citep{Treguier1996}. An example of the effect of such a choice is shown in 
-Fig.~\ref{Fig_zgr_e3}.
+\autoref{fig:zgr_e3}.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]     \begin{center}
 \includegraphics[width=0.90\textwidth]{Fig_zgr_e3}
-\caption{ \protect\label{Fig_zgr_e3}    
+\caption{ \protect\label{fig:zgr_e3}    
 Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
 and (b) analytically derived grid-point position and scale factors. 
@@ -401,5 +401,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Choice of horizontal grid}
-\label{DOM_hgr_msh_choice}
+\label{subsec:DOM_hgr_msh_choice}
 
 
@@ -408,5 +408,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Output grid files}
-\label{DOM_hgr_files}
+\label{subsec:DOM_hgr_files}
 
 All the arrays relating to a particular ocean model configuration (grid-point 
@@ -426,5 +426,5 @@
 % ================================================================
 \section{Vertical grid (\protect\mdl{domzgr})}
-\label{DOM_zgr}
+\label{sec:DOM_zgr}
 %-----------------------------------------nam_zgr & namdom-------------------------------------------
 %\forfile{../namelists/namzgr} 
@@ -444,5 +444,5 @@
 \begin{figure}[!tb]    \begin{center}
 \includegraphics[width=1.0\textwidth]{Fig_z_zps_s_sps}
-\caption{  \protect\label{Fig_z_zps_s_sps}   
+\caption{  \protect\label{fig:z_zps_s_sps}   
 The ocean bottom as seen by the model: 
 (a) $z$-coordinate with full step, 
@@ -451,5 +451,5 @@
 (d) hybrid $s-z$ coordinate, 
 (e) hybrid $s-z$ coordinate with partial step, and 
-(f) same as (e) but in the non-linear free surface (\np{ln\_linssh}\forcode{ = .false.}). 
+(f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}\forcode{ = .false.}). 
 Note that the non-linear free surface can be used with any of the 
 5 coordinates (a) to (e).}
@@ -460,11 +460,11 @@
 must be done once of all at the beginning of an experiment. It is not intended as an 
 option which can be enabled or disabled in the middle of an experiment. Three main 
-choices are offered (Fig.~\ref{Fig_z_zps_s_sps}a to c): $z$-coordinate with full step 
+choices are offered (\autoref{fig:z_zps_s_sps}a to c): $z$-coordinate with full step 
 bathymetry (\np{ln\_zco}\forcode{ = .true.}), $z$-coordinate with partial step bathymetry 
 (\np{ln\_zps}\forcode{ = .true.}), or generalized, $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}). 
 Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate 
-(Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). By default a non-linear free surface is used:
+(\autoref{fig:z_zps_s_sps}d and \autoref{fig:z_zps_s_sps}e). By default a non-linear free surface is used:
 the coordinate follow the time-variation of the free surface so that the transformation is time dependent: 
-$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). When a linear free surface is assumed (\np{ln\_linssh}\forcode{ = .true.}), 
+$z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}f). When a linear free surface is assumed (\np{ln\_linssh}\forcode{ = .true.}), 
 the vertical coordinate are fixed in time, but the seawater can move up and down across the z=0 surface 
 (in other words, the top of the ocean in not a rigid-lid). 
@@ -513,5 +513,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Meter bathymetry}
-\label{DOM_bathy}
+\label{subsec:DOM_bathy}
 
 Three options are possible for defining the bathymetry, according to the 
@@ -541,5 +541,5 @@
 This is unnecessary when the ocean is forced by fixed atmospheric conditions, 
 so these seas can be removed from the ocean domain. The user has the option 
-to set the bathymetry in closed seas to zero (see \S\ref{MISC_closea}), but the 
+to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}), but the 
 code has to be adapted to the user's configuration. 
 
@@ -549,13 +549,13 @@
 \subsection[$Z$-coordinate (\protect\np{ln\_zco}\forcode{ = .true.}) and ref. coordinate]
 				{$Z$-coordinate (\protect\np{ln\_zco}\forcode{ = .true.}) and reference coordinate}
-\label{DOM_zco}
+\label{subsec:DOM_zco}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!tb]    \begin{center}
 \includegraphics[width=0.90\textwidth]{Fig_zgr}
-\caption{ \protect\label{Fig_zgr}    
+\caption{ \protect\label{fig:zgr}    
 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 \eqref{DOM_zgr_ana} using \eqref{DOM_zgr_coef} in $z$-coordinate.}
+from \autoref{eq:DOM_zgr_ana} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -563,5 +563,5 @@
 The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ 
 and $gdepw_0$ for $t$- and $w$-points, respectively. As indicated on 
-Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the 
+\autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the 
 ocean surface. There are at most \jp{jpk}-1 $t$-points inside the ocean, the 
 additional $t$-point at $jk=jpk$ is below the sea floor and is not used. 
@@ -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{DOM_zgr_ana}
+\begin{equation} \label{eq:DOM_zgr_ana}
 \begin{split}
  z_0 (k) 	&= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\ 
@@ -588,9 +588,9 @@
 expression allows us to define a nearly uniform vertical location of levels at the 
 ocean top and bottom with a smooth hyperbolic tangent transition in between 
-(Fig.~\ref{Fig_zgr}).
+(\autoref{fig:zgr}).
 
 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{DOM_zgr_ana}
+\begin{equation} \label{eq:DOM_zgr_ana}
 \begin{split}
  e_3^T(k) &= z_W (k+1) - z_W (k)   \\
@@ -605,5 +605,5 @@
 surface (bottom) layers and a depth which varies from 0 at the sea surface to a 
 minimum of $-5000~m$. This leads to the following conditions:
-\begin{equation} \label{DOM_zgr_coef}
+\begin{equation} \label{eq:DOM_zgr_coef}
 \begin{split}
  e_3 (1+1/2)		&=10. \\ 
@@ -616,10 +616,10 @@
 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 
-\eqref{DOM_zgr_ana} have been determined such that \eqref{DOM_zgr_coef} is 
+\autoref{eq:DOM_zgr_ana} have been determined such that \autoref{eq:DOM_zgr_coef} is 
 satisfied, through an optimisation procedure using a bisection method. For the first 
 standard ORCA2 vertical grid this led to the following values: $h_{sur} =4762.96$, 
 $h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$. The resulting depths and 
-scale factors as a function of the model levels are shown in Fig.~\ref{Fig_zgr} and 
-given in Table \ref{Tab_orca_zgr}. Those values correspond to the parameters 
+scale factors as a function of the model levels are shown in \autoref{fig:zgr} and 
+given in \autoref{tab:orca_zgr}. Those values correspond to the parameters 
 \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. 
 
@@ -675,7 +675,7 @@
 31	& 	\textbf{5250.23}&	5000.00 &	\textbf{500.56} &	500.33 \\ \hline
 \end{tabular} \end{center} 
-\caption{ \protect\label{Tab_orca_zgr}   
+\caption{ \protect\label{tab:orca_zgr}   
 Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed 
-from \eqref{DOM_zgr_ana} using the coefficients given in \eqref{DOM_zgr_coef}}
+from \autoref{eq:DOM_zgr_ana} using the coefficients given in \autoref{eq:DOM_zgr_coef}}
 \end{table}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -685,5 +685,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}\forcode{ = .true.})}
-\label{DOM_zps}
+\label{subsec:DOM_zps}
 %--------------------------------------------namdom-------------------------------------------------------
 \forfile{../namelists/namdom} 
@@ -717,5 +717,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{$S$-coordinate (\protect\np{ln\_sco}\forcode{ = .true.})}
-\label{DOM_sco}
+\label{subsec:DOM_sco}
 %------------------------------------------nam_zgr_sco---------------------------------------------------
 %\forfile{../namelists/namzgr_sco} 
@@ -726,5 +726,5 @@
 function or its derivative, respectively:
 
-\begin{equation} \label{DOM_sco_ana}
+\begin{equation} \label{eq:DOM_sco_ana}
 \begin{split}
  z(k) 		&= h(i,j) \; z_0(k)	\\
@@ -737,5 +737,5 @@
 surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean 
 depth, since a mixed step-like and bottom-following representation of the 
-topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e) or an envelop bathymetry can be defined (Fig.~\ref{Fig_z_zps_s_sps}f).
+topography can be used (\autoref{fig:z_zps_s_sps}d-e) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}f).
 The namelist parameter \np{rn\_rmax} determines the slope at which the terrain-following coordinate intersects 
 the sea bed and becomes a pseudo z-coordinate. 
@@ -764,5 +764,5 @@
 \end{equation}
 
-\begin{equation} \label{DOM_sco_function}
+\begin{equation} \label{eq:DOM_sco_function}
 \begin{split}
 C(s)	&=  \frac{ \left[	  \tanh{ \left( \theta \, (s+b) \right)} 
@@ -784,5 +784,5 @@
 \begin{figure}[!ht]    \begin{center}
 \includegraphics[width=1.0\textwidth]{Fig_sco_function}
-\caption{  \protect\label{Fig_sco_function}   
+\caption{  \protect\label{fig:sco_function}   
 Examples of the stretching function applied to a seamount; from left to right: 
 surface, surface and bottom, and bottom intensified resolutions}
@@ -794,5 +794,5 @@
 are the surface and bottom control parameters such that $0\leqslant \theta \leqslant 20$, and 
 $0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom 
-increase of the vertical resolution (Fig.~\ref{Fig_sco_function}).
+increase of the vertical resolution (\autoref{fig:sco_function}).
 
 Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows 
@@ -807,10 +807,10 @@
 The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate:
 
-\begin{equation} \label{DOM_gamma_deriv}
+\begin{equation} \label{eq:DOM_gamma_deriv}
 \gamma= A\left(\sigma-\frac{1}{2}\left(\sigma^{2}+f\left(\sigma\right)\right)\right)+B\left(\sigma^{3}-f\left(\sigma\right)\right)+f\left(\sigma\right)
 \end{equation}
 
 Where:
-\begin{equation} \label{DOM_gamma}
+\begin{equation} \label{eq:DOM_gamma}
 f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2} \quad \text{ and } \quad \sigma = \frac{k}{n-1} 
 \end{equation}
@@ -821,5 +821,5 @@
 and bottom depths. The bottom cell depth in this example is given as a function of water depth:
 
-\begin{equation} \label{DOM_zb}
+\begin{equation} \label{eq:DOM_zb}
 Z_b= h a + b
 \end{equation}
@@ -831,5 +831,5 @@
    \includegraphics[width=1.0\textwidth]{FIG_DOM_compare_coordinates_surface}
         \caption{A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), a 50 level $Z$-coordinate (contoured surfaces) and the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface 100m for a idealised bathymetry that goes from 50m to 5500m depth. For clarity every third coordinate surface is shown.}
-    \label{fig_compare_coordinates_surface}
+    \label{fig:fig_compare_coordinates_surface}
 \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -845,5 +845,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{$Z^*$- or $S^*$-coordinate (\protect\np{ln\_linssh}\forcode{ = .false.}) }
-\label{DOM_zgr_star}
+\label{subsec:DOM_zgr_star}
 
 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 
@@ -855,5 +855,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Level bathymetry and mask}
-\label{DOM_msk}
+\label{subsec:DOM_msk}
 
 Whatever the vertical coordinate used, the model offers the possibility of 
@@ -892,5 +892,5 @@
 
 Note that, without ice shelves cavities, masks at $t-$ and $w-$points are identical with 
-the numerical indexing used (\S~\ref{DOM_Num_Index}). Nevertheless, $wmask$ are required 
+the numerical indexing used (\autoref{subsec:DOM_Num_Index}). Nevertheless, $wmask$ are required 
 with oceean cavities to deal with the top boundary (ice shelf/ocean interface) 
 exactly in the same way as for the bottom boundary. 
@@ -900,5 +900,5 @@
 case of an east-west cyclical boundary condition, \textit{mbathy} has its last 
 column equal to the second one and its first column equal to the last but one 
-(and so too the mask arrays) (see \S~\ref{LBC_jperio}).
+(and so too the mask arrays) (see \autoref{fig:LBC_jperio}).
 
 
@@ -907,5 +907,5 @@
 % ================================================================
 \section{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})}
-\label{DTA_tsd}
+\label{sec:DTA_tsd}
 %-----------------------------------------namtsd-------------------------------------------
 \forfile{../namelists/namtsd} 
@@ -918,5 +918,5 @@
 \item[\np{ln\_tsd\_init}\forcode{ = .true.}] use a T and S input files that can be given on the model grid itself or 
 on their native input data grid. In the latter case, the data will be interpolated on-the-fly both in the 
-horizontal and the vertical to the model grid (see \S~\ref{SBC_iof}). The information relative to the 
+horizontal and the vertical to the model grid (see \autoref{subsec:SBC_iof}). The information relative to the 
 input files are given in the \np{sn\_tem} and \np{sn\_sal} structures. 
 The computation is done in the \mdl{dtatsd} module.
