Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_TRA.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_TRA.tex	(revision 9394)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_TRA.tex	(revision 9407)
@@ -5,5 +5,5 @@
 % ================================================================
 \chapter{Ocean Tracers (TRA)}
-\label{TRA}
+\label{chap:TRA}
 \minitoc
 
@@ -17,5 +17,5 @@
 %$\ $\newline    % force a new ligne
 
-Using the representation described in Chap.~\ref{DOM}, several semi-discrete 
+Using the representation described in \autoref{chap:DOM}, several semi-discrete 
 space forms of the tracer equations are available depending on the vertical 
 coordinate used and on the physics used. In all the equations presented 
@@ -40,5 +40,5 @@
 require complex inputs and complex calculations ($e.g.$ bulk formulae, estimation 
 of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and 
-described in chapters \S\ref{SBC}, \S\ref{LDF} and  \S\ref{ZDF}, respectively. 
+described in \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. 
 Note that \mdl{tranpc}, the non-penetrative convection module, although 
 located in the NEMO/OPA/TRA directory as it directly modifies the tracer fields, 
@@ -57,5 +57,5 @@
 
 The user has the option of extracting each tendency term on the RHS of the tracer 
-equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}\forcode{ = .true.}), as described in Chap.~\ref{DIA}.
+equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}\forcode{ = .true.}), as described in \autoref{chap:DIA}.
 
 $\ $\newline    % force a new ligne
@@ -64,5 +64,5 @@
 % ================================================================
 \section{Tracer advection (\protect\mdl{traadv})}
-\label{TRA_adv}
+\label{sec:TRA_adv}
 %------------------------------------------namtra_adv-----------------------------------------------------
 \forfile{../namelists/namtra_adv}
@@ -72,5 +72,5 @@
 the advection tendency of a tracer is expressed in flux form, 
 $i.e.$ as the divergence of the advective fluxes. Its discrete expression is given by :
-\begin{equation} \label{Eq_tra_adv}
+\begin{equation} \label{eq:tra_adv}
 ADV_\tau =-\frac{1}{b_t} \left( 
 \;\delta _i \left[ e_{2u}\,e_{3u} \;  u\; \tau _u  \right]
@@ -79,5 +79,5 @@
 \end{equation}
 where $\tau$ is either T or S, and $b_t= e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells. 
-The flux form in \eqref{Eq_tra_adv} 
+The flux form in \autoref{eq:tra_adv} 
 implicitly requires the use of the continuity equation. Indeed, it is obtained
 by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 
@@ -87,10 +87,10 @@
 advection tendency so that it is consistent with the continuity equation in order to 
 enforce the conservation properties of the continuous equations. In other words, 
-by setting $\tau = 1$ in (\ref{Eq_tra_adv}) we recover the discrete form of 
+by setting $\tau = 1$ in (\autoref{eq:tra_adv}) we recover the discrete form of 
 the continuity equation which is used to calculate the vertical velocity.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t] 	 \begin{center}
 \includegraphics[width=0.9\textwidth]{Fig_adv_scheme}
-\caption{	\protect\label{Fig_adv_scheme} 
+\caption{	\protect\label{fig:adv_scheme} 
 Schematic representation of some ways used to evaluate the tracer value 
 at $u$-point and the amount of tracer exchanged between two neighbouring grid 
@@ -107,5 +107,5 @@
 The key difference between the advection schemes available in \NEMO is the choice 
 made in space and time interpolation to define the value of the tracer at the 
-velocity points (Fig.~\ref{Fig_adv_scheme}). 
+velocity points (\autoref{fig:adv_scheme}). 
 
 Along solid lateral and bottom boundaries a zero tracer flux is automatically 
@@ -131,9 +131,9 @@
 height, two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}.
 
-The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco}) 
+The velocity field that appears in (\autoref{eq:tra_adv}) and (\autoref{eq:tra_adv_zco}) 
 is the centred (\textit{now}) \textit{effective} ocean velocity, $i.e.$ the \textit{eulerian} velocity
-(see Chap.~\ref{DYN}) plus the eddy induced velocity (\textit{eiv}) 
+(see \autoref{chap:DYN}) plus the eddy induced velocity (\textit{eiv}) 
 and/or the mixed layer eddy induced velocity (\textit{eiv})
-when those parameterisations are used (see Chap.~\ref{LDF}).
+when those parameterisations are used (see \autoref{chap:LDF}).
 
 Several tracer advection scheme are proposed, namely 
@@ -174,5 +174,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{CEN: Centred scheme (\protect\np{ln\_traadv\_cen}\forcode{ = .true.})}
-\label{TRA_adv_cen}
+\label{subsec:TRA_adv_cen}
 
 %        2nd order centred scheme  
@@ -186,5 +186,5 @@
 is evaluated as the mean of the two neighbouring $T$-point values. 
 For example, in the $i$-direction :
-\begin{equation} \label{Eq_tra_adv_cen2}
+\begin{equation} \label{eq:tra_adv_cen2}
 \tau _u^{cen2} =\overline T ^{i+1/2}
 \end{equation}
@@ -195,8 +195,8 @@
 produce a sensible solution. The associated time-stepping is performed using 
 a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in 
-(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value. 
+(\autoref{eq:tra_adv_cen2}) is the \textit{now} tracer value. 
 
 Note that using the CEN2, the overall tracer advection is of second 
-order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2}) 
+order accuracy since both (\autoref{eq:tra_adv}) and (\autoref{eq:tra_adv_cen2}) 
 have this order of accuracy.
 
@@ -206,5 +206,5 @@
 a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points. 
 For example, in the $i$-direction:
-\begin{equation} \label{Eq_tra_adv_cen4}
+\begin{equation} \label{eq:tra_adv_cen4}
 \tau _u^{cen4} 
 =\overline{   T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right]   }^{\,i+1/2}
@@ -219,5 +219,5 @@
 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme 
 but a $4^{th}$ order evaluation of advective fluxes, since the divergence of 
-advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order. 
+advective fluxes \autoref{eq:tra_adv} is kept at $2^{nd}$ order. 
 The expression \textit{$4^{th}$ order scheme} used in oceanographic literature 
 is usually associated with the scheme presented here. 
@@ -232,5 +232,5 @@
 to produce a sensible solution. 
 As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction 
-with an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
+with an Asselin time-filter, so $T$ in (\autoref{eq:tra_adv_cen4}) is the \textit{now} tracer.
 
 At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), 
@@ -245,5 +245,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}\forcode{ = .true.})}
-\label{TRA_adv_tvd}
+\label{subsec:TRA_adv_tvd}
 
 The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}\forcode{ = .true.}. 
@@ -254,5 +254,5 @@
 In FCT formulation, the tracer at velocity points is evaluated using a combination of 
 an upstream and a centred scheme. For example, in the $i$-direction :
-\begin{equation} \label{Eq_tra_adv_fct}
+\begin{equation} \label{eq:tra_adv_fct}
 \begin{split}
 \tau _u^{ups}&= \begin{cases}
@@ -280,7 +280,7 @@
 by vertical advection \citep{Lemarie_OM2015}. Note that in this case, a similar split-explicit 
 time stepping should be used on vertical advection of momentum to insure a better stability
-(see \S\ref{DYN_zad}).
-
-For stability reasons (see \S\ref{STP}), $\tau _u^{cen}$ is evaluated in (\ref{Eq_tra_adv_fct}) 
+(see \autoref{subsec:DYN_zad}).
+
+For stability reasons (see \autoref{chap:STP}), $\tau _u^{cen}$ is evaluated in (\autoref{eq:tra_adv_fct}) 
 using the \textit{now} tracer while $\tau _u^{ups}$ is evaluated using the \textit{before} tracer. In other words, 
 the advective part of the scheme is time stepped with a leap-frog scheme 
@@ -291,5 +291,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}\forcode{ = .true.})}
-\label{TRA_adv_mus}
+\label{subsec:TRA_adv_mus}
 
 The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}\forcode{ = .true.}. 
@@ -298,6 +298,6 @@
 MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}. In its formulation, the tracer at velocity points 
 is evaluated assuming a linear tracer variation between two $T$-points 
-(Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction :
-\begin{equation} \label{Eq_tra_adv_mus}
+(\autoref{fig:adv_scheme}). For example, in the $i$-direction :
+\begin{equation} \label{eq:tra_adv_mus}
    \tau _u^{mus} = \left\{      \begin{aligned}
          &\tau _i  &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right)
@@ -323,5 +323,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}\forcode{ = .true.})}
-\label{TRA_adv_ubs}
+\label{subsec:TRA_adv_ubs}
 
 The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}\forcode{ = .true.}. 
@@ -332,5 +332,5 @@
 third order scheme based on an upstream-biased parabolic interpolation.  
 For example, in the $i$-direction :
-\begin{equation} \label{Eq_tra_adv_ubs}
+\begin{equation} \label{eq:tra_adv_ubs}
    \tau _u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{      
    \begin{aligned}
@@ -355,6 +355,6 @@
 or a $4^th$ order COMPACT scheme (\np{nn\_cen\_v}\forcode{ = 2 or 4}).
 
-For stability reasons  (see \S\ref{STP}),
-the first term  in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order 
+For stability reasons  (see \autoref{chap:STP}),
+the first term  in \autoref{eq:tra_adv_ubs} (which corresponds to a second order 
 centred scheme) is evaluated using the \textit{now} tracer (centred in time) 
 while the second term (which is the diffusive part of the scheme), is 
@@ -362,5 +362,5 @@
 This choice is discussed by \citet{Webb_al_JAOT98} in the context of the 
 QUICK advection scheme. UBS and QUICK schemes only differ 
-by one coefficient. Replacing 1/6 with 1/8 in \eqref{Eq_tra_adv_ubs} 
+by one coefficient. Replacing 1/6 with 1/8 in \autoref{eq:tra_adv_ubs} 
 leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 
 This option is not available through a namelist parameter, since the 
@@ -368,6 +368,6 @@
 substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme.
 
-Note that it is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
-\begin{equation} \label{Eq_traadv_ubs2}
+Note that it is straightforward to rewrite \autoref{eq:tra_adv_ubs} as follows:
+\begin{equation} \label{eq:traadv_ubs2}
 \tau _u^{ubs} = \tau _u^{cen4} + \frac{1}{12} \left\{	 
    \begin{aligned}
@@ -377,5 +377,5 @@
 \end{equation}
 or equivalently 
-\begin{equation} \label{Eq_traadv_ubs2b}
+\begin{equation} \label{eq:traadv_ubs2b}
 u_{i+1/2} \ \tau _u^{ubs} 
 =u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\delta _i\left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2}
@@ -383,13 +383,13 @@
 \end{equation}
 
-\eqref{Eq_traadv_ubs2} has several advantages. Firstly, it clearly reveals 
+\autoref{eq:traadv_ubs2} has several advantages. Firstly, it clearly reveals 
 that the UBS scheme is based on the fourth order scheme to which an 
 upstream-biased diffusion term is added. Secondly, this emphasises that the 
 $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has 
-to be evaluated at the \emph{now} time step using \eqref{Eq_tra_adv_ubs}. 
+to be evaluated at the \emph{now} time step using \autoref{eq:tra_adv_ubs}. 
 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy 
 coefficient which is simply proportional to the velocity:
  $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note the current version of NEMO uses 
-the computationally more efficient formulation \eqref{Eq_tra_adv_ubs}.
+the computationally more efficient formulation \autoref{eq:tra_adv_ubs}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -397,5 +397,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}\forcode{ = .true.})}
-\label{TRA_adv_qck}
+\label{subsec:TRA_adv_qck}
 
 The Quadratic Upstream Interpolation for Convective Kinematics with 
@@ -423,5 +423,5 @@
 % ================================================================
 \section{Tracer lateral diffusion (\protect\mdl{traldf})}
-\label{TRA_ldf}
+\label{sec:TRA_ldf}
 %-----------------------------------------nam_traldf------------------------------------------------------
 \forfile{../namelists/namtra_ldf}
