Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex	(revision 11543)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex	(revision 11544)
@@ -18,5 +18,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \section{Upstream Biased Scheme (UBS) (\protect\np{ln\_traadv\_ubs}\forcode{ = .true.})}
-\label{sec:TRA_adv_ubs}
+\label{sec:ALGOS_tra_adv_ubs}
 
 The UBS advection scheme is an upstream biased third order scheme based on
@@ -25,5 +25,5 @@
 For example, in the $i$-direction:
 \begin{equation}
-  \label{eq:tra_adv_ubs2}
+  \label{eq:ALGOS_tra_adv_ubs2}
   \tau_u^{ubs} = \left\{
 	 \begin{aligned}
@@ -35,5 +35,5 @@
 or equivalently, the advective flux is
 \begin{equation}
-  \label{eq:tra_adv_ubs2}
+  \label{eq:ALGOS_tra_adv_ubs2}
   U_{i+1/2} \ \tau_u^{ubs}
   =U_{i+1/2} \ \overline{ T_i - \frac{1}{6}\,\tau"_i }^{\,i+1/2}
@@ -85,5 +85,5 @@
 NB 3: It is straight forward to rewrite \autoref{eq:TRA_adv_ubs} as follows:
 \begin{equation}
-  \label{eq:tra_adv_ubs2}
+  \label{eq:ALGOS_tra_adv_ubs2}
   \tau_u^{ubs} = \left\{
 	 \begin{aligned}
@@ -95,5 +95,5 @@
 or equivalently
 \begin{equation}
-  \label{eq:tra_adv_ubs2}
+  \label{eq:ALGOS_tra_adv_ubs2}
   \begin{split}
     e_{2u} e_{3u}\,u_{i+1/2} \ \tau_u^{ubs}
@@ -112,5 +112,5 @@
 laplacian diffusion:
 \begin{equation}
-  \label{eq:tra_ldf_lap}
+  \label{eq:ALGOS_tra_ldf_lap}
   \begin{split}
     D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\;  e_{3T} } &\left[ {\quad \delta_i
@@ -125,5 +125,5 @@
 bilaplacian:
 \begin{equation}
-  \label{eq:tra_ldf_lap}
+  \label{eq:ALGOS_tra_ldf_lap}
   \begin{split}
     D_T^{lT} =&-\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
@@ -138,5 +138,5 @@
 it comes:
 \begin{equation}
-  \label{eq:tra_ldf_lap}
+  \label{eq:ALGOS_tra_ldf_lap}
   \begin{split}
     D_T^{lT} =&-\frac{1}{12}\,\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
@@ -149,5 +149,5 @@
 if the velocity is uniform (\ie\ $|u|=cst$) then the diffusive flux is
 \begin{equation}
-  \label{eq:tra_ldf_lap}
+  \label{eq:ALGOS_tra_ldf_lap}
   \begin{split}
     F_u^{lT} = - \frac{1}{12}
@@ -161,5 +161,5 @@
 
 \begin{equation}
-  \label{eq:tra_adv_ubs2}
+  \label{eq:ALGOS_tra_adv_ubs2}
   \begin{split}
     F_u^{lT} &= - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
@@ -171,5 +171,5 @@
 sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$):
 \begin{equation}
-  \label{eq:tra_adv_ubs2}
+  \label{eq:ALGOS_tra_adv_ubs2}
   \begin{split}
     F_u^{lT} &= - \frac{1}{12} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{e_{1T}^3\,|u|}{e_{1T}e_{2T}\,e_{3T}}\,\delta_i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta_{i+1/2}[\tau] \right] \right]
@@ -180,5 +180,5 @@
 sol 2 coefficient at u-point: split $|u|$ into $\sqrt{|u|}$ and $e_{1T}$ into $\sqrt{e_{1u}}$
 \begin{equation}
-  \label{eq:tra_adv_ubs2}
+  \label{eq:ALGOS_tra_adv_ubs2}
   \begin{split}
     F_u^{lT} &= - \frac{1}{12} {e_{1u}}^1 \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{1}{e_{2T}\,e_{3T}}\,\delta_i \left[ \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u} }{e_{1u} } \delta_{i+1/2}[\tau] \right] \right] \\
@@ -192,5 +192,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \section{Leapfrog energetic}
-\label{sec:LF}
+\label{sec:ALGOS_LF}
 
 We adopt the following semi-discrete notation for time derivative.
