Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_ZDF.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_ZDF.tex	(revision 9393)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_ZDF.tex	(revision 9407)
@@ -5,5 +5,5 @@
 % ================================================================
 \chapter{Vertical Ocean Physics (ZDF)}
-\label{ZDF}
+\label{chap:ZDF}
 \minitoc
 
@@ -19,19 +19,19 @@
 % ================================================================
 \section{Vertical mixing}
-\label{ZDF_zdf}
+\label{sec:ZDF_zdf}
 
 The discrete form of the ocean subgrid scale physics has been presented in 
-\S\ref{TRA_zdf} and \S\ref{DYN_zdf}. At the surface and bottom boundaries, 
+\autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}. At the surface and bottom boundaries, 
 the turbulent fluxes of momentum, heat and salt have to be defined. At the 
-surface they are prescribed from the surface forcing (see Chap.~\ref{SBC}), 
+surface they are prescribed from the surface forcing (see \autoref{chap:SBC}), 
 while at the bottom they are set to zero for heat and salt, unless a geothermal 
 flux forcing is prescribed as a bottom boundary condition ($i.e.$ \key{trabbl} 
-defined, see \S\ref{TRA_bbc}), and specified through a bottom friction 
-parameterisation for momentum (see \S\ref{ZDF_bfr}).
+defined, see \autoref{subsec:TRA_bbc}), and specified through a bottom friction 
+parameterisation for momentum (see \autoref{sec:ZDF_bfr}).
 
 In this section we briefly discuss the various choices offered to compute 
 the vertical eddy viscosity and diffusivity coefficients, $A_u^{vm}$ , 
 $A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$- 
-points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These 
+points, respectively (see \autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}). These 
 coefficients can be assumed to be either constant, or a function of the local 
 Richardson number, or computed from a turbulent closure model (either 
@@ -44,5 +44,5 @@
 (namelist parameter \np{ln\_zdfexp}\forcode{ = .true.}) or a backward time stepping 
 scheme (\np{ln\_zdfexp}\forcode{ = .false.}) depending on the magnitude of the mixing 
-coefficients, and thus of the formulation used (see \S\ref{STP}).
+coefficients, and thus of the formulation used (see \autoref{chap:STP}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -50,5 +50,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Constant (\protect\key{zdfcst})}
-\label{ZDF_cst}
+\label{subsec:ZDF_cst}
 %--------------------------------------------namzdf---------------------------------------------------------
 \forfile{../namelists/namzdf}
@@ -75,5 +75,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Richardson number dependent (\protect\key{zdfric})}
-\label{ZDF_ric}
+\label{subsec:ZDF_ric}
 
 %--------------------------------------------namric---------------------------------------------------------
@@ -91,5 +91,5 @@
 ratio of stratification to vertical shear). Following \citet{Pacanowski_Philander_JPO81}, the following 
 formulation has been implemented:
-\begin{equation} \label{Eq_zdfric}
+\begin{equation} \label{eq:zdfric}
    \left\{      \begin{aligned}
          A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT}       \\
@@ -98,7 +98,7 @@
 \end{equation}
 where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson 
-number, $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 
+number, $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 
 $A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the 
-constant case (see \S\ref{ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ 
+constant case (see \autoref{subsec:ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ 
 is the maximum value that can be reached by the coefficient when $Ri\leq 0$, 
 $a=5$ and $n=2$. The last three values can be modified by setting the 
@@ -133,5 +133,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{TKE turbulent closure scheme (\protect\key{zdftke})}
-\label{ZDF_tke}
+\label{subsec:ZDF_tke}
 
 %--------------------------------------------namzdf_tke--------------------------------------------------
@@ -150,5 +150,5 @@
 $\bar{e}$ through vertical shear, its destruction through stratification, its vertical 
 diffusion, and its dissipation of \citet{Kolmogorov1942} type:
-\begin{equation} \label{Eq_zdftke_e}
+\begin{equation} \label{eq:zdftke_e}
 \frac{\partial \bar{e}}{\partial t} = 
 \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
@@ -159,5 +159,5 @@
 - c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }
 \end{equation}
-\begin{equation} \label{Eq_zdftke_kz}
+\begin{equation} \label{eq:zdftke_kz}
    \begin{split}
          K_m &= C_k\  l_k\  \sqrt {\bar{e}\; }  	\\
@@ -165,5 +165,5 @@
    \end{split}
 \end{equation}
-where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 
+where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 
 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 
 $P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity 
@@ -173,5 +173,5 @@
 $P_{rt}$ can be set to unity or, following \citet{Blanke1993}, be a function 
 of the local Richardson number, $R_i$:
-\begin{align*} \label{Eq_prt}
+\begin{align*} \label{eq:prt}
 P_{rt} = \begin{cases}
                     \ \ \ 1 &      \text{if $\ R_i \leq 0.2$} 	\\
@@ -187,5 +187,5 @@
 namelist parameter. The default value of $e_{bb}$ is 3.75. \citep{Gaspar1990}), 
 however a much larger value can be used when taking into account the 
-surface wave breaking (see below Eq. \eqref{ZDF_Esbc}). 
+surface wave breaking (see below Eq. \autoref{eq:ZDF_Esbc}). 
 The bottom value of TKE is assumed to be equal to the value of the level just above. 
 The time integration of the $\bar{e}$ equation may formally lead to negative values 
@@ -199,5 +199,5 @@
 instabilities associated with too weak vertical diffusion. They must be 
 specified at least larger than the molecular values, and are set through 
-\np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \S\ref{ZDF_cst}).
+\np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \autoref{subsec:ZDF_cst}).
 
 \subsubsection{Turbulent length scale}