@@ -434,12 +434,12 @@
 $(iii)$ some specific options related to the rotated operators ($i.e.$ non-iso-level operator), and 
 $(iv)$  the specification of eddy diffusivity coefficient (either constant or variable in space and time).
-Item $(iv)$ will be described in Chap.\ref{LDF} .
+Item $(iv)$ will be described in \autoref{chap:LDF} .
 The direction along which the operators act is defined through the slope between this direction and the iso-level surfaces.
-The slope is computed in the \mdl{ldfslp} module and will also be described in Chap.~\ref{LDF}. 
+The slope is computed in the \mdl{ldfslp} module and will also be described in \autoref{chap:LDF}. 
 
 The lateral diffusion of tracers is evaluated using a forward scheme, 
 $i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time, 
 except for the pure vertical component that appears when a rotation tensor is used. 
-This latter component is solved implicitly together with the vertical diffusion term (see \S\ref{STP}). 
+This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}). 
 When \np{ln\_traldf\_msc}\forcode{ = .true.}, a Method of Stabilizing Correction is used in which 
 the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}.
@@ -450,5 +450,5 @@
 \subsection[Type of operator (\protect\np{ln\_traldf}\{\_NONE,\_lap,\_blp\}\})]
 			     {Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) } 
-\label{TRA_ldf_op}
+\label{subsec:TRA_ldf_op}
 
 Three operator options are proposed and, one and only one of them must be selected:
@@ -459,11 +459,11 @@
 \item [\np{ln\_traldf\_lap}\forcode{ = .true.}]: a laplacian operator is selected. This harmonic operator 
 takes the following expression:  $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $, 
-where the gradient operates along the selected direction (see \S\ref{TRA_ldf_dir}),
-and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see Chap.~\ref{LDF}).
+where the gradient operates along the selected direction (see \autoref{subsec:TRA_ldf_dir}),
+and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}).
 \item [\np{ln\_traldf\_blp}\forcode{ = .true.}]: a bilaplacian operator is selected. This biharmonic operator 
 takes the following expression:  
 $\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$ 
 where the gradient operats along the selected direction,
-and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$  (see Chap.~\ref{LDF}).
+and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$  (see \autoref{chap:LDF}).
 In the code, the bilaplacian operator is obtained by calling the laplacian twice.
 \end{description}
@@ -483,5 +483,5 @@
 \subsection[Action direction (\protect\np{ln\_traldf}\{\_lev,\_hor,\_iso,\_triad\})]
 			     {Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) } 
-\label{TRA_ldf_dir}
+\label{subsec:TRA_ldf_dir}
 
 The choice of a direction of action determines the form of operator used. 
@@ -509,8 +509,8 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Iso-level (bi-)laplacian operator ( \protect\np{ln\_traldf\_iso}) }
-\label{TRA_ldf_lev}
+\label{subsec:TRA_ldf_lev}
 
 The laplacian diffusion operator acting along the model (\textit{i,j})-surfaces is given by: 
-\begin{equation} \label{Eq_tra_ldf_lap}
+\begin{equation} \label{eq:tra_ldf_lap}
 D_t^{lT} =\frac{1}{b_t} \left( \;
    \delta _{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta _{i+1/2} [T] \right] 
@@ -533,7 +533,7 @@
 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}), tracers in horizontally 
 adjacent cells are located at different depths in the vicinity of the bottom. 
-In this case, horizontal derivatives in (\ref{Eq_tra_ldf_lap}) at the bottom level 
+In this case, horizontal derivatives in (\autoref{eq:tra_ldf_lap}) at the bottom level 
 require a specific treatment. They are calculated in the \mdl{zpshde} module, 
-described in \S\ref{TRA_zpshde}.
+described in \autoref{sec:TRA_zpshde}.
 
 
@@ -542,13 +542,13 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Standard and triad (bi-)laplacian operator}
-\label{TRA_ldf_iso_triad}
+\label{subsec:TRA_ldf_iso_triad}
 
 %&&    Standard rotated (bi-)laplacian operator
 %&& ----------------------------------------------
 \subsubsection{Standard rotated (bi-)laplacian operator (\protect\mdl{traldf\_iso})}
-\label{TRA_ldf_iso}
+\label{subsec:TRA_ldf_iso}
 The general form of the second order lateral tracer subgrid scale physics 
-(\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and $s$-coordinates:
-\begin{equation} \label{Eq_tra_ldf_iso}
+(\autoref{eq:PE_zdf}) takes the following semi-discrete space form in $z$- and $s$-coordinates:
+\begin{equation} \label{eq:tra_ldf_iso}
 \begin{split}
  D_T^{lT} = \frac{1}{b_t}   & \left\{   \,\;\delta_i \left[   A_u^{lT}  \left( 
@@ -576,12 +576,12 @@
 in addition to \np{ln\_traldf\_lap}\forcode{ = .true.}, we have \np{ln\_traldf\_iso}\forcode{ = .true.}, 
 or both \np{ln\_traldf\_hor}\forcode{ = .true.} and \np{ln\_zco}\forcode{ = .true.}. The way these 
-slopes are evaluated is given in \S\ref{LDF_slp}. At the surface, bottom 
+slopes are evaluated is given in \autoref{sec:LDF_slp}. At the surface, bottom 
 and lateral boundaries, the turbulent fluxes of heat and salt are set to zero 
-using the mask technique (see \S\ref{LBC_coast}). 
-
-The operator in \eqref{Eq_tra_ldf_iso} involves both lateral and vertical 
+using the mask technique (see \autoref{sec:LBC_coast}). 
+
+The operator in \autoref{eq:tra_ldf_iso} involves both lateral and vertical 
 derivatives. For numerical stability, the vertical second derivative must 
 be solved using the same implicit time scheme as that used in the vertical 
-physics (see \S\ref{TRA_zdf}). For computer efficiency reasons, this term 
+physics (see \autoref{sec:TRA_zdf}). For computer efficiency reasons, this term 
 is not computed in the \mdl{traldf\_iso} module, but in the \mdl{trazdf} module 
 where, if iso-neutral mixing is used, the vertical mixing coefficient is simply 
@@ -590,29 +590,29 @@
 This formulation conserves the tracer but does not ensure the decrease 
 of the tracer variance. Nevertheless the treatment performed on the slopes 
-(see \S\ref{LDF}) allows the model to run safely without any additional 
+(see \autoref{chap:LDF}) allows the model to run safely without any additional 
 background horizontal diffusion \citep{Guilyardi_al_CD01}. 
 
 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}), the horizontal derivatives 
-at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific treatment. 
-They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
+at the bottom level in \autoref{eq:tra_ldf_iso} require a specific treatment. 
+They are calculated in module zpshde, described in \autoref{sec:TRA_zpshde}.
 
 %&&     Triad rotated (bi-)laplacian operator
 %&&  -------------------------------------------
 \subsubsection{Triad rotated (bi-)laplacian operator (\protect\np{ln\_traldf\_triad})}
-\label{TRA_ldf_triad}
-
-If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}\forcode{ = .true.} ; see App.\ref{sec:triad}) 
+\label{subsec:TRA_ldf_triad}
+
+If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}\forcode{ = .true.} ; see \autoref{apdx:triad}) 
 
 An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases 
 is also available in \NEMO (\np{ln\_traldf\_grif}\forcode{ = .true.}). A complete description of 