@@ -198,5 +198,5 @@
 the time derivation and averaging operators at the mid time step are:
 \[
-  % \label{eq:dt_mt}
+  % \label{eq:ALGOS_dt_mt}
   \begin{split}
     \delta_{t+\rdt/2} [q]     &=  \  \ \,   q^{t+\rdt}  - q^{t}		\\
@@ -210,5 +210,5 @@
 The Leap-frog time stepping given by \autoref{eq:DOM_nxt} can be defined as:
 \[
-  % \label{eq:LF}
+  % \label{eq:ALGOS_LF}
   \frac{\partial q}{\partial t}
   \equiv \frac{1}{\rdt} \overline{ \delta_{t+\rdt/2}[q]}^{\,t}
@@ -220,5 +220,5 @@
 As such it respects the quadratic invariant in integral forms, \ie\ the following continuous property,
 \[
-  % \label{eq:Energy}
+  % \label{eq:ALGOS_Energy}
   \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}
   =\int_{t_0}^{t_1} {\frac{1}{2}\, \frac{\partial q^2}{\partial t} \;dt}
@@ -278,5 +278,5 @@
 For example in the (\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows:
 \begin{equation}
-  \label{eq:Gf_triads}
+  \label{eq:ALGOS_Gf_triads}
   _i^k \mathbb{T}_{i_p}^{k_p} (T)
   = \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k  	\left(
@@ -291,5 +291,5 @@
 and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad:
 \begin{equation}
-  \label{eq:Gf_slopes}
+  \label{eq:ALGOS_Gf_slopes}
   _i^k \mathbb{R}_{i_p}^{k_p}
   =\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} \ \frac
@@ -297,11 +297,11 @@
   {\left(\alpha / \beta \right)_i^k  \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }
 \end{equation}
-Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of
+Note that in \autoref{eq:ALGOS_Gf_slopes} we use the ratio $\alpha / \beta$ instead of
 multiplying the temperature derivative by $\alpha$ and the salinity derivative by $\beta$.
 This is more efficient as the ratio $\alpha / \beta$ can to be evaluated directly.
 
-Note that in \autoref{eq:Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of ${b_{uw}}_{\,i+i_p}^{\,k+k_p}$.
+Note that in \autoref{eq:ALGOS_Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of ${b_{uw}}_{\,i+i_p}^{\,k+k_p}$.
 This choice has been motivated by the decrease of tracer variance and
-the presence of partial cell at the ocean bottom (see \autoref{apdx:Gf_operator}).
+the presence of partial cell at the ocean bottom (see \autoref{subsec:ALGOS_Gf_operator}).
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -310,5 +310,5 @@
     \includegraphics[width=\textwidth]{Fig_ISO_triad}
     \caption{
-      \protect\label{fig:ISO_triad}
+      \protect\label{fig:ALGOS_ISO_triad}
       Triads used in the Griffies's like iso-neutral diffision scheme for
       $u$-component (upper panel) and $w$-component (lower panel).
@@ -321,5 +321,5 @@
 They take the following expression:
 \begin{flalign*}
-  % \label{eq:Gf_fluxes}
+  % \label{eq:ALGOS_Gf_fluxes}
   \begin{split}
     {_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)
@@ -334,5 +334,5 @@
 the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:TRIADS_ISO_triad}):
 \begin{flalign}
-  \label{eq:iso_flux}
+  \label{eq:ALGOS_iso_flux}
   \textbf{F}_{iso}(T)
   &\equiv  \sum_{\substack{i_p,\,k_p}}
@@ -364,5 +364,5 @@
 the divergence of the sum of all the four triad fluxes:
 \begin{equation}
-  \label{eq:Gf_operator}
+  \label{eq:ALGOS_Gf_operator}
   D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{
     \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
@@ -377,5 +377,5 @@
   the limit of flat iso-neutral direction:
   \[
-    % \label{eq:Gf_property1a}
+    % \label{eq:ALGOS_Gf_property1a}
     D_l^T = \frac{1}{b_T}  \ \delta_{i}
     \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right]
@@ -401,5 +401,5 @@
   The iso-neutral flux of locally referenced potential density is zero, \ie
   \begin{align*}
-    % \label{eq:Gf_property2}
+    % \label{eq:ALGOS_Gf_property2}
     \begin{matrix}
       &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} (\rho)}
@@ -411,11 +411,11 @@
     \end{matrix}
   \end{align*}
-  This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$ and
-  the definition of the triads' slopes \autoref{eq:Gf_slopes}.
+  This result is trivially obtained using the \autoref{eq:ALGOS_Gf_triads} applied to $T$ and $S$ and
+  the definition of the triads' slopes \autoref{eq:ALGOS_Gf_slopes}.
 
 \item[$\bullet$ conservation of tracer]
   The iso-neutral diffusion term conserve the total tracer content, \ie
   \[
-    % \label{eq:Gf_property1}
+    % \label{eq:ALGOS_Gf_property1}
     \sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0
   \]
@@ -425,8 +425,8 @@
   The iso-neutral diffusion term does not increase the total tracer variance, \ie
   \[
-    % \label{eq:Gf_property1}
+    % \label{eq:ALGOS_Gf_property1}
     \sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0
   \]
-The property is demonstrated in the \autoref{apdx:Gf_operator}.
+The property is demonstrated in the \autoref{subsec:ALGOS_Gf_operator}.
 It is a key property for a diffusion term.
 It means that the operator is also a dissipation term,
@@ -438,5 +438,5 @@
   The iso-neutral diffusion operator is self-adjoint, \ie
   \[
-    % \label{eq:Gf_property1}
+    % \label{eq:ALGOS_Gf_property1}
     \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}
   \]
@@ -444,5 +444,5 @@
 We just have to apply the same routine.
 This properties can be demonstrated quite easily in a similar way the "non increase of tracer variance" property
-has been proved (see \autoref{apdx:Gf_operator}).
+has been proved (see \autoref{apdx:ALGOS_Gf_operator}).
 \end{description}
 
@@ -462,5 +462,5 @@
 The eddy induced velocity is given by:
 \begin{equation}
-  \label{eq:eiv_v}
+  \label{eq:ALGOS_eiv_v}
   \begin{split}
     u^* & = - \frac{1}{e_2\,e_{3}}          \;\partial_k \left( e_2 \, A_e \; r_i  \right)
@@ -487,5 +487,5 @@
 % give here the expression using the triads. It is different from the one given in \autoref{eq:LDF_eiv}
 % see just below a copy of this equation:
-%\begin{equation} \label{eq:ldfeiv}
+%\begin{equation} \label{eq:ALGOS_ldfeiv}
 %\begin{split}
 % u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\
@@ -495,5 +495,5 @@
 %\end{equation}
 \[
-  % \label{eq:eiv_vd}
+  % \label{eq:ALGOS_eiv_vd}
   \textbf{F}_{eiv}^T   \equiv   \left(
     \begin{aligned}
@@ -540,5 +540,5 @@
 we end up with the skew form of the eddy induced advective fluxes:
 \begin{equation}
-  \label{eq:eiv_skew_continuous}
+  \label{eq:ALGOS_eiv_skew_continuous}
   \textbf{F}_{eiv}^T =
   \begin{pmatrix}
@@ -548,11 +548,11 @@
 \end{equation}
 The tendency associated with eddy induced velocity is then simply the divergence of
-the \autoref{eq:eiv_skew_continuous} fluxes.
+the \autoref{eq:ALGOS_eiv_skew_continuous} fluxes.
 It naturally conserves the tracer content, as it is expressed in flux form and,
 as the advective form, it preserves the tracer variance.
-Another interesting property of \autoref{eq:eiv_skew_continuous} form is that when $A=A_e$,
+Another interesting property of \autoref{eq:ALGOS_eiv_skew_continuous} form is that when $A=A_e$,
 a simplification occurs in the sum of the iso-neutral diffusion and eddy induced velocity terms:
 \begin{flalign*}
-  % \label{eq:eiv_skew+eiv_continuous}
+  % \label{eq:ALGOS_eiv_skew+eiv_continuous}
   \textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &=
   \begin{pmatrix}
@@ -576,10 +576,10 @@
 This property has been used to reduce the computational time \citep{griffies_JPO98},
 but it is not of practical use as usually $A \neq A_e$.
-Nevertheless this property can be used to choose a discret form of \autoref{eq:eiv_skew_continuous} which
-is consistent with the iso-neutral operator \autoref{eq:Gf_operator}.
-Using the slopes \autoref{eq:Gf_slopes} and defining $A_e$ at $T$-point(\ie\ as $A$,
+Nevertheless this property can be used to choose a discret form of \autoref{eq:ALGOS_eiv_skew_continuous} which
+is consistent with the iso-neutral operator \autoref{eq:ALGOS_Gf_operator}.
+Using the slopes \autoref{eq:ALGOS_Gf_slopes} and defining $A_e$ at $T$-point(\ie\ as $A$,
 the eddy diffusivity coefficient), the resulting discret form is given by:
 \begin{equation}
-  \label{eq:eiv_skew}
+  \label{eq:ALGOS_eiv_skew}
   \textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left(
     \begin{aligned}
@@ -593,5 +593,5 @@
   \right)
 \end{equation}
-Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells.
+Note that \autoref{eq:ALGOS_eiv_skew} is valid in $z$-coordinate with or without partial cells.
 In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces must be added to
 $\mathbb{R}$ for the discret form to be exact.
@@ -599,5 +599,5 @@
 Such a choice of discretisation is consistent with the iso-neutral operator as
 it uses the same definition for the slopes.
-It also ensures the conservation of the tracer variance (see Appendix \autoref{apdx:eiv_skew}),
+It also ensures the conservation of the tracer variance (see \autoref{subsec:ALGOS_eiv_skew}),
 \ie\ it does not include a diffusive component but is a "pure" advection term.
 
@@ -607,5 +607,5 @@
 % ================================================================
 \subsection{Discrete invariants of the iso-neutral diffrusion}
-\label{subsec:Gf_operator}
+\label{subsec:ALGOS_Gf_operator}
 
 Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane.
@@ -702,5 +702,5 @@
   \intertext{
   Then outing in factor the triad in each of the four terms of the summation and
-  substituting the triads by their expression given in \autoref{eq:Gf_triads}.
+  substituting the triads by their expression given in \autoref{eq:ALGOS_Gf_triads}.
   It becomes:
   }
@@ -774,5 +774,5 @@
 % ================================================================
 \subsection{Discrete invariants of the skew flux formulation}
-\label{subsec:eiv_skew}
+\label{subsec:ALGOS_eiv_skew}
 
 Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.
@@ -784,5 +784,5 @@
   \int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0
 \end{align*}
-The discrete form of its left hand side is obtained using \autoref{eq:eiv_skew}
+The discrete form of its left hand side is obtained using \autoref{eq:ALGOS_eiv_skew}
 \begin{align*}
   \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex	(revision 11543)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex	(revision 11544)
@@ -57,5 +57,5 @@
 %===============================================================
 \section{Warm layer model}
-\label{sec:warm_layer_sec}
+\label{sec:DIU_warm_layer_sec}
 %===============================================================
 
@@ -65,9 +65,9 @@
 \frac{\partial{\Delta T_{\mathrm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p
 \nu}-\frac{(\nu+1)ku^*_{w}f(L_a)\Delta T}{D_T\Phi\!\left(\frac{D_T}{L}\right)} \mbox{,}
-\label{eq:ecmwf1} \\
-L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:ecmwf2}
+\label{eq:DIU_ecmwf1} \\
+L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:DIU_ecmwf2}
 \end{align}
 where $\Delta T_{\mathrm{wl}}$ is the temperature difference between the top of the warm layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal.
-In equation (\autoref{eq:ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion coefficient of water,
+In equation (\autoref{eq:DIU_ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion coefficient of water,
 $\kappa=0.4$ is von K\'{a}rm\'{a}n's constant, $c_p$ is the heat capacity at constant pressure of sea water,
 $\rho_w$ is the water density, and $L$ is the Monin-Obukhov length.
@@ -79,14 +79,14 @@
 the relationship $u^*_{w} = u_{10}\sqrt{\frac{C_d\rho_a}{\rho_w}}$, where $C_d$ is the drag coefficient,
 and $\rho_a$ is the density of air.
-The symbol $Q$ in equation (\autoref{eq:ecmwf1}) is the instantaneous total thermal energy flux into
+The symbol $Q$ in equation (\autoref{eq:DIU_ecmwf1}) is the instantaneous total thermal energy flux into
 the diurnal layer, \ie
 \[
   Q = Q_{\mathrm{sol}} + Q_{\mathrm{lw}} + Q_{\mathrm{h}}\mbox{,}
-  % \label{eq:e_flux_eqn}
+  % \label{eq:DIU_e_flux_eqn}
 \]
 where $Q_{\mathrm{h}}$ is the sensible and latent heat flux, $Q_{\mathrm{lw}}$ is the long wave flux,
 and $Q_{\mathrm{sol}}$ is the solar flux absorbed within the diurnal warm layer.
 For $Q_{\mathrm{sol}}$ the 9 term representation of \citet{gentemann.minnett.ea_JGR09} is used.
-In equation \autoref{eq:ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$,
+In equation \autoref{eq:DIU_ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$,
 where $L_a=0.3$\footnote{
   This is a global average value, more accurately $L_a$ could be computed as $L_a=(u^*_{w}/u_s)^{\frac{1}{2}}$,
@@ -99,11 +99,11 @@
 4\zeta^2}{1+3\zeta+0.25\zeta^2} &(\zeta \ge 0) \\
                                     (1 - 16\zeta)^{-\frac{1}{2}} & (\zeta < 0) \mbox{,}
-                                    \end{array} \right. \label{eq:stab_func_eqn}
+                                    \end{array} \right. \label{eq:DIU_stab_func_eqn}
 \end{equation}
-where $\zeta=\frac{D_T}{L}$.  It is clear that the first derivative of (\autoref{eq:stab_func_eqn}),
-and thus of (\autoref{eq:ecmwf1}), is discontinuous at $\zeta=0$ (\ie\ $Q\rightarrow0$ in
-equation (\autoref{eq:ecmwf2})).
+where $\zeta=\frac{D_T}{L}$.  It is clear that the first derivative of (\autoref{eq:DIU_stab_func_eqn}),
+and thus of (\autoref{eq:DIU_ecmwf1}), is discontinuous at $\zeta=0$ (\ie\ $Q\rightarrow0$ in
+equation (\autoref{eq:DIU_ecmwf2})).
 
-The two terms on the right hand side of (\autoref{eq:ecmwf1}) represent different processes.
+The two terms on the right hand side of (\autoref{eq:DIU_ecmwf1}) represent different processes.
 The first term is simply the diabatic heating or cooling of the diurnal warm layer due to
 thermal energy fluxes into and out of the layer.
@@ -114,5 +114,5 @@
 
 \section{Cool skin model}
-\label{sec:cool_skin_sec}
+\label{sec:DIU_cool_skin_sec}
 
 %===============================================================
@@ -121,5 +121,5 @@
 As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\mathrm{cs}}$ becomes
 \[
-  % \label{eq:sunders_eqn}
+  % \label{eq:DIU_sunders_eqn}
   \Delta T_{\mathrm{cs}}=\frac{Q_{\mathrm{ns}}\delta}{k_t} \mbox{,}
 \]
@@ -128,5 +128,5 @@
 $\delta$ is the thickness of the skin layer and is given by
 \begin{equation}
-\label{eq:sunders_thick_eqn}
+\label{eq:DIU_sunders_thick_eqn}
 \delta=\frac{\lambda \mu}{u^*_{w}} \mbox{,}
 \end{equation}
@@ -134,10 +134,10 @@
 \citet{saunders_JAS67} suggested varied between 5 and 10.
 
-The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{artale.iudicone.ea_JGR02},
+The value of $\lambda$ used in equation (\autoref{eq:DIU_sunders_thick_eqn}) is that of \citet{artale.iudicone.ea_JGR02},
 which is shown in \citet{tu.tsuang_GRL05} to outperform a number of other parametrisations at
 both low and high wind speeds.
 Specifically,
 \[
-  % \label{eq:artale_lambda_eqn}
+  % \label{eq:DIU_artale_lambda_eqn}
   \lambda = \frac{ 8.64\times10^4 u^*_{w} k_t }{ \rho c_p h \mu \gamma }\mbox{,}
 \]
@@ -145,5 +145,5 @@
 $\gamma$ is a dimensionless function of wind speed $u$:
 \[
-  % \label{eq:artale_gamma_eqn}
+  % \label{eq:DIU_artale_gamma_eqn}
   \gamma =
   \begin{cases}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex	(revision 11543)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex	(revision 11544)
@@ -30,5 +30,5 @@
 As detailed in \cite{brankart_OM13}, the stochastic formulation of the equation of state can be written as:
 \begin{equation}
-  \label{eq:eos_sto}
+  \label{eq:STO_eos_sto}
   \rho = \frac{1}{2} \sum_{i=1}^m\{ \rho[T+\Delta T_i,S+\Delta S_i,p_o(z)] + \rho[T-\Delta T_i,S-\Delta S_i,p_o(z)] \}
 \end{equation}
@@ -37,5 +37,5 @@
 the scalar product of the respective local T/S gradients with random walks $\mathbf{\xi}$:
 \begin{equation}
-  \label{eq:sto_pert}
+  \label{eq:STO_sto_pert}
   \Delta T_i = \mathbf{\xi}_i \cdot \nabla T \qquad \hbox{and} \qquad \Delta S_i = \mathbf{\xi}_i \cdot \nabla S
 \end{equation}
@@ -59,5 +59,5 @@
 
 \begin{equation}
-  \label{eq:autoreg}
+  \label{eq:STO_autoreg}
   \xi^{(i)}_{k+1} = a^{(i)} \xi^{(i)}_k + b^{(i)} w^{(i)} + c^{(i)}
 \end{equation}
@@ -74,5 +74,5 @@
 
   \[
-    % \label{eq:ord1}
+    % \label{eq:STO_ord1}
     \left\{
       \begin{array}{l}
@@ -91,5 +91,5 @@
 
   \begin{equation}
-    \label{eq:ord2}
+    \label{eq:STO_ord2}
     \left\{
       \begin{array}{l}
@@ -107,7 +107,7 @@
 \noindent
 In this way, higher order processes can be easily generated recursively using the same piece of code implementing
-\autoref{eq:autoreg}, and using successive processes from order $0$ to~$n-1$ as~$w^{(i)}$.
-The parameters in \autoref{eq:ord2} are computed so that this recursive application of
-\autoref{eq:autoreg} leads to processes with the required standard deviation and correlation timescale,
+\autoref{eq:STO_autoreg}, and using successive processes from order $0$ to~$n-1$ as~$w^{(i)}$.
+The parameters in \autoref{eq:STO_ord2} are computed so that this recursive application of
+\autoref{eq:STO_autoreg} leads to processes with the required standard deviation and correlation timescale,
 with the additional condition that the $n-1$ first derivatives of the autocorrelation function are equal to
 zero at~$t=0$, so that the resulting processes become smoother and smoother as $n$ increases.
@@ -135,11 +135,11 @@
 
 \mdl{stopts} : stochastic parametrisation associated with the non-linearity of the equation of
- seawater, implementing \autoref{eq:sto_pert} so as specifics in the equation of state
- implementing \autoref{eq:eos_sto}.
+ seawater, implementing \autoref{eq:STO_sto_pert} so as specifics in the equation of state
+ implementing \autoref{eq:STO_eos_sto}.
 % \end{description}
 
 The \mdl{stopar} module includes three public routines called in the model:
 
-(\rou{sto\_par}) is a direct implementation of \autoref{eq:autoreg},
+(\rou{sto\_par}) is a direct implementation of \autoref{eq:STO_autoreg},
 applied at each model grid point (in 2D or 3D), and called at each model time step ($k$) to
 update every autoregressive process ($i=1,\ldots,m$).
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex	(revision 11543)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex	(revision 11544)
@@ -7,5 +7,6 @@
 % ================================================================
 \chapter{Invariants of the Primitive Equations}
-\label{chap:Invariant}
+\label{chap:CONS}
+
 \chaptertoc
 
@@ -45,5 +46,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \section{Conservation properties on ocean dynamics}
-\label{sec:Invariant_dyn}
+\label{sec:CONS_Invariant_dyn}
 
 The non linear term of the momentum equations has been split into a vorticity term,
@@ -63,5 +64,5 @@
 The continuous formulation of the vorticity term satisfies following integral constraints:
 \[
-  % \label{eq:vor_vorticity}
+  % \label{eq:CONS_vor_vorticity}
   \int_D {{\textbf {k}}\cdot \frac{1}{e_3 }\nabla \times \left( {\varsigma
         \;{\mathrm {\mathbf k}}\times {\textbf {U}}_h } \right)\;dv} =0
@@ -69,5 +70,5 @@
 
 \[
-  % \label{eq:vor_enstrophy}
+  % \label{eq:CONS_vor_enstrophy}
   if\quad \chi =0\quad \quad \int\limits_D {\varsigma \;{\textbf{k}}\cdot
     \frac{1}{e_3 }\nabla \times \left( {\varsigma {\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} =-\int\limits_D {\frac{1}{2}\varsigma ^2\,\chi \;dv}
@@ -76,5 +77,5 @@
 
 \[
-  % \label{eq:vor_energy}
+  % \label{eq:CONS_vor_energy}
   \int_D {{\textbf{U}}_h \times \left( {\varsigma \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} =0
 \]
@@ -88,5 +89,5 @@
 Using the symmetry or anti-symmetry properties of the operators (Eqs II.1.10 and 11),
 it can be shown that the scheme (II.2.11) satisfies (II.4.1b) but not (II.4.1c),
-while scheme (II.2.12) satisfies (II.4.1c) but not (II.4.1b) (see appendix C). 
+while scheme (II.2.12) satisfies (II.4.1c) but not (II.4.1b) (see appendix C).
 Note that the enstrophy conserving scheme on total vorticity has been chosen as the standard discrete form of
 the vorticity term.
@@ -102,12 +103,12 @@
 the horizontal gradient of horizontal kinetic energy:
 
-\begin{equation} \label{eq:keg_zad}
-\int_D {{\textbf{U}}_h \cdot \nabla _h \left( {1/2\;{\textbf{U}}_h ^2} \right)\;dv} =-\int_D {{\textbf{U}}_h \cdot \frac{w}{e_3 }\;\frac{\partial 
+\begin{equation} \label{eq:CONS_keg_zad}
+\int_D {{\textbf{U}}_h \cdot \nabla _h \left( {1/2\;{\textbf{U}}_h ^2} \right)\;dv} =-\int_D {{\textbf{U}}_h \cdot \frac{w}{e_3 }\;\frac{\partial
 {\textbf{U}}_h }{\partial k}\;dv}
 \end{equation}
 
 Using the discrete form given in {\S}II.2-a and the symmetry or anti-symmetry properties of
-the mean and difference operators, \autoref{eq:keg_zad} is demonstrated in the Appendix C.
-The main point here is that satisfying \autoref{eq:keg_zad} links the choice of the discrete forms of
+the mean and difference operators, \autoref{eq:CONS_keg_zad} is demonstrated in the Appendix C.
+The main point here is that satisfying \autoref{eq:CONS_keg_zad} links the choice of the discrete forms of
 the vertical advection and of the horizontal gradient of horizontal kinetic energy.
 Choosing one imposes the other.
@@ -127,5 +128,5 @@
 
 \[
-  % \label{eq:hpg_pe}
+  % \label{eq:CONS_hpg_pe}
   \int_D {-\frac{1}{\rho_o }\left. {\nabla p^h} \right|_z \cdot {\textbf {U}}_h \;dv} \;=\;\int_D {\nabla .\left( {\rho \,{\textbf{U}}} \right)\;g\;z\;\;dv}
 \]
@@ -133,5 +134,5 @@
 Using the discrete form given in {\S}~II.2-a and the symmetry or anti-symmetry properties of
 the mean and difference operators, (II.4.3) is demonstrated in the Appendix C.
-The main point here is that satisfying (II.4.3) strongly constraints the discrete expression of the depth of 
+The main point here is that satisfying (II.4.3) strongly constraints the discrete expression of the depth of
 $T$-points and of the term added to the pressure gradient in $s-$coordinates: the depth of a $T$-point, $z_T$,
 is defined as the sum the vertical scale factors at $w$-points starting from the surface.
@@ -145,8 +146,8 @@
 Nevertheless, the $\psi$-equation is solved numerically by an iterative solver (see {\S}~III.5),
 thus the property is only satisfied with the accuracy required on the solver.
-In addition, with the rigid-lid approximation, the change of horizontal kinetic energy due to the work of 
+In addition, with the rigid-lid approximation, the change of horizontal kinetic energy due to the work of
 surface pressure forces is exactly zero:
 \[
-  % \label{eq:spg}
+  % \label{eq:CONS_spg}
   \int_D {-\frac{1}{\rho_o }\nabla _h } \left( {p_s } \right)\cdot {\textbf{U}}_h \;dv=0
 \]
@@ -161,10 +162,10 @@
 % -------------------------------------------------------------------------------------------------------------
 \section{Conservation properties on ocean thermodynamics}
-\label{sec:Invariant_tra}
+\label{sec:CONS_Invariant_tra}
 
 In continuous formulation, the advective terms of the tracer equations conserve the tracer content and
 the quadratic form of the tracer, \ie
 \[
-  % \label{eq:tra_tra2}
+  % \label{eq:CONS_tra_tra2}
   \int_D {\nabla .\left( {T\;{\textbf{U}}} \right)\;dv} =0
   \;\text{and}
@@ -176,5 +177,5 @@
 Note that in both continuous and discrete formulations, there is generally no strict conservation of mass,
 since the equation of state is non linear with respect to $T$ and $S$.
-In practice, the mass is conserved with a very good accuracy. 
+In practice, the mass is conserved with a very good accuracy.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -182,5 +183,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Conservation properties on momentum physics}
-\label{subsec:Invariant_dyn_physics}
+\label{subsec:CONS_Invariant_dyn_physics}
 
 \textbf{* lateral momentum diffusion term}
@@ -188,5 +189,5 @@
 The continuous formulation of the horizontal diffusion of momentum satisfies the following integral constraints~:
 \[
-  % \label{eq:dynldf_dyn}
+  % \label{eq:CONS_dynldf_dyn}
   \int\limits_D {\frac{1}{e_3 }{\mathrm {\mathbf k}}\cdot \nabla \times \left[ {\nabla
         _h \left( {A^{lm}\;\chi } \right)-\nabla _h \times \left( {A^{lm}\;\zeta
@@ -195,5 +196,5 @@
 
 \[
-  % \label{eq:dynldf_div}
+  % \label{eq:CONS_dynldf_div}
   \int\limits_D {\nabla _h \cdot \left[ {\nabla _h \left( {A^{lm}\;\chi }
         \right)-\nabla _h \times \left( {A^{lm}\;\zeta \;{\mathrm {\mathbf k}}} \right)}
@@ -202,5 +203,5 @@
 
 \[
-  % \label{eq:dynldf_curl}
+  % \label{eq:CONS_dynldf_curl}
   \int_D {{\mathrm {\mathbf U}}_h \cdot \left[ {\nabla _h \left( {A^{lm}\;\chi }
         \right)-\nabla _h \times \left( {A^{lm}\;\zeta \;{\mathrm {\mathbf k}}} \right)}
@@ -209,5 +210,5 @@
 
 \[
-  % \label{eq:dynldf_curl2}
+  % \label{eq:CONS_dynldf_curl2}
   \mbox{if}\quad A^{lm}=cste\quad \quad \int_D {\zeta \;{\mathrm {\mathbf k}}\cdot
     \nabla \times \left[ {\nabla _h \left( {A^{lm}\;\chi } \right)-\nabla _h
@@ -217,5 +218,5 @@
 
 \[
-  % \label{eq:dynldf_div2}
+  % \label{eq:CONS_dynldf_div2}
   \mbox{if}\quad A^{lm}=cste\quad \quad \int_D {\chi \;\nabla _h \cdot \left[
       {\nabla _h \left( {A^{lm}\;\chi } \right)-\nabla _h \times \left(
@@ -250,5 +251,5 @@
 
 \[
-  % \label{eq:dynzdf_dyn}
+  % \label{eq:CONS_dynzdf_dyn}
   \begin{aligned}
     & \int_D {\frac{1}{e_3 }}  \frac{\partial }{\partial k}\left( \frac{A^{vm}}{e_3 }\frac{\partial {\textbf{U}}_h }{\partial k} \right) \;dv = \overrightarrow{\textbf{0}} \\
@@ -258,5 +259,5 @@
 conservation of vorticity, dissipation of enstrophy
 \[
-  % \label{eq:dynzdf_vor}
+  % \label{eq:CONS_dynzdf_vor}
   \begin{aligned}
     & \int_D {\frac{1}{e_3 }{\mathrm {\mathbf k}}\cdot \nabla \times \left( {\frac{1}{e_3
@@ -270,5 +271,5 @@
 conservation of horizontal divergence, dissipation of square of the horizontal divergence
 \[
-  % \label{eq:dynzdf_div}
+  % \label{eq:CONS_dynzdf_div}
   \begin{aligned}
     &\int_D {\nabla \cdot \left( {\frac{1}{e_3 }\frac{\partial }{\partial
@@ -289,5 +290,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Conservation properties on tracer physics}
-\label{subsec:Invariant_tra_physics}
+\label{subsec:CONS_Invariant_tra_physics}
 
 The numerical schemes used for tracer subgridscale physics are written in such a way that
@@ -296,11 +297,11 @@
 the quadratic form of these quantities (\ie\ their variance) globally tends to diminish.
 As for the advective term, there is generally no strict conservation of mass even if,
-in practice, the mass is conserved with a very good accuracy. 
-
-\textbf{* lateral physics: }conservation of tracer, dissipation of tracer 
+in practice, the mass is conserved with a very good accuracy.
+
+\textbf{* lateral physics: }conservation of tracer, dissipation of tracer
 variance, i.e.
 
 \[
-  % \label{eq:traldf_t_t2}
+  % \label{eq:CONS_traldf_t_t2}
   \begin{aligned}
     &\int_D \nabla\, \cdot\, \left( A^{lT} \,\Re \,\nabla \,T \right)\;dv = 0 \\
@@ -312,5 +313,5 @@
 
 \[
-  % \label{eq:trazdf_t_t2}
+  % \label{eq:CONS_trazdf_t_t2}
   \begin{aligned}
     & \int_D \frac{1}{e_3 } \frac{\partial }{\partial k}\left( \frac{A^{vT}}{e_3 }  \frac{\partial T}{\partial k}  \right)\;dv = 0 \\
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex	(revision 11543)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex	(revision 11544)
@@ -26,5 +26,5 @@
 To overcome problems with vanishing surface and/or bottom cells, we consider the zstar coordinate
 \[
-  % \label{eq:PE_}
+  % \label{eq:MBZ_PE_}
   z^\star = H \left( \frac{z-\eta}{H+\eta} \right)
 \]
@@ -75,5 +75,5 @@
 \section[Surface pressure gradient and sea surface heigth (\textit{dynspg.F90})]
 {Surface pressure gradient and sea surface heigth (\protect\mdl{dynspg})}
-\label{sec:DYN_hpg_spg}
+\label{sec:MBZ_dyn_hpg_spg}
 %-----------------------------------------nam_dynspg----------------------------------------------------
 
@@ -100,5 +100,5 @@
 \subsubsection[Explicit (\texttt{\textbf{key\_dynspg\_exp}})]
 {Explicit (\protect\key{dynspg\_exp})}
-\label{subsec:DYN_spg_exp}
+\label{subsec:MBZ_dyn_spg_exp}
 
 In the explicit free surface formulation, the model time step is chosen small enough to
@@ -106,5 +106,5 @@
 The sea surface height is given by:
 \begin{equation}
-  \label{eq:dynspg_ssh}
+  \label{eq:MBZ_dynspg_ssh}
   \frac{\partial \eta }{\partial t}\equiv \frac{\text{EMP}}{\rho_w }+\frac{1}{e_{1T}
     e_{2T} }\sum\limits_k {\left( {\delta_i \left[ {e_{2u} e_{3u} u}
@@ -120,5 +120,5 @@
 The surface pressure gradient, also evaluated using a leap-frog scheme, is then simply given by:
 \begin{equation}
-  \label{eq:dynspg_exp}
+  \label{eq:MBZ_dynspg_exp}
   \left\{
     \begin{aligned}
@@ -137,5 +137,5 @@
 \subsubsection[Split-explicit time-stepping (\texttt{\textbf{key\_dynspg\_ts}})]
 {Split-explicit time-stepping (\protect\key{dynspg\_ts})}
-\label{subsec:DYN_spg_ts}
+\label{subsec:MBZ_dyn_spg_ts}
 %--------------------------------------------namdom----------------------------------------------------
 
@@ -152,5 +152,5 @@
     \includegraphics[width=\textwidth]{Fig_DYN_dynspg_ts}
     \caption{
-      \protect\label{fig:DYN_dynspg_ts}
+      \protect\label{fig:MBZ_dyn_dynspg_ts}
       Schematic of the split-explicit time stepping scheme for the barotropic and baroclinic modes,
       after \citet{Griffies2004?}.
@@ -186,10 +186,10 @@
 We have
 \[
-  % \label{eq:DYN_spg_ts_eta}
+  % \label{eq:MBZ_dyn_spg_ts_eta}
   \eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1})
   = 2 \Delta t \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right]
 \]
 \begin{multline*}
-  % \label{eq:DYN_spg_ts_u}
+  % \label{eq:MBZ_dyn_spg_ts_u}
   \textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1})  \\
   = 2\Delta t \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n})
@@ -207,5 +207,5 @@
 This is also the time that sets the barotropic time steps via
 \[
-  % \label{eq:DYN_spg_ts_t}
+  % \label{eq:MBZ_dyn_spg_ts_t}
   t_n=\tau+n\Delta t
 \]
@@ -213,5 +213,5 @@
 The density scaled surface pressure is evaluated via
 \[
-  % \label{eq:DYN_spg_ts_ps}
+  % \label{eq:MBZ_dyn_spg_ts_ps}
   p_s^{(b)}(\tau,t_{n}) =
   \begin{cases}
@@ -222,5 +222,5 @@
 To get started, we assume the following initial conditions
 \[
-  % \label{eq:DYN_spg_ts_eta}
+  % \label{eq:MBZ_dyn_spg_ts_eta}
   \begin{split}
     \eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)} \\
@@ -230,5 +230,5 @@
 with
 \[
-  % \label{eq:DYN_spg_ts_etaF}
+  % \label{eq:MBZ_dyn_spg_ts_etaF}
   \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\Delta t,t_{n})
 \]
@@ -236,5 +236,5 @@
 Likewise,
 \[
-  % \label{eq:DYN_spg_ts_u}
+  % \label{eq:MBZ_dyn_spg_ts_u}
   \textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)}	\\ \\
   \textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \Delta t \ \text{RHS}_{n=0}
@@ -242,5 +242,5 @@
 with
 \[
-  % \label{eq:DYN_spg_ts_u}
+  % \label{eq:MBZ_dyn_spg_ts_u}
   \overline{\textbf{U}^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\Delta t,t_{n})
 \]
@@ -251,5 +251,5 @@
 produce the updated vertically integrated velocity at baroclinic time $\tau + \Delta \tau$
 \[
-  % \label{eq:DYN_spg_ts_u}
+  % \label{eq:MBZ_dyn_spg_ts_u}
   \textbf{U}(\tau+\Delta t) = \overline{\textbf{U}^{(b)}(\tau+\Delta t)}
   = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n})
@@ -258,5 +258,5 @@
 a baroclinic leap-frog using the following form
 \begin{equation}
-  \label{eq:DYN_spg_ts_ssh}
+  \label{eq:MBZ_dyn_spg_ts_ssh}
   \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\Delta t \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right]
 \end{equation}
@@ -267,5 +267,5 @@
 
 In general, some form of time filter is needed to maintain integrity of the surface height field due to
-the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}.
+the leap-frog splitting mode in equation \autoref{eq:MBZ_dyn_spg_ts_ssh}.
 We have tried various forms of such filtering,
 with the following method discussed in Griffies et al. (2001) chosen due to its stability and
@@ -273,5 +273,5 @@
 
 \begin{equation}
-  \label{eq:DYN_spg_ts_sshf}
+  \label{eq:MBZ_dyn_spg_ts_sshf}
   \eta^{F}(\tau-\Delta) =  \overline{\eta^{(b)}(\tau)}
 \end{equation}
@@ -279,5 +279,5 @@
 
 \[
-  % \label{eq:DYN_spg_ts_sshf2}
+  % \label{eq:MBZ_dyn_spg_ts_sshf2}
   \eta^{F}(\tau-\Delta) = \eta(\tau)
   + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\Delta t)
@@ -288,5 +288,5 @@
 This isolation allows for an easy check that tracer conservation is exact when eliminating tracer and
 surface height time filtering (see ?? for more complete discussion).
-However, in the general case with a non-zero $\alpha$, the filter \autoref{eq:DYN_spg_ts_sshf} was found to
+However, in the general case with a non-zero $\alpha$, the filter \autoref{eq:MBZ_dyn_spg_ts_sshf} was found to
 be more conservative, and so is recommended.
 
@@ -296,5 +296,5 @@
 \subsubsection[Filtered formulation (\texttt{\textbf{key\_dynspg\_flt}})]
 {Filtered formulation (\protect\key{dynspg\_flt})}
-\label{subsec:DYN_spg_flt}
+\label{subsec:MBZ_dyn_spg_flt}
 
 The filtered formulation follows the \citet{Roullet2000?} implementation.
@@ -311,5 +311,5 @@
 \subsection[Non-linear free surface formulation (\texttt{\textbf{key\_vvl}})]
 {Non-linear free surface formulation (\protect\key{vvl})}
-\label{subsec:DYN_spg_vvl}
+\label{subsec:MBZ_dyn_spg_vvl}
 
 In the non-linear free surface formulation, the variations of volume are fully taken into account.