@@ -207,5 +207,5 @@
 parameter. The first two are based on the following first order approximation 
 \citep{Blanke1993}:
-\begin{equation} \label{Eq_tke_mxl0_1}
+\begin{equation} \label{eq:tke_mxl0_1}
 l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N
 \end{equation}
@@ -219,10 +219,10 @@
 \np{nn\_mxl}\forcode{ = 2..3} cases, which add an extra assumption concerning the vertical 
 gradient of the computed length scale. So, the length scales are first evaluated 
-as in \eqref{Eq_tke_mxl0_1} and then bounded such that:
-\begin{equation} \label{Eq_tke_mxl_constraint}
+as in \autoref{eq:tke_mxl0_1} and then bounded such that:
+\begin{equation} \label{eq:tke_mxl_constraint}
 \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1
 \qquad \text{with }\  l =  l_k = l_\epsilon
 \end{equation}
-\eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length 
+\autoref{eq:tke_mxl_constraint} means that the vertical variations of the length 
 scale cannot be larger than the variations of depth. It provides a better 
 approximation of the \citet{Gaspar1990} formulation while being much less 
@@ -230,5 +230,5 @@
 by the distance to the surface or to the ocean bottom but also by the distance 
 to a strongly stratified portion of the water column such as the thermocline 
-(Fig.~\ref{Fig_mixing_length}). In order to impose the \eqref{Eq_tke_mxl_constraint} 
+(\autoref{fig:mixing_length}). In order to impose the \autoref{eq:tke_mxl_constraint} 
 constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, 
 the upward and downward length scales, and evaluate the dissipation and 
@@ -237,10 +237,10 @@
 \begin{figure}[!t] \begin{center}
 \includegraphics[width=1.00\textwidth]{Fig_mixing_length}
-\caption{ \protect\label{Fig_mixing_length} 
+\caption{ \protect\label{fig:mixing_length} 
 Illustration of the mixing length computation. }
 \end{center}  
 \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{equation} \label{Eq_tke_mxl2}
+\begin{equation} \label{eq:tke_mxl2}
 \begin{aligned}
  l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \;  \right)
@@ -250,5 +250,5 @@
 \end{aligned}
 \end{equation}
-where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1}, 
+where $l^{(k)}$ is computed using \autoref{eq:tke_mxl0_1}, 
 $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
 
@@ -257,5 +257,5 @@
 \np{nn\_mxl}\forcode{ = 3} case, the dissipation and mixing turbulent length scales are give 
 as in \citet{Gaspar1990}:
-\begin{equation} \label{Eq_tke_mxl_gaspar}
+\begin{equation} \label{eq:tke_mxl_gaspar}
 \begin{aligned}
 & l_k          = \sqrt{\  l_{up} \ \ l_{dwn}\ }  	\\
@@ -282,5 +282,5 @@
 
 Following \citet{Craig_Banner_JPO94}, the boundary condition on surface TKE value is :
-\begin{equation}  \label{ZDF_Esbc}
+\begin{equation}  \label{eq:ZDF_Esbc}
 \bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
 \end{equation}
@@ -289,5 +289,5 @@
 younger waves \citep{Mellor_Blumberg_JPO04}. 
 The boundary condition on the turbulent length scale follows the Charnock's relation:
-\begin{equation} \label{ZDF_Lsbc}
+\begin{equation} \label{eq:ZDF_Lsbc}
 l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o}
 \end{equation}
@@ -297,5 +297,5 @@
 As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, 
 with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds 
-to $\alpha_{CB} = 100$. Further setting  \np{ln\_mxl0} to true applies \eqref{ZDF_Lsbc} 
+to $\alpha_{CB} = 100$. Further setting  \np{ln\_mxl0} to true applies \autoref{eq:ZDF_Lsbc} 
 as surface boundary condition on length scale, with $\beta$ hard coded to the Stacey's value.
 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) 
@@ -316,5 +316,5 @@
 The parameterization, tuned against large-eddy simulation, includes the whole effect
 of LC in an extra source terms of TKE, $P_{LC}$.
-The presence of $P_{LC}$ in \eqref{Eq_zdftke_e}, the TKE equation, is controlled 
+The presence of $P_{LC}$ in \autoref{eq:zdftke_e}, the TKE equation, is controlled 
 by setting \np{ln\_lc} to \forcode{.true.} in the namtke namelist.
  
@@ -368,6 +368,6 @@
 When using this parameterization ($i.e.$ when \np{nn\_etau}\forcode{ = 1}), the TKE input to the ocean ($S$) 
 imposed by the winds in the form of near-inertial oscillations, swell and waves is parameterized 