-the algorithm is given in App.\ref{sec:triad}.
+the algorithm is given in \autoref{apdx:triad}.
 
 The lateral fourth order bilaplacian operator on tracers is obtained by 
-applying (\ref{Eq_tra_ldf_lap}) twice. The operator requires an additional assumption 
+applying (\autoref{eq:tra_ldf_lap}) twice. The operator requires an additional assumption 
 on boundary conditions: both first and third derivative terms normal to the 
 coast are set to zero.
 
 The lateral fourth order operator formulation on tracers is obtained by 
-applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption 
+applying (\autoref{eq:tra_ldf_iso}) twice. It requires an additional assumption 
 on boundary conditions: first and third derivative terms normal to the 
 coast, normal to the bottom and normal to the surface are set to zero. 
@@ -621,5 +621,5 @@
 %&& ----------------------------------------------
 \subsubsection{Option for the rotated operators}
-\label{TRA_ldf_options}
+\label{subsec:TRA_ldf_options}
 
 \np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators)
@@ -637,5 +637,5 @@
 % ================================================================
 \section{Tracer vertical diffusion (\protect\mdl{trazdf})}
-\label{TRA_zdf}
+\label{sec:TRA_zdf}
 %--------------------------------------------namzdf---------------------------------------------------------
 \forfile{../namelists/namzdf}
@@ -645,7 +645,7 @@
 The formulation of the vertical subgrid scale tracer physics is the same 
 for all the vertical coordinates, and is based on a laplacian operator. 
-The vertical diffusion operator given by (\ref{Eq_PE_zdf}) takes the 
+The vertical diffusion operator given by (\autoref{eq:PE_zdf}) takes the 
 following semi-discrete space form:
-\begin{equation} \label{Eq_tra_zdf}
+\begin{equation} \label{eq:tra_zdf}
 \begin{split}
 D^{vT}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vT}_w}{e_{3w}}  \delta_{k+1/2}[T] \;\right] 
@@ -658,15 +658,15 @@
 $A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is 
 parameterised ($i.e.$ \key{zdfddm} is defined). The way these coefficients 
-are evaluated is given in \S\ref{ZDF} (ZDF). Furthermore, when 
+are evaluated is given in \autoref{chap:ZDF} (ZDF). Furthermore, when 
 iso-neutral mixing is used, both mixing coefficients are increased 
 by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ 
-to account for the vertical second derivative of \eqref{Eq_tra_ldf_iso}. 
+to account for the vertical second derivative of \autoref{eq:tra_ldf_iso}. 
 
 At the surface and bottom boundaries, the turbulent fluxes of 
 heat and salt must be specified. At the surface they are prescribed 
-from the surface forcing and added in a dedicated routine (see \S\ref{TRA_sbc}), 
+from the surface forcing and added in a dedicated routine (see \autoref{subsec:TRA_sbc}), 
 whilst at the bottom they are set to zero for heat and salt unless 
 a geothermal flux forcing is prescribed as a bottom boundary 
-condition (see \S\ref{TRA_bbc}). 
+condition (see \autoref{subsec:TRA_bbc}). 
 
 The large eddy coefficient found in the mixed layer together with high 
@@ -684,5 +684,5 @@
 % ================================================================
 \section{External forcing}
-\label{TRA_sbc_qsr_bbc}
+\label{sec:TRA_sbc_qsr_bbc}
 
 % -------------------------------------------------------------------------------------------------------------
@@ -690,5 +690,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Surface boundary condition (\protect\mdl{trasbc})}
-\label{TRA_sbc}
+\label{subsec:TRA_sbc}
 
 The surface boundary condition for tracers is implemented in a separate 
@@ -703,13 +703,13 @@
 of the ocean is due both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) 
 and to the heat and salt content of the mass exchange. They are both included directly in $Q_{ns}$, 
-the surface heat flux, and $F_{salt}$, the surface salt flux (see \S\ref{SBC} for further details).
+the surface heat flux, and $F_{salt}$, the surface salt flux (see \autoref{chap:SBC} for further details).
 By doing this, the forcing formulation is the same for any tracer (including temperature and salinity).
 
-The surface module (\mdl{sbcmod}, see \S\ref{SBC}) provides the following 
+The surface module (\mdl{sbcmod}, see \autoref{chap:SBC}) provides the following 
 forcing fields (used on tracers):
 
 $\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 
 (i.e. the difference between the total surface heat flux and the fraction of the short wave flux that 
-penetrates into the water column, see \S\ref{TRA_qsr}) plus the heat content associated with 
+penetrates into the water column, see \autoref{subsec:TRA_qsr}) plus the heat content associated with 
 of the mass exchange with the atmosphere and lands.
 