-by \eqref{ZDF_Esbc} the standard TKE surface boundary condition, plus a depth depend one given by:
-\begin{equation}  \label{ZDF_Ehtau}
+by \autoref{eq:ZDF_Esbc} the standard TKE surface boundary condition, plus a depth depend one given by:
+\begin{equation}  \label{eq:ZDF_Ehtau}
 S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau} 
 \end{equation}
@@ -385,5 +385,5 @@
 
 Note that two other option existe, \np{nn\_etau}\forcode{ = 2..3}. They correspond to applying 
-\eqref{ZDF_Ehtau} only at the base of the mixed layer, or to using the high frequency part 
+\autoref{eq:ZDF_Ehtau} only at the base of the mixed layer, or to using the high frequency part 
 of the stress to evaluate the fraction of TKE that penetrate the ocean. 
 Those two options are obsolescent features introduced for test purposes.
@@ -406,10 +406,10 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{TKE discretization considerations (\protect\key{zdftke})}
-\label{ZDF_tke_ene}
+\label{subsec:ZDF_tke_ene}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]   \begin{center}
 \includegraphics[width=1.00\textwidth]{Fig_ZDF_TKE_time_scheme}
-\caption{ \protect\label{Fig_TKE_time_scheme} 
+\caption{ \protect\label{fig:TKE_time_scheme} 
 Illustration of the TKE time integration and its links to the momentum and tracer time integration. }
 \end{center}  
@@ -418,15 +418,15 @@
 
 The production of turbulence by vertical shear (the first term of the right hand side 
-of \eqref{Eq_zdftke_e}) should balance the loss of kinetic energy associated with
-the vertical momentum diffusion (first line in \eqref{Eq_PE_zdf}). To do so a special care 
+of \autoref{eq:zdftke_e}) should balance the loss of kinetic energy associated with
+the vertical momentum diffusion (first line in \autoref{eq:PE_zdf}). To do so a special care 
 have to be taken for both the time and space discretization of the TKE equation 
 \citep{Burchard_OM02,Marsaleix_al_OM08}.
 
-Let us first address the time stepping issue. Fig.~\ref{Fig_TKE_time_scheme} shows 
+Let us first address the time stepping issue. \autoref{fig:TKE_time_scheme} shows 
 how the two-level Leap-Frog time stepping of the momentum and tracer equations interplays 
 with the one-level forward time stepping of TKE equation. With this framework, the total loss 
 of kinetic energy (in 1D for the demonstration) due to the vertical momentum diffusion is 
 obtained by multiplying this quantity by $u^t$ and summing the result vertically:   
-\begin{equation} \label{Eq_energ1}
+\begin{equation} \label{eq:energ1}
 \begin{split}
 \int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\
@@ -436,10 +436,10 @@
 \end{equation}
 Here, the vertical diffusion of momentum is discretized backward in time 
-with a coefficient, $K_m$, known at time $t$ (Fig.~\ref{Fig_TKE_time_scheme}), 
-as it is required when using the TKE scheme (see \S\ref{STP_forward_imp}). 
-The first term of the right hand side of \eqref{Eq_energ1} represents the kinetic energy 
+with a coefficient, $K_m$, known at time $t$ (\autoref{fig:TKE_time_scheme}), 
+as it is required when using the TKE scheme (see \autoref{sec:STP_forward_imp}). 
+The first term of the right hand side of \autoref{eq:energ1} represents the kinetic energy 
 transfer at the surface (atmospheric forcing) and at the bottom (friction effect). 
 The second term is always negative. It is the dissipation rate of kinetic energy, 
-and thus minus the shear production rate of $\bar{e}$. \eqref{Eq_energ1} 
+and thus minus the shear production rate of $\bar{e}$. \autoref{eq:energ1} 
 implies that, to be energetically consistent, the production rate of $\bar{e}$ 
 used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as 
@@ -448,10 +448,10 @@
 
 A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification 
-(second term of the right hand side of \eqref{Eq_zdftke_e}). This term 
+(second term of the right hand side of \autoref{eq:zdftke_e}). This term 
 must balance the input of potential energy resulting from vertical mixing. 
 The rate of change of potential energy (in 1D for the demonstration) due vertical 
 mixing is obtained by multiplying vertical density diffusion 
 tendency by $g\,z$ and and summing the result vertically:
-\begin{equation} \label{Eq_energ2}
+\begin{equation} \label{eq:energ2}
 \begin{split}
 \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\
@@ -463,15 +463,15 @@
 \end{equation}
 where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. 
-The first term of the right hand side of \eqref{Eq_energ2}  is always zero 
+The first term of the right hand side of \autoref{eq:energ2}  is always zero 
 because there is no diffusive flux through the ocean surface and bottom). 
 The second term is minus the destruction rate of  $\bar{e}$ due to stratification. 
-Therefore \eqref{Eq_energ1} implies that, to be energetically consistent, the product 
-${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \eqref{Eq_zdftke_e}, the TKE equation.
+Therefore \autoref{eq:energ1} implies that, to be energetically consistent, the product 
+${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:zdftke_e}, the TKE equation.
 
 Let us now address the space discretization issue. 
 The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity 
 components are in the centre of the side faces of a $t$-box in staggered C-grid 
-(Fig.\ref{Fig_cell}). A space averaging is thus required to obtain the shear TKE production term.
-By redoing the \eqref{Eq_energ1} in the 3D case, it can be shown that the product of 
+(\autoref{fig:cell}). A space averaging is thus required to obtain the shear TKE production term.
+By redoing the \autoref{eq:energ1} in the 3D case, it can be shown that the product of 
 eddy coefficient by the shear at $t$ and $t-\rdt$ must be performed prior to the averaging.
 Furthermore, the possible time variation of $e_3$ (\key{vvl} case) have to be taken into 
@@ -480,5 +480,5 @@
 The above energetic considerations leads to 
 the following final discrete form for the TKE equation:
-\begin{equation} \label{Eq_zdftke_ene}
+\begin{equation} \label{eq:zdftke_ene}
 \begin{split}
 \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv  
@@ -497,10 +497,10 @@
 \end{split}
 \end{equation}
-where the last two terms in \eqref{Eq_zdftke_ene} (vertical diffusion and Kolmogorov dissipation) 
-are time stepped using a backward scheme (see\S\ref{STP_forward_imp}). 
+where the last two terms in \autoref{eq:zdftke_ene} (vertical diffusion and Kolmogorov dissipation) 
+are time stepped using a backward scheme (see\autoref{sec:STP_forward_imp}). 
 Note that the Kolmogorov term has been linearized in time in order to render 
 the implicit computation possible. The restart of the TKE scheme 
 requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as they all appear in 
-the right hand side of \eqref{Eq_zdftke_ene}. For the latter, it is in fact 
+the right hand side of \autoref{eq:zdftke_ene}. For the latter, it is in fact 
 the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 
 
@@ -509,5 +509,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{GLS: Generic Length Scale (\protect\key{zdfgls})}
-\label{ZDF_gls}
+\label{subsec:ZDF_gls}
 
 %--------------------------------------------namzdf_gls---------------------------------------------------------
@@ -519,10 +519,10 @@
 for the generic length scale, $\psi$ \citep{Umlauf_Burchard_JMS03, Umlauf_Burchard_CSR05}. 
 This later variable is defined as : $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 
-where the triplet $(p, m, n)$ value given in Tab.\ref{Tab_GLS} allows to recover 
+where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:GLS} allows to recover 
 a number of well-known turbulent closures ($k$-$kl$ \citep{Mellor_Yamada_1982}, 
 $k$-$\epsilon$ \citep{Rodi_1987}, $k$-$\omega$ \citep{Wilcox_1988} 
 among others \citep{Umlauf_Burchard_JMS03,Kantha_Carniel_CSR05}). 
 The GLS scheme is given by the following set of equations:
-\begin{equation} \label{Eq_zdfgls_e}
+\begin{equation} \label{eq:zdfgls_e}
 \frac{\partial \bar{e}}{\partial t} = 
 \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
@@ -533,5 +533,5 @@
 \end{equation}
 
-\begin{equation} \label{Eq_zdfgls_psi}
+\begin{equation} \label{eq:zdfgls_psi}
    \begin{split}
 \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{
@@ -544,5 +544,5 @@
 \end{equation}
 
-\begin{equation} \label{Eq_zdfgls_kz}
+\begin{equation} \label{eq:zdfgls_kz}
    \begin{split}
          K_m    &= C_{\mu} \ \sqrt {\bar{e}} \ l         \\
@@ -551,12 +551,12 @@
 \end{equation}
 
-\begin{equation} \label{Eq_zdfgls_eps}
+\begin{equation} \label{eq:zdfgls_eps}
 {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \;
 \end{equation}
-where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}) 
+where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) 
 and $\epsilon$ the dissipation rate. 
 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) 
 depends of the choice of the turbulence model. Four different turbulent models are pre-defined 
-(Tab.\ref{Tab_GLS}). They are made available through the \np{nn\_clo} namelist parameter. 
+(Tab.\autoref{tab:GLS}). They are made available through the \np{nn\_clo} namelist parameter. 
 
 %--------------------------------------------------TABLE--------------------------------------------------
@@ -579,5 +579,5 @@
 \hline
 \end{tabular}
-\caption{   \protect\label{Tab_GLS} 
+\caption{   \protect\label{tab:GLS} 
 Set of predefined GLS parameters, or equivalently predefined turbulence models available 
 with \protect\key{zdfgls} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls} .}
@@ -596,5 +596,5 @@
 As for TKE closure , the wave effect on the mixing is considered when \np{ln\_crban}\forcode{ = .true.}
 \citep{Craig_Banner_JPO94, Mellor_Blumberg_JPO04}. The \np{rn\_crban} namelist parameter 