@@ -720,11 +720,11 @@
 
 $\bullet$ \textit{rnf}, the mass flux associated with runoff 
-(see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
+(see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
 
 $\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, 
-(see \S\ref{SBC_isf} for further details on how the ice shelf melt is computed and applied).
+(see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied).
 
 The surface boundary condition on temperature and salinity is applied as follows:
-\begin{equation} \label{Eq_tra_sbc}
+\begin{equation} \label{eq:tra_sbc}
 \begin{aligned}
  &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns}       }^t  & \\ 
@@ -734,5 +734,5 @@
 where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps 
 ($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the 
-divergence of odd and even time step (see \S\ref{STP}).
+divergence of odd and even time step (see \autoref{chap:STP}).
 
 In the linear free surface case (\np{ln\_linssh}\forcode{ = .true.}), 
@@ -742,5 +742,5 @@
 would have resulted from a change in the volume of the first level.
 The resulting surface boundary condition is applied as follows:
-\begin{equation} \label{Eq_tra_sbc_lin}
+\begin{equation} \label{eq:tra_sbc_lin}
 \begin{aligned}
  &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }   
@@ -754,5 +754,5 @@
 In the linear free surface case, there is a small imbalance. The imbalance is larger 
 than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}. 
-This is the reason why the modified filter is not applied in the linear free surface case (see \S\ref{STP}).
+This is the reason why the modified filter is not applied in the linear free surface case (see \autoref{chap:STP}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -760,18 +760,18 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Solar radiation penetration (\protect\mdl{traqsr})}
-\label{TRA_qsr}
+\label{subsec:TRA_qsr}
 %--------------------------------------------namqsr--------------------------------------------------------
 \forfile{../namelists/namtra_qsr}
 %--------------------------------------------------------------------------------------------------------------
 
-Options are defined through the  \ngn{namtra\_qsr} namelist variables.
+Options are defined through the \ngn{namtra\_qsr} namelist variables.
 When the penetrative solar radiation option is used (\np{ln\_flxqsr}\forcode{ = .true.}), 
 the solar radiation penetrates the top few tens of meters of the ocean. If it is not used 
 (\np{ln\_flxqsr}\forcode{ = .false.}) all the heat flux is absorbed in the first ocean level. 
 Thus, in the former case a term is added to the time evolution equation of 
-temperature \eqref{Eq_PE_tra_T} and the surface boundary condition is 
+temperature \autoref{eq:PE_tra_T} and the surface boundary condition is 
 modified to take into account only the non-penetrative part of the surface 
 heat flux:
-\begin{equation} \label{Eq_PE_qsr}
+\begin{equation} \label{eq:PE_qsr}
 \begin{split}
 \frac{\partial T}{\partial t} &= {\ldots} + \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k}	\\
@@ -781,24 +781,24 @@
 where $Q_{sr}$ is the penetrative part of the surface heat flux ($i.e.$ the shortwave radiation) 
 and $I$ is the downward irradiance ($\left. I \right|_{z=\eta}=Q_{sr}$). 
-The additional term in \eqref{Eq_PE_qsr} is discretized as follows:
-\begin{equation} \label{Eq_tra_qsr}
+The additional term in \autoref{eq:PE_qsr} is discretized as follows:
+\begin{equation} \label{eq:tra_qsr}
 \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k} \equiv \frac{1}{\rho_o\, C_p\, e_{3t}} \delta_k \left[ I_w \right]
 \end{equation}
 
-The shortwave radiation,  $Q_{sr}$, consists of energy distributed across a wide spectral range. 
+The shortwave radiation, $Q_{sr}$, consists of energy distributed across a wide spectral range. 
 The ocean is strongly absorbing for wavelengths longer than 700~nm and these 
 wavelengths contribute to heating the upper few tens of centimetres. The fraction of $Q_{sr}$ 
 that resides in these almost non-penetrative wavebands, $R$, is $\sim 58\%$ (specified 
-through namelist parameter \np{rn\_abs}).  It is assumed to penetrate the ocean 
+through namelist parameter \np{rn\_abs}). It is assumed to penetrate the ocean 
 with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 
-of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the namtra\_qsr namelist).
+of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the \ngn{namtra\_qsr} namelist).
 For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy 
 propagates to larger depths where it contributes to 
 local heating. 
 The way this second part of the solar energy penetrates into the ocean depends on 
-which formulation is chosen. In the simple 2-waveband light penetration scheme  (\np{ln\_qsr\_2bd}\forcode{ = .true.}) 
+which formulation is chosen. In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}\forcode{ = .true.}) 
 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths, 
-leading to the following expression  \citep{Paulson1977}:
-\begin{equation} \label{Eq_traqsr_iradiance}
+leading to the following expression \citep{Paulson1977}:
+\begin{equation} \label{eq:traqsr_iradiance}
 I(z) = Q_{sr} \left[Re^{-z / \xi_0} + \left( 1-R\right) e^{-z / \xi_1} \right]
 \end{equation}
@@ -810,5 +810,5 @@
 Such assumptions have been shown to provide a very crude and simplistic 
 representation of observed light penetration profiles (\cite{Morel_JGR88}, see also 
-Fig.\ref{Fig_traqsr_irradiance}). Light absorption in the ocean depends on 
+\autoref{fig:traqsr_irradiance}). Light absorption in the ocean depends on 
 particle concentration and is spectrally selective. \cite{Morel_JGR88} has shown 
 that an accurate representation of light penetration can be provided by a 61 waveband 
@@ -819,5 +819,5 @@
 attenuation coefficient is fitted to the coefficients computed from the full spectral model 
 of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}), assuming 
-the same power-law relationship. As shown in Fig.\ref{Fig_traqsr_irradiance}, 
+the same power-law relationship. As shown in \autoref{fig:traqsr_irradiance}, 
 this formulation, called RGB (Red-Green-Blue), reproduces quite closely 
 the light penetration profiles predicted by the full spectal model, but with much greater 
@@ -842,5 +842,5 @@
 light limitation in PISCES or LOBSTER and the oceanic heating rate. 
 \end{description} 