-is $\alpha_{CB}$ in \eqref{ZDF_Esbc} and \np{rn\_charn} provides the value of $\beta$ in \eqref{ZDF_Lsbc}. 
+is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and \np{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}. 
 
 The $\psi$ equation is known to fail in stably stratified flows, and for this reason 
@@ -609,5 +609,5 @@
 
 The time and space discretization of the GLS equations follows the same energetic 
-consideration as for the TKE case described in \S\ref{ZDF_tke_ene}  \citep{Burchard_OM02}. 
+consideration as for the TKE case described in \autoref{subsec:ZDF_tke_ene}  \citep{Burchard_OM02}. 
 Examples of performance of the 4 turbulent closure scheme can be found in \citet{Warner_al_OM05}.
 
@@ -616,5 +616,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{OSM: OSMOSIS boundary layer scheme (\protect\key{zdfosm})}
-\label{ZDF_osm}
+\label{subsec:ZDF_osm}
 
 %--------------------------------------------namzdf_osm---------------------------------------------------------
@@ -628,5 +628,5 @@
 % ================================================================
 \section{Convection}
-\label{ZDF_conv}
+\label{sec:ZDF_conv}
 
 %--------------------------------------------namzdf--------------------------------------------------------
@@ -648,5 +648,5 @@
 \subsection[Non-penetrative convective adjmt (\protect\np{ln\_tranpc}\forcode{ = .true.})]
 				{Non-penetrative convective adjustment (\protect\np{ln\_tranpc}\forcode{ = .true.})}
-\label{ZDF_npc}
+\label{subsec:ZDF_npc}
 
 %--------------------------------------------namzdf--------------------------------------------------------
@@ -657,5 +657,5 @@
 \begin{figure}[!htb]  	\begin{center}
 \includegraphics[width=0.90\textwidth]{Fig_npc}
-\caption{  \protect\label{Fig_npc} 
+\caption{  \protect\label{fig:npc} 
 Example of an unstable density profile treated by the non penetrative 
 convective adjustment algorithm. $1^{st}$ step: the initial profile is checked from 
@@ -677,5 +677,5 @@
 column has \textit{exactly} the density of the water just below) \citep{Madec_al_JPO91}. 
 The associated algorithm is an iterative process used in the following way 
-(Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is 
+(\autoref{fig:npc}): starting from the top of the ocean, the first instability is 
 found. Assume in the following that the instability is located between levels 
 $k$ and $k+1$. The temperature and salinity in the two levels are 
@@ -714,5 +714,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Enhanced vertical diffusion (\protect\np{ln\_zdfevd}\forcode{ = .true.})}
-\label{ZDF_evd}
+\label{subsec:ZDF_evd}
 
 %--------------------------------------------namzdf--------------------------------------------------------
@@ -739,5 +739,5 @@
 Note that the stability test is performed on both \textit{before} and \textit{now} 
 values of $N^2$. This removes a potential source of divergence of odd and
-even time step in a leapfrog environment \citep{Leclair_PhD2010} (see \S\ref{STP_mLF}).
+even time step in a leapfrog environment \citep{Leclair_PhD2010} (see \autoref{sec:STP_mLF}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -745,15 +745,15 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection[Turbulent closure scheme (\protect\key{zdf}\{tke,gls,osm\})]{Turbulent Closure Scheme (\protect\key{zdftke}, \protect\key{zdfgls} or \protect\key{zdfosm})}
-\label{ZDF_tcs}
-
-The turbulent closure scheme presented in \S\ref{ZDF_tke} and \S\ref{ZDF_gls} 
+\label{subsec:ZDF_tcs}
+
+The turbulent closure scheme presented in \autoref{subsec:ZDF_tke} and \autoref{subsec:ZDF_gls} 
 (\key{zdftke} or \key{zdftke} is defined) in theory solves the problem of statically 
 unstable density profiles. In such a case, the term corresponding to the 
-destruction of turbulent kinetic energy through stratification in \eqref{Eq_zdftke_e} 
-or \eqref{Eq_zdfgls_e} becomes a source term, since $N^2$ is negative. 
+destruction of turbulent kinetic energy through stratification in \autoref{eq:zdftke_e} 
+or \autoref{eq:zdfgls_e} becomes a source term, since $N^2$ is negative. 
 It results in large values of $A_T^{vT}$ and  $A_T^{vT}$, and also the four neighbouring 
 $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1})$. These large values 
 restore the static stability of the water column in a way similar to that of the 
-enhanced vertical diffusion parameterisation (\S\ref{ZDF_evd}). However, 
+enhanced vertical diffusion parameterisation (\autoref{subsec:ZDF_evd}). However, 
 in the vicinity of the sea surface (first ocean layer), the eddy coefficients 
 computed by the turbulent closure scheme do not usually exceed $10^{-2}m.s^{-1}$, 
@@ -772,5 +772,5 @@
 % ================================================================
 \section{Double diffusion mixing (\protect\key{zdfddm})}
-\label{ZDF_ddm}
+\label{sec:ZDF_ddm}
 
 %-------------------------------------------namzdf_ddm-------------------------------------------------
@@ -789,5 +789,5 @@
 
 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients 
-\begin{align*} % \label{Eq_zdfddm_Kz}
+\begin{align*} % \label{eq:zdfddm_Kz}
     &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT}	\\
     &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS}
@@ -797,12 +797,12 @@
 mixing depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, 
 where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline 
-contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt 
+contraction (see \autoref{subsec:TRA_eos}). To represent mixing of $S$ and $T$ by salt 
 fingering, we adopt the diapycnal diffusivities suggested by Schmitt (1981):
-\begin{align} \label{Eq_zdfddm_f}
+\begin{align} \label{eq:zdfddm_f}
 A_f^{vS} &= 	\begin{cases}
 	\frac{A^{\ast v}}{1+(R_\rho / R_c)^n   } &\text{if  $R_\rho > 1$ and $N^2>0$ } \\
 	0 				  					    &\text{otherwise} 
 				\end{cases}   
-\\ 		    \label{Eq_zdfddm_f_T}
+\\ 		    \label{eq:zdfddm_f_T}
 A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 
 \end{align}
@@ -811,5 +811,5 @@
 \begin{figure}[!t]   \begin{center}
 \includegraphics[width=0.99\textwidth]{Fig_zdfddm}
-\caption{  \protect\label{Fig_zdfddm}
+\caption{  \protect\label{fig:zdfddm}
 From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$ 
 and $A_f^{vS}$ for temperature and salt in regions of salt fingering. Heavy 
@@ -822,5 +822,5 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-The factor 0.7 in \eqref{Eq_zdfddm_f_T} reflects the measured ratio 
+The factor 0.7 in \autoref{eq:zdfddm_f_T} reflects the measured ratio 
 $\alpha F_T /\beta F_S \approx  0.7$ of buoyancy flux of heat to buoyancy 
 flux of salt ($e.g.$, \citet{McDougall_Taylor_JMR84}). Following  \citet{Merryfield1999}, 
@@ -828,5 +828,5 @@
 
 To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by Federov (1988) is used: 
-\begin{align} 	\label{Eq_zdfddm_d}
+\begin{align} 	\label{eq:zdfddm_d}
 A_d^{vT} &= 	\begin{cases}
 	1.3635 \, \exp{\left( 4.6\, \exp{ \left[  -0.54\,( R_{\rho}^{-1} - 1 )  \right] }    \right)}
@@ -834,5 +834,5 @@
 	0 								&\text{otherwise} 
 				\end{cases}   
-\\   			\label{Eq_zdfddm_d_S}
+\\   			\label{eq:zdfddm_d_S}
 A_d^{vS} &= 	\begin{cases}
 	A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right)
@@ -844,6 +844,6 @@
 \end{align}
 
-The dependencies of \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d_S} on $R_\rho$ 
-are illustrated in Fig.~\ref{Fig_zdfddm}. Implementing this requires computing 
+The dependencies of \autoref{eq:zdfddm_f} to \autoref{eq:zdfddm_d_S} on $R_\rho$ 
+are illustrated in \autoref{fig:zdfddm}. Implementing this requires computing 
 $R_\rho$ at each grid point on every time step. This is done in \mdl{eosbn2} at the 
 same time as $N^2$ is computed. This avoids duplication in the computation of 
@@ -854,5 +854,5 @@
 % ================================================================
 \section{Bottom and top friction (\protect\mdl{zdfbfr})}
-\label{ZDF_bfr}
+\label{sec:ZDF_bfr}
 
 %--------------------------------------------nambfr--------------------------------------------------------
@@ -870,5 +870,5 @@
 flux (bottom friction) enter the equations as a condition on the vertical 
 diffusive flux. For the bottom boundary layer, one has:
-\begin{equation} \label{Eq_zdfbfr_flux}
+\begin{equation} \label{eq:zdfbfr_flux}
 A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U}
 \end{equation}
@@ -886,5 +886,5 @@
 as a body force over the depth of the top or bottom model layer. To illustrate 
 this, consider the equation for $u$ at $k$, the last ocean level:
-\begin{equation} \label{Eq_zdfbfr_flux2}
+\begin{equation} \label{eq:zdfbfr_flux2}
 \frac{\partial u_k}{\partial t} = \frac{1}{e_{3u}} \left[ \frac{A_{uw}^{vm}}{e_{3uw}} \delta_{k+1/2}\;[u] - {\cal F}^u_h \right] \approx - \frac{{\cal F}^u_{h}}{e_{3u}}
 \end{equation}
@@ -907,5 +907,5 @@
 These coefficients are computed in \mdl{zdfbfr} and generally take the form 
 $c_b^{\textbf U}$ where:
-\begin{equation} \label{Eq_zdfbfr_bdef}
+\begin{equation} \label{eq:zdfbfr_bdef}
 \frac{\partial {\textbf U_h}}{\partial t} = 
   - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b