-The trend in \eqref{Eq_tra_qsr} associated with the penetration of the solar radiation 
+The trend in \autoref{eq:tra_qsr} associated with the penetration of the solar radiation 
 is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. 
 
@@ -857,5 +857,5 @@
 \begin{figure}[!t] 	  \begin{center}
 \includegraphics[width=1.0\textwidth]{Fig_TRA_Irradiance}
-\caption{ 	 \protect\label{Fig_traqsr_irradiance}
+\caption{ 	 \protect\label{fig:traqsr_irradiance}
 Penetration profile of the downward solar irradiance calculated by four models. 
 Two waveband chlorophyll-independent formulation (blue), a chlorophyll-dependent 
@@ -870,5 +870,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Bottom boundary condition (\protect\mdl{trabbc})}
-\label{TRA_bbc}
+\label{subsec:TRA_bbc}
 %--------------------------------------------nambbc--------------------------------------------------------
 \forfile{../namelists/nambbc}
@@ -877,5 +877,5 @@
 \begin{figure}[!t] 	  \begin{center}
 \includegraphics[width=1.0\textwidth]{Fig_TRA_geoth}
-\caption{ 	\protect\label{Fig_geothermal}
+\caption{ 	\protect\label{fig:geothermal}
 Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}.
 It is inferred from the age of the sea floor and the formulae of \citet{Stein_Stein_Nat92}.}
@@ -902,5 +902,5 @@
 When  \np{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is 
 introduced which is provided in the \ifile{geothermal\_heating} NetCDF file 
-(Fig.\ref{Fig_geothermal}) \citep{Emile-Geay_Madec_OS09}.
+(\autoref{fig:geothermal}) \citep{Emile-Geay_Madec_OS09}.
 
 % ================================================================
@@ -908,5 +908,5 @@
 % ================================================================
 \section{Bottom boundary layer (\protect\mdl{trabbl} - \protect\key{trabbl})}
-\label{TRA_bbl}
+\label{sec:TRA_bbl}
 %--------------------------------------------nambbl---------------------------------------------------------
 \forfile{../namelists/nambbl}
@@ -943,9 +943,9 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}\forcode{ = 1})}
-\label{TRA_bbl_diff}
+\label{subsec:TRA_bbl_diff}
 
 When applying sigma-diffusion (\key{trabbl} defined and \np{nn\_bbl\_ldf} set to 1), 
 the diffusive flux between two adjacent cells at the ocean floor is given by 
-\begin{equation} \label{Eq_tra_bbl_diff}
+\begin{equation} \label{eq:tra_bbl_diff}
 {\rm {\bf F}}_\sigma=A_l^\sigma \; \nabla_\sigma T
 \end{equation} 
@@ -953,5 +953,5 @@
 and  $A_l^\sigma$ the lateral diffusivity in the BBL. Following \citet{Beckmann_Doscher1997}, 
 the latter is prescribed with a spatial dependence, $i.e.$ in the conditional form
-\begin{equation} \label{Eq_tra_bbl_coef}
+\begin{equation} \label{eq:tra_bbl_coef}
 A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l}
  A_{bbl}  \quad \quad   \mbox{if}  \quad   \nabla_\sigma \rho  \cdot  \nabla H<0 \\ 
@@ -962,10 +962,10 @@
 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist 
 parameter \np{rn\_ahtbbl} and usually set to a value much larger 
-than the one used for lateral mixing in the open ocean. The constraint in \eqref{Eq_tra_bbl_coef} 
+than the one used for lateral mixing in the open ocean. The constraint in \autoref{eq:tra_bbl_coef} 
 implies that sigma-like diffusion only occurs when the density above the sea floor, at the top of 
-the slope, is larger than in the deeper ocean (see green arrow in Fig.\ref{Fig_bbl}). 
+the slope, is larger than in the deeper ocean (see green arrow in \autoref{fig:bbl}). 
 In practice, this constraint is applied separately in the two horizontal directions, 
-and the density gradient in \eqref{Eq_tra_bbl_coef} is evaluated with the log gradient formulation: 
-\begin{equation} \label{Eq_tra_bbl_Drho}
+and the density gradient in \autoref{eq:tra_bbl_coef} is evaluated with the log gradient formulation: 
+\begin{equation} \label{eq:tra_bbl_Drho}
 	\nabla_\sigma \rho / \rho = \alpha \,\nabla_\sigma T + \beta   \,\nabla_\sigma S
 \end{equation} 