@@ -917,5 +917,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Linear bottom friction (\protect\np{nn\_botfr}\forcode{ = 0..1})}
-\label{ZDF_bfr_linear}
+\label{subsec:ZDF_bfr_linear}
 
 The linear bottom friction parameterisation (including the special case 
@@ -923,5 +923,5 @@
 is proportional to the interior velocity (i.e. the velocity of the last 
 model level):
-\begin{equation} \label{Eq_zdfbfr_linear}
+\begin{equation} \label{eq:zdfbfr_linear}
 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b
 \end{equation}
@@ -941,6 +941,6 @@
 
 For the linear friction case the coefficients defined in the general 
-expression \eqref{Eq_zdfbfr_bdef} are: 
-\begin{equation} \label{Eq_zdfbfr_linbfr_b}
+expression \autoref{eq:zdfbfr_bdef} are: 
+\begin{equation} \label{eq:zdfbfr_linbfr_b}
 \begin{split}
  c_b^u &= - r\\
@@ -961,9 +961,9 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Non-linear bottom friction (\protect\np{nn\_botfr}\forcode{ = 2})}
-\label{ZDF_bfr_nonlinear}
+\label{subsec:ZDF_bfr_nonlinear}
 
 The non-linear bottom friction parameterisation assumes that the bottom 
 friction is quadratic: 
-\begin{equation} \label{Eq_zdfbfr_nonlinear}
+\begin{equation} \label{eq:zdfbfr_nonlinear}
 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 
 }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b 
@@ -983,5 +983,5 @@
 For the non-linear friction case the terms
 computed in \mdl{zdfbfr}  are: 
-\begin{equation} \label{Eq_zdfbfr_nonlinbfr}
+\begin{equation} \label{eq:zdfbfr_nonlinbfr}
 \begin{split}
  c_b^u &= - \; C_D\;\left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^{1/2}\\
@@ -1003,5 +1003,5 @@
 \subsection[Log-layer btm frict enhncmnt (\protect\np{nn\_botfr}\forcode{ = 2}, \protect\np{ln\_loglayer}\forcode{ = .true.})]
 				{Log-layer bottom friction enhancement (\protect\np{nn\_botfr}\forcode{ = 2}, \protect\np{ln\_loglayer}\forcode{ = .true.})}
-\label{ZDF_bfr_loglayer}
+\label{subsec:ZDF_bfr_loglayer}
 
 In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally
@@ -1033,11 +1033,11 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Bottom friction stability considerations}
-\label{ZDF_bfr_stability}
+\label{subsec:ZDF_bfr_stability}
 
 Some care needs to exercised over the choice of parameters to ensure that the
 implementation of bottom friction does not induce numerical instability. For 
-the purposes of stability analysis, an approximation to \eqref{Eq_zdfbfr_flux2}
+the purposes of stability analysis, an approximation to \autoref{eq:zdfbfr_flux2}
 is:
-\begin{equation} \label{Eqn_bfrstab}
+\begin{equation} \label{eq:Eqn_bfrstab}
 \begin{split}
  \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\
@@ -1050,5 +1050,5 @@
  |\Delta u| < \;|u| 
 \end{equation}
-\noindent which, using \eqref{Eqn_bfrstab}, gives:
+\noindent which, using \autoref{eq:Eqn_bfrstab}, gives:
 \begin{equation}
 r\frac{2\rdt}{e_{3u}} < 1 \qquad  \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\
@@ -1075,5 +1075,5 @@
 
 Limits on the bottom friction coefficient are not imposed if the user has elected to
-handle the bottom friction implicitly (see \S\ref{ZDF_bfr_imp}). The number of potential
+handle the bottom friction implicitly (see \autoref{subsec:ZDF_bfr_imp}). The number of potential
 breaches of the explicit stability criterion are still reported for information purposes.
 
@@ -1082,5 +1082,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Implicit bottom friction (\protect\np{ln\_bfrimp}\forcode{ = .true.})}
-\label{ZDF_bfr_imp}
+\label{subsec:ZDF_bfr_imp}
 
 An optional implicit form of bottom friction has been implemented to improve
@@ -1093,5 +1093,5 @@
 bottom boundary condition is implemented implicitly.
 
-\begin{equation} \label{Eq_dynzdf_bfr}
+\begin{equation} \label{eq:dynzdf_bfr}
 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk}
 	 = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}}
@@ -1112,9 +1112,9 @@
 following:
 
-\begin{equation} \label{Eq_dynspg_ts_bfr1}
+\begin{equation} \label{eq:dynspg_ts_bfr1}
 \frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b}
 \left(\textbf{U}_{med}-\textbf{U}^{m-1}\right)
 \end{equation}
-\begin{equation} \label{Eq_dynspg_ts_bfr2}
+\begin{equation} \label{eq:dynspg_ts_bfr2}
 \frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+
 \left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)-
@@ -1136,5 +1136,5 @@
 \subsection[Bottom friction w/ split-explicit time splitting (\protect\np{ln\_bfrimp})]
 				{Bottom friction with split-explicit time splitting (\protect\np{ln\_bfrimp})}
-\label{ZDF_bfr_ts}
+\label{subsec:ZDF_bfr_ts}
 
 When calculating the momentum trend due to bottom friction in \mdl{dynbfr}, the
@@ -1175,9 +1175,9 @@
 the barotropic component which uses the unrestricted value of the coefficient. However, if the
 limiting is thought to be having a major effect (a more likely prospect in coastal and shelf seas
-applications) then the fully implicit form of the bottom friction should be used (see \S\ref{ZDF_bfr_imp} ) 
+applications) then the fully implicit form of the bottom friction should be used (see \autoref{subsec:ZDF_bfr_imp} ) 
 which can be selected by setting \np{ln\_bfrimp} $=$ \forcode{.true.}.
 
 Otherwise, the implicit formulation takes the form:
-\begin{equation} \label{Eq_zdfbfr_implicitts}
+\begin{equation} \label{eq:zdfbfr_implicitts}
  \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ]  
 \end{equation}
@@ -1193,5 +1193,5 @@
 % ================================================================
 \section{Tidal mixing (\protect\key{zdftmx})}
-\label{ZDF_tmx}
+\label{sec:ZDF_tmx}
 
 %--------------------------------------------namzdf_tmx--------------------------------------------------
@@ -1204,5 +1204,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Bottom intensified tidal mixing}
-\label{ZDF_tmx_bottom}
+\label{subsec:ZDF_tmx_bottom}
 
 Options are defined through the  \ngn{namzdf\_tmx} namelist variables.
@@ -1213,9 +1213,9 @@
 $A^{vT}_{tides}$ is expressed as a function of $E(x,y)$, the energy transfer from barotropic 
 tides to baroclinic tides : 
-\begin{equation} \label{Eq_Ktides}
+\begin{equation} \label{eq:Ktides}
 A^{vT}_{tides} =  q \,\Gamma \,\frac{ E(x,y) \, F(z) }{ \rho \, N^2 }
 \end{equation}
 where $\Gamma$ is the mixing efficiency, $N$ the Brunt-Vais\"{a}l\"{a} frequency 
-(see \S\ref{TRA_bn2}), $\rho$ the density, $q$ the tidal dissipation efficiency, 
+(see \autoref{subsec:TRA_bn2}), $\rho$ the density, $q$ the tidal dissipation efficiency, 
 and $F(z)$ the vertical structure function. 
 
@@ -1230,5 +1230,5 @@
 It is implemented as a simple exponential decaying upward away from the bottom, 
 with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter, with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04}, 
-\begin{equation} \label{Eq_Fz}
+\begin{equation} \label{eq:Fz}
 F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) }
 \end{equation}
@@ -1241,5 +1241,5 @@
 usually set to $10^{-8} s^{-2}$. These bounds are usually rarely encountered.
 
-The internal wave energy map, $E(x, y)$ in \eqref{Eq_Ktides}, is derived 
+The internal wave energy map, $E(x, y)$ in \autoref{eq:Ktides}, is derived 
 from a barotropic model of the tides utilizing a parameterization of the 
 conversion of barotropic tidal energy into internal waves. 
@@ -1250,5 +1250,5 @@
 the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}.
 This model provides the dissipation associated with internal wave energy for the M2 and K1 
-tides component (Fig.~\ref{Fig_ZDF_M2_K1_tmx}). The S2 dissipation is simply approximated
+tides component (\autoref{fig:ZDF_M2_K1_tmx}). The S2 dissipation is simply approximated
 as being $1/4$ of the M2 one. The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$. 
 Its global mean value is $1.1$ TW, in agreement with independent estimates 
@@ -1258,5 +1258,5 @@
 \begin{figure}[!t] 	\begin{center}
 \includegraphics[width=0.90\textwidth]{Fig_ZDF_M2_K1_tmx}
-\caption{  \protect\label{Fig_ZDF_M2_K1_tmx} 
+\caption{  \protect\label{fig:ZDF_M2_K1_tmx} 
 (a) M2 and (b) K1 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$). }
 \end{center}  	\end{figure}
@@ -1267,5 +1267,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Indonesian area specific treatment (\protect\np{ln\_zdftmx\_itf})}
-\label{ZDF_tmx_itf}
+\label{subsec:ZDF_tmx_itf}
 
 When the Indonesian Through Flow (ITF) area is included in the model domain,
@@ -1294,5 +1294,5 @@
 proportional to $N^2$ below the core of the thermocline and to $N$ above. 
 The resulting $F(z)$ is:
-\begin{equation} \label{Eq_Fz_itf}
+\begin{equation} \label{eq:Fz_itf}
 F(i,j,k) \sim     \left\{ \begin{aligned}
 \frac{q\,\Gamma E(i,j) } {\rho N \, \int N     dz}    \qquad \text{when $\partial_z N < 0$} \\
@@ -1315,5 +1315,5 @@
 % ================================================================
 \section{Internal wave-driven mixing (\protect\key{zdftmx\_new})}
-\label{ZDF_tmx_new}
+\label{sec:ZDF_tmx_new}
 
 %--------------------------------------------namzdf_tmx_new------------------------------------------
@@ -1325,5 +1325,5 @@
 A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed, 
 and the resulting diffusivity is obtained as 
-\begin{equation} \label{Eq_Kwave}
+\begin{equation} \label{eq:Kwave}
 A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 }
 \end{equation}