@@ -978,5 +978,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Advective bottom boundary layer  (\protect\np{nn\_bbl\_adv}\forcode{ = 1..2})}
-\label{TRA_bbl_adv}
+\label{subsec:TRA_bbl_adv}
 
 \sgacomment{"downsloping flow" has been replaced by "downslope flow" in the following
@@ -986,5 +986,5 @@
 \begin{figure}[!t] 	\begin{center}
 \includegraphics[width=0.7\textwidth]{Fig_BBL_adv}
-\caption{ 	\protect\label{Fig_bbl}  
+\caption{ 	\protect\label{fig:bbl}  
 Advective/diffusive Bottom Boundary Layer. The BBL parameterisation is 
 activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$. 
@@ -1011,5 +1011,5 @@
 
 \np{nn\_bbl\_adv}\forcode{ = 1} : the downslope velocity is chosen to be the Eulerian
-ocean velocity just above the topographic step (see black arrow in Fig.\ref{Fig_bbl}) 
+ocean velocity just above the topographic step (see black arrow in \autoref{fig:bbl}) 
 \citep{Beckmann_Doscher1997}. It is a \textit{conditional advection}, that is, advection
 is allowed only if dense water overlies less dense water on the slope ($i.e.$ 
@@ -1021,7 +1021,7 @@
 The advection is allowed only  if dense water overlies less dense water on the slope ($i.e.$ 
 $\nabla_\sigma \rho  \cdot  \nabla H<0$). For example, the resulting transport of the 
-downslope flow, here in the $i$-direction (Fig.\ref{Fig_bbl}), is simply given by the 
+downslope flow, here in the $i$-direction (\autoref{fig:bbl}), is simply given by the 
 following expression:
-\begin{equation} \label{Eq_bbl_Utr}
+\begin{equation} \label{eq:bbl_Utr}
  u^{tr}_{bbl} = \gamma \, g \frac{\Delta \rho}{\rho_o}  e_{1u} \; min \left( {e_{3u}}_{kup},{e_{3u}}_{kdwn} \right)
 \end{equation}
@@ -1039,21 +1039,21 @@
 water at intermediate depths. The entrainment is replaced by the vertical mixing 
 implicit in the advection scheme. Let us consider as an example the 
-case displayed in Fig.\ref{Fig_bbl} where the density at level $(i,kup)$ is 
+case displayed in \autoref{fig:bbl} where the density at level $(i,kup)$ is 
 larger than the one at level $(i,kdwn)$. The advective BBL scheme
 modifies the tracer time tendency of the ocean cells near the 
-topographic step by the downslope flow \eqref{Eq_bbl_dw}, 
-the horizontal \eqref{Eq_bbl_hor}  and the upward \eqref{Eq_bbl_up} 
+topographic step by the downslope flow \autoref{eq:bbl_dw}, 
+the horizontal \autoref{eq:bbl_hor} and the upward \autoref{eq:bbl_up} 
 return flows as follows: 
 \begin{align} 
 \partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw}
-                                     +  \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}}  \left( T^{sh}_{kup} - T^{do}_{kdw} \right)  \label{Eq_bbl_dw} \\
+                                     +  \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}}  \left( T^{sh}_{kup} - T^{do}_{kdw} \right)  \label{eq:bbl_dw} \\
 %
 \partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} 
-				   + \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}}   \left( T^{do}_{kup} - T^{sh}_{kup} \right)   \label{Eq_bbl_hor} \\
+				   + \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}}   \left( T^{do}_{kup} - T^{sh}_{kup} \right)   \label{eq:bbl_hor} \\
 %
 \intertext{and for $k =kdw-1,\;..., \; kup$ :} 
 %
 \partial_t T^{do}_{k} &\equiv \partial_t S^{do}_{k}
-				   + \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   \left( T^{do}_{k+1} - T^{sh}_{k} \right)   \label{Eq_bbl_up}
+				   + \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   \left( T^{do}_{k+1} - T^{sh}_{k} \right)   \label{eq:bbl_up}
 \end{align}
 where $b_t$ is the $T$-cell volume. 
@@ -1067,5 +1067,5 @@
 % ================================================================
 \section{Tracer damping (\protect\mdl{tradmp})}
-\label{TRA_dmp}
+\label{sec:TRA_dmp}
 %--------------------------------------------namtra_dmp-------------------------------------------------
 \forfile{../namelists/namtra_dmp}
@@ -1074,5 +1074,5 @@
 In some applications it can be useful to add a Newtonian damping term 
 into the temperature and salinity equations:
-\begin{equation} \label{Eq_tra_dmp}
+\begin{equation} \label{eq:tra_dmp}
 \begin{split}
  \frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right)  \\
@@ -1087,8 +1087,8 @@
 in \textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are 
 correctly set  ($i.e.$ that $T_o$ and $S_o$ are provided in input files and read 
-using \mdl{fldread}, see \S\ref{SBC_fldread}). 
+using \mdl{fldread}, see \autoref{subsec:SBC_fldread}). 
 The restoring coefficient $\gamma$ is a three-dimensional array read in during the \rou{tra\_dmp\_init} routine. The file name is specified by the namelist variable \np{cn\_resto}. The DMP\_TOOLS tool is provided to allow users to generate the netcdf file.
 
-The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)} 
+The two main cases in which \autoref{eq:tra_dmp} is used are \textit{(a)} 
 the specification of the boundary conditions along artificial walls of a 
 limited domain basin and \textit{(b)} the computation of the velocity 
@@ -1151,5 +1151,5 @@
 % ================================================================
 \section{Tracer time evolution (\protect\mdl{tranxt})}
-\label{TRA_nxt}
+\label{sec:TRA_nxt}
 %--------------------------------------------namdom-----------------------------------------------------
 \forfile{../namelists/namdom}
@@ -1159,6 +1159,6 @@
 The general framework for tracer time stepping is a modified leap-frog scheme 
 \citep{Leclair_Madec_OM09}, $i.e.$ a three level centred time scheme associated 
-with a Asselin time filter (cf. \S\ref{STP_mLF}):
-\begin{equation} \label{Eq_tra_nxt}
+with a Asselin time filter (cf. \autoref{sec:STP_mLF}):
+\begin{equation} \label{eq:tra_nxt}
 \begin{aligned}
 (e_{3t}T)^{t+\rdt} &= (e_{3t}T)_f^{t-\rdt} &+ 2 \, \rdt  \,e_{3t}^t\ \text{RHS}^t &	\\
@@ -1174,5 +1174,5 @@
 $\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter). 
 Its default value is \np{rn\_atfp}\forcode{ = 10.e-3}. Note that the forcing correction term in the filter
-is not applied in linear free surface (\jp{lk\_vvl}\forcode{ = .false.}) (see \S\ref{TRA_sbc}.
+is not applied in linear free surface (\jp{lk\_vvl}\forcode{ = .false.}) (see \autoref{subsec:TRA_sbc}.
 Not also that in constant volume case, the time stepping is performed on $T$, 
 not on its content, $e_{3t}T$.
@@ -1189,5 +1189,5 @@
 % ================================================================
 \section{Equation of state (\protect\mdl{eosbn2}) }
-\label{TRA_eosbn2}
+\label{sec:TRA_eosbn2}
 %--------------------------------------------nameos-----------------------------------------------------
 \forfile{../namelists/nameos}
@@ -1198,5 +1198,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Equation of seawater (\protect\np{nn\_eos}\forcode{ = -1..1})}
-\label{TRA_eos}
+\label{subsec:TRA_eos}
 
 The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship 
@@ -1272,5 +1272,5 @@
 and \textit{practical} salinity.
 S-EOS takes the following expression:
-\begin{equation} \label{Eq_tra_S-EOS}
+\begin{equation} \label{eq:tra_S-EOS}
 \begin{split}
   d_a(T,S,z)  =  ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a  \\
@@ -1280,5 +1280,5 @@
 \end{split}
 \end{equation} 
-where the computer name of the coefficients as well as their standard value are given in \ref{Tab_SEOS}.
+where the computer name of the coefficients as well as their standard value are given in \autoref{tab:SEOS}.
 In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing 
 the associated coefficients. 
@@ -1303,5 +1303,5 @@
 $\mu_2$     & \np{rn\_mu2} 	& 1.1090 $10^{-5}$ &  thermobaric coeff. in S   	      \\ \hline
 \end{tabular}
-\caption{ \protect\label{Tab_SEOS}
+\caption{ \protect\label{tab:SEOS}
 Standard value of S-EOS coefficients. }
 \end{center}
@@ -1314,5 +1314,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos}\forcode{ = 0..2})}
-\label{TRA_bn2}
+\label{subsec:TRA_bn2}
 
 An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a}
@@ -1323,5 +1323,5 @@
  (pressure in decibar being approximated by the depth in meters). The expression for $N^2$ 
  is given by: 
-\begin{equation} \label{Eq_tra_bn2}
+\begin{equation} \label{eq:tra_bn2}
 N^2 =	\frac{g}{e_{3w}} \left(   \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T]   \right)
 \end{equation} 
@@ -1336,8 +1336,8 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Freezing point of seawater}
-\label{TRA_fzp}
+\label{subsec:TRA_fzp}
 
 The freezing point of seawater is a function of salinity and pressure \citep{UNESCO1983}:
-\begin{equation} \label{Eq_tra_eos_fzp}
+\begin{equation} \label{eq:tra_eos_fzp}
    \begin{split}
 T_f (S,p) = \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} 
@@ -1347,7 +1347,7 @@
 \end{equation}
 
-\eqref{Eq_tra_eos_fzp} is only used to compute the potential freezing point of 
+\autoref{eq:tra_eos_fzp} is only used to compute the potential freezing point of 
 sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent 
-terms in \eqref{Eq_tra_eos_fzp} (last term) have been dropped. The freezing
+terms in \autoref{eq:tra_eos_fzp} (last term) have been dropped. The freezing
 point is computed through \textit{eos\_fzp}, a \textsc{Fortran} function that can be found 
 in \mdl{eosbn2}.  
@@ -1358,5 +1358,5 @@
 % -------------------------------------------------------------------------------------------------------------
 %\subsection{Potential Energy anomalies}
-%\label{TRA_bn2}
+%\label{subsec:TRA_bn2}
 
 %    =====>>>>> TO BE written
@@ -1368,5 +1368,5 @@
 % ================================================================
 \section{Horizontal derivative in \textit{zps}-coordinate (\protect\mdl{zpshde})}
-\label{TRA_zpshde}
+\label{sec:TRA_zpshde}
 
 \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, 
@@ -1382,5 +1382,5 @@
 Before taking horizontal gradients between the tracers next to the bottom, a linear 
 interpolation in the vertical is used to approximate the deeper tracer as if it actually 
-lived at the depth of the shallower tracer point (Fig.~\ref{Fig_Partial_step_scheme}). 
+lived at the depth of the shallower tracer point (\autoref{fig:Partial_step_scheme}). 
 For example, for temperature in the $i$-direction the needed interpolated 
 temperature, $\widetilde{T}$, is:
@@ -1389,7 +1389,7 @@
 \begin{figure}[!p] 	 \begin{center}
 \includegraphics[width=0.9\textwidth]{Partial_step_scheme}
-\caption{ 	\protect\label{Fig_Partial_step_scheme} 
+\caption{ 	\protect\label{fig:Partial_step_scheme} 
 Discretisation of the horizontal difference and average of tracers in the $z$-partial 
-step coordinate (\np{ln\_zps}\forcode{ = .true.}) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. 
+step coordinate (\protect\np{ln\_zps}\forcode{ = .true.}) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. 
 A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value 
 at the depth of the shallower tracer point of the two adjacent bottom $T$-points. 
@@ -1409,5 +1409,5 @@
 and the resulting forms for the horizontal difference and the horizontal average 
 value of $T$ at a $U$-point are: 
-\begin{equation} \label{Eq_zps_hde}
+\begin{equation} \label{eq:zps_hde}
 \begin{aligned}
  \delta _{i+1/2} T= 	\begin{cases}
@@ -1432,6 +1432,6 @@
 of density, we compute $\widetilde{\rho }$ from the interpolated values of $T$ 
 and $S$, and the pressure at a $u$-point (in the equation of state pressure is 
-approximated by depth, see \S\ref{TRA_eos} ) : 
-\begin{equation} \label{Eq_zps_hde_rho}
+approximated by depth, see \autoref{subsec:TRA_eos} ) : 
+\begin{equation} \label{eq:zps_hde_rho}
 \widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u }) 
 \quad \text{where }\  z_u = \min \left( {z_T^{i+1} ,z_T^i } \right)
@@ -1441,9 +1441,9 @@
 thus pressure) is highly non-linear with a true equation of state and thus is badly 
 approximated with a linear interpolation. This approximation is used to compute 
-both the horizontal pressure gradient (\S\ref{DYN_hpg}) and the slopes of neutral 
-surfaces (\S\ref{LDF_slp})
+both the horizontal pressure gradient (\autoref{sec:DYN_hpg}) and the slopes of neutral 
+surfaces (\autoref{sec:LDF_slp})
 
 Note that in almost all the advection schemes presented in this Chapter, both 
-averaging and differencing operators appear. Yet \eqref{Eq_zps_hde} has not 
+averaging and differencing operators appear. Yet \autoref{eq:zps_hde} has not 
 been used in these schemes: in contrast to diffusion and pressure gradient 
 computations, no correction for partial steps is applied for advection. The main 
