Changeset 11799 for NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex/NEMO/subfiles/chap_DYN.tex

Timestamp:

2019-10-25T16:27:34+02:00 (5 years ago)

Author:

mocavero

Message:

Update the branch to v4.0.1 of the trunk

Location:

NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc

Files:

• . (modified) (1 prop)
• latex (modified) (1 prop)
• latex/NEMO (modified) (1 prop)
• latex/NEMO/subfiles (modified) (1 prop)
• latex/NEMO/subfiles/chap_DYN.tex (modified) (84 diffs)

NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex/NEMO/subfiles

@@ -1 @@
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
+*.ind
 *.ilg
-*.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*

@@ -1 @@
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
+*.ind
 *.ilg
-*.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*

NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex/NEMO/subfiles/chap_DYN.tex

@@ -2 +2 @@
 
 \begin{document}
-% ================================================================
-% Chapter ——— Ocean Dynamics (DYN)
-% ================================================================
+
 \chapter{Ocean Dynamics (DYN)}
 \label{chap:DYN}
 
+\thispagestyle{plain}
+
 \chaptertoc
+
+\paragraph{Changes record} ~\\
+
+{\footnotesize
+  \begin{tabularx}{\textwidth}{l||X|X}
+    Release & Author(s) & Modifications \\
+    \hline
+    {\em   4.0} & {\em ...} & {\em ...} \\
+    {\em   3.6} & {\em ...} & {\em ...} \\
+    {\em   3.4} & {\em ...} & {\em ...} \\
+    {\em <=3.4} & {\em ...} & {\em ...}
+  \end{tabularx}
+}
+
+\clearpage
 
 Using the representation described in \autoref{chap:DOM},
@@ -52 +67 @@
 Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \texttt{trdvor?} is defined)
 can be derived from the 3D terms.
-%%%
 \gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does
 MISC correspond to "extracting tendency terms" or "vorticity balance"?}
 
-% ================================================================
-% Sea Surface Height evolution & Diagnostics variables
-% ================================================================
+%% =================================================================================================
 \section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)}
 \label{sec:DYN_divcur_wzv}
 
-%--------------------------------------------------------------------------------------------------------------
-%           Horizontal divergence and relative vorticity
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]
-{Horizontal divergence and relative vorticity (\protect\mdl{divcur})}
+%% =================================================================================================
+\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})}
 \label{subsec:DYN_divcur}
 
 The vorticity is defined at an $f$-point (\ie\ corner point) as follows:
 \begin{equation}
-  \label{eq:divcur_cur}
+  \label{eq:DYN_divcur_cur}
   \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right]
       -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right)
@@ -79 +88 @@
 It is given by:
 \[
-  % \label{eq:divcur_div}
+  % \label{eq:DYN_divcur_div}
   \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
   \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right]
@@ -99 +108 @@
 the nonlinear advection and of the vertical velocity respectively.
 
-%--------------------------------------------------------------------------------------------------------------
-%           Sea Surface Height evolution
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]
-{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}
+%% =================================================================================================
+\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}
 \label{subsec:DYN_sshwzv}
 
 The sea surface height is given by:
 \begin{equation}
-  \label{eq:dynspg_ssh}
+  \label{eq:DYN_spg_ssh}
   \begin{aligned}
     \frac{\partial \eta }{\partial t}
@@ -123 +129 @@
 \textit{emp} can be written as the evaporation minus precipitation, minus the river runoff.
 The sea-surface height is evaluated using exactly the same time stepping scheme as
-the tracer equation \autoref{eq:tra_nxt}:
+the tracer equation \autoref{eq:TRA_nxt}:
 a leapfrog scheme in combination with an Asselin time filter,
-\ie\ the velocity appearing in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity).
+\ie\ the velocity appearing in \autoref{eq:DYN_spg_ssh} is centred in time (\textit{now} velocity).
 This is of paramount importance.
 Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to
@@ -134 +140 @@
 taking into account the change of the thickness of the levels:
 \begin{equation}
-  \label{eq:wzv}
+  \label{eq:DYN_wzv}
   \left\{
     \begin{aligned}
@@ -148 +154 @@
 re-orientated downward.
 \gmcomment{not sure of this...  to be modified with the change in emp setting}
-In the case of a linear free surface, the time derivative in \autoref{eq:wzv} disappears.
+In the case of a linear free surface, the time derivative in \autoref{eq:DYN_wzv} disappears.
 The upper boundary condition applies at a fixed level $z=0$.
 The top vertical velocity is thus equal to the divergence of the barotropic transport
-(\ie\ the first term in the right-hand-side of \autoref{eq:dynspg_ssh}).
+(\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}).
 
 Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates,
 its physical meaning is not the same:
 in the second case, $w$ is the velocity normal to the $s$-surfaces.
-Note also that the $k$-axis is re-orientated downwards in the \fortran code compared to
-the indexing used in the semi-discrete equations such as \autoref{eq:wzv}
+Note also that the $k$-axis is re-orientated downwards in the \fortran\ code compared to
+the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv}
 (see \autoref{subsec:DOM_Num_Index_vertical}).
 
-
-% ================================================================
-% Coriolis and Advection terms: vector invariant form
-% ================================================================
+%% =================================================================================================
 \section{Coriolis and advection: vector invariant form}
 \label{sec:DYN_adv_cor_vect}
-%-----------------------------------------nam_dynadv----------------------------------------------------
-
-\nlst{namdyn_adv}
-%-------------------------------------------------------------------------------------------------------------
+
+\begin{listing}
+  \nlst{namdyn_adv}
+  \caption{\forcode{&namdyn_adv}}
+  \label{lst:namdyn_adv}
+\end{listing}
 
 The vector invariant form of the momentum equations is the one most often used in
 applications of the \NEMO\ ocean model.
 The flux form option (see next section) has been present since version $2$.
-Options are defined through the \nam{dyn\_adv} namelist variables Coriolis and
+Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables Coriolis and
 momentum advection terms are evaluated using a leapfrog scheme,
 \ie\ the velocity appearing in these expressions is centred in time (\textit{now} velocity).
@@ -180 +185 @@
 \autoref{chap:LBC}.
 
-% -------------------------------------------------------------------------------------------------------------
-%        Vorticity term
-% -------------------------------------------------------------------------------------------------------------
-\subsection[Vorticity term (\textit{dynvor.F90})]
-{Vorticity term (\protect\mdl{dynvor})}
+%% =================================================================================================
+\subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})}
 \label{subsec:DYN_vor}
-%------------------------------------------nam_dynvor----------------------------------------------------
-
-\nlst{namdyn_vor}
-%-------------------------------------------------------------------------------------------------------------
-
-Options are defined through the \nam{dyn\_vor} namelist variables.
-Four discretisations of the vorticity term (\texttt{ln\_dynvor\_xxx}\forcode{ = .true.}) are available:
+
+\begin{listing}
+  \nlst{namdyn_vor}
+  \caption{\forcode{&namdyn_vor}}
+  \label{lst:namdyn_vor}
+\end{listing}
+
+Options are defined through the \nam{dyn_vor}{dyn\_vor} namelist variables.
+Four discretisations of the vorticity term (\texttt{ln\_dynvor\_xxx}\forcode{=.true.}) are available:
 conserving potential enstrophy of horizontally non-divergent flow (ENS scheme);
 conserving horizontal kinetic energy (ENE scheme);
@@ -198 +202 @@
 horizontal kinetic energy for the planetary vorticity term (MIX scheme);
 or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy
-(EEN scheme) (see \autoref{subsec:C_vorEEN}).
+(EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}).
 In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of
-vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{ = .true.}).
+vorticity term with analytical equations (\np[=.true.]{ln_dynvor_con}{ln\_dynvor\_con}).
 The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module.
 
-%-------------------------------------------------------------
 %                 enstrophy conserving scheme
-%-------------------------------------------------------------
-\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens = .true.})]
-{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})}
+%% =================================================================================================
+\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens})]{Enstrophy conserving scheme (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens})}
 \label{subsec:DYN_vor_ens}
 
@@ -216 +218 @@
 It is given by:
 \begin{equation}
-  \label{eq:dynvor_ens}
+  \label{eq:DYN_vor_ens}
   \left\{
     \begin{aligned}
@@ -227 +229 @@
 \end{equation}
 
-%-------------------------------------------------------------
 %                 energy conserving scheme
-%-------------------------------------------------------------
-\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene = .true.})]
-{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})}
+%% =================================================================================================
+\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene})]{Energy conserving scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene})}
 \label{subsec:DYN_vor_ene}
 
@@ -237 +237 @@
 It is given by:
 \begin{equation}
-  \label{eq:dynvor_ene}
+  \label{eq:DYN_vor_ene}
   \left\{
     \begin{aligned}
@@ -248 +248 @@
 \end{equation}
 
-%-------------------------------------------------------------
 %                 mix energy/enstrophy conserving scheme
-%-------------------------------------------------------------
-\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix = .true.})]
-{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.})}
+%% =================================================================================================
+\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln_dynvor_mix}{ln\_dynvor\_mix})}
 \label{subsec:DYN_vor_mix}
 
 For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used.
-It consists of the ENS scheme (\autoref{eq:dynvor_ens}) for the relative vorticity term,
-and of the ENE scheme (\autoref{eq:dynvor_ene}) applied to the planetary vorticity term.
-\[
-  % \label{eq:dynvor_mix}
+It consists of the ENS scheme (\autoref{eq:DYN_vor_ens}) for the relative vorticity term,
+and of the ENE scheme (\autoref{eq:DYN_vor_ene}) applied to the planetary vorticity term.
+\[
+  % \label{eq:DYN_vor_mix}
   \left\{ {
       \begin{aligned}
@@ -274 +272 @@
 \]
 
-%-------------------------------------------------------------
 %                 energy and enstrophy conserving scheme
-%-------------------------------------------------------------
-\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een = .true.})]
-{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})}
+%% =================================================================================================
+\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een})]{Energy and enstrophy conserving scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een})}
 \label{subsec:DYN_vor_een}
 
@@ -297 +293 @@
 The idea is to get rid of the double averaging by considering triad combinations of vorticity.
 It is noteworthy that this solution is conceptually quite similar to the one proposed by
-\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}).
+\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:INVARIANTS}).
 
 The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified
@@ -303 +299 @@
 First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:
 \[
-  % \label{eq:pot_vor}
+  % \label{eq:DYN_pot_vor}
   q  = \frac{\zeta +f} {e_{3f} }
 \]
-where the relative vorticity is defined by (\autoref{eq:divcur_cur}),
+where the relative vorticity is defined by (\autoref{eq:DYN_divcur_cur}),
 the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is:
 \begin{equation}
-  \label{eq:een_e3f}
+  \label{eq:DYN_een_e3f}
   e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2}
 \end{equation}
 
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht]
-  \begin{center}
-    \includegraphics[width=\textwidth]{Fig_DYN_een_triad}
-    \caption{
-      \protect\label{fig:DYN_een_triad}
-      Triads used in the energy and enstrophy conserving scheme (een) for
-      $u$-component (upper panel) and $v$-component (lower panel).
-    }
-  \end{center}
+  \centering
+  \includegraphics[width=0.66\textwidth]{Fig_DYN_een_triad}
+  \caption[Triads used in the energy and enstrophy conserving scheme (EEN)]{
+    Triads used in the energy and enstrophy conserving scheme (EEN) for
+    $u$-component (upper panel) and $v$-component (lower panel).}
+  \label{fig:DYN_een_triad}
 \end{figure}
-% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-A key point in \autoref{eq:een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.
+
+A key point in \autoref{eq:DYN_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.
 It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks
-(\np{nn\_een\_e3f}\forcode{ = 1}), or just by $4$ (\np{nn\_een\_e3f}\forcode{ = .true.}).
+(\np[=1]{nn_een_e3f}{nn\_een\_e3f}), or just by $4$ (\np[=.true.]{nn_een_e3f}{nn\_een\_e3f}).
 The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and
 extends by continuity the value of $e_{3f}$ into the land areas.
@@ -340 +332 @@
 (\autoref{fig:DYN_een_triad}):
 \begin{equation}
-  \label{eq:Q_triads}
+  \label{eq:DYN_Q_triads}
   _i^j \mathbb{Q}^{i_p}_{j_p}
   = \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right)
@@ -348 +340 @@
 Finally, the vorticity terms are represented as:
 \begin{equation}
-  \label{eq:dynvor_een}
+  \label{eq:DYN_vor_een}
   \left\{ {
       \begin{aligned}
@@ -361 +353 @@
 This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes.
 It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow
-(\ie\ $\chi$=$0$) (see \autoref{subsec:C_vorEEN}).
+(\ie\ $\chi$=$0$) (see \autoref{subsec:INVARIANTS_vorEEN}).
 Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of
 the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}.
@@ -368 +360 @@
 leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}.
 
-%--------------------------------------------------------------------------------------------------------------
-%           Kinetic Energy Gradient term
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]
-{Kinetic energy gradient term (\protect\mdl{dynkeg})}
+%% =================================================================================================
+\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})}
 \label{subsec:DYN_keg}
 
-As demonstrated in \autoref{apdx:C},
+As demonstrated in \autoref{apdx:INVARIANTS},
 there is a single discrete formulation of the kinetic energy gradient term that,
 together with the formulation chosen for the vertical advection (see below),
 conserves the total kinetic energy:
 \[
-  % \label{eq:dynkeg}
+  % \label{eq:DYN_keg}
   \left\{
     \begin{aligned}
@@ -389 +378 @@
 \]
 
-%--------------------------------------------------------------------------------------------------------------
-%           Vertical advection term
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Vertical advection term (\textit{dynzad.F90})]
-{Vertical advection term (\protect\mdl{dynzad})}
+%% =================================================================================================
+\subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})}
 \label{subsec:DYN_zad}
 
@@ -400 +386 @@
 conserves the total kinetic energy.
 Indeed, the change of KE due to the vertical advection is exactly balanced by
-the change of KE due to the gradient of KE (see \autoref{apdx:C}).
-\[
-  % \label{eq:dynzad}
+the change of KE due to the gradient of KE (see \autoref{apdx:INVARIANTS}).
+\[
+  % \label{eq:DYN_zad}
   \left\{
     \begin{aligned}
@@ -410 +396 @@
   \right.
 \]
-When \np{ln\_dynzad\_zts}\forcode{ = .true.},
+When \np[=.true.]{ln_dynzad_zts}{ln\_dynzad\_zts},
 a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term.
 This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}.
 Note that in this case,
 a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability,
-an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}).
-
-
-% ================================================================
-% Coriolis and Advection : flux form
-% ================================================================
+an option which is only available with a TVD scheme (see \np{ln_traadv_tvd_zts}{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}).
+
+%% =================================================================================================
 \section{Coriolis and advection: flux form}
 \label{sec:DYN_adv_cor_flux}
-%------------------------------------------nam_dynadv----------------------------------------------------
-
-\nlst{namdyn_adv}
-%-------------------------------------------------------------------------------------------------------------
-
-Options are defined through the \nam{dyn\_adv} namelist variables.
+
+Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables.
 In the flux form (as in the vector invariant form),
 the Coriolis and momentum advection terms are evaluated using a leapfrog scheme,
@@ -435 +414 @@
 no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}.
 
-
-%--------------------------------------------------------------------------------------------------------------
-%           Coriolis plus curvature metric terms
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]
-{Coriolis plus curvature metric terms (\protect\mdl{dynvor})}
+%% =================================================================================================
+\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})}
 \label{subsec:DYN_cor_flux}
 
@@ -447 +422 @@
 It is given by:
 \begin{multline*}
-  % \label{eq:dyncor_metric}
+  % \label{eq:DYN_cor_metric}
   f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right)  \\
   \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right]
@@ -453 +428 @@
 \end{multline*}
 
-Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) schemes can be used to
+Any of the (\autoref{eq:DYN_vor_ens}), (\autoref{eq:DYN_vor_ene}) and (\autoref{eq:DYN_vor_een}) schemes can be used to
 compute the product of the Coriolis parameter and the vorticity.
-However, the energy-conserving scheme (\autoref{eq:dynvor_een}) has exclusively been used to date.
+However, the energy-conserving scheme (\autoref{eq:DYN_vor_een}) has exclusively been used to date.
 This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity).
 
-%--------------------------------------------------------------------------------------------------------------
-%           Flux form Advection term
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Flux form advection term (\textit{dynadv.F90})]
-{Flux form advection term (\protect\mdl{dynadv})}
+%% =================================================================================================
+\subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})}
 \label{subsec:DYN_adv_flux}
 
 The discrete expression of the advection term is given by:
 \[
-  % \label{eq:dynadv}
+  % \label{eq:DYN_adv}
   \left\{
     \begin{aligned}
@@ -487 +459 @@
 or a $3^{rd}$ order upstream biased scheme, UBS.
 The latter is described in \citet{shchepetkin.mcwilliams_OM05}.
-The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}.
+The schemes are selected using the namelist logicals \np{ln_dynadv_cen2}{ln\_dynadv\_cen2} and \np{ln_dynadv_ubs}{ln\_dynadv\_ubs}.
 In flux form, the schemes differ by the choice of a space and time interpolation to define the value of
 $u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie\ at the $T$-, $f$-,
 and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$.
 
-%-------------------------------------------------------------
 %                 2nd order centred scheme
-%-------------------------------------------------------------
-\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2 = .true.})]
-{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})}
+%% =================================================================================================
+\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln_dynadv_cen2}{ln\_dynadv\_cen2})}
 \label{subsec:DYN_adv_cen2}
 
 In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points:
 \begin{equation}
-  \label{eq:dynadv_cen2}
+  \label{eq:DYN_adv_cen2}
   \left\{
     \begin{aligned}
@@ -516 +486 @@
 so $u$ and $v$ are the \emph{now} velocities.
 
-%-------------------------------------------------------------
 %                 UBS scheme
-%-------------------------------------------------------------
-\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs = .true.})]
-{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})}
+%% =================================================================================================
+\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs})]{UBS: Upstream Biased Scheme (\protect\np{ln_dynadv_ubs}{ln\_dynadv\_ubs})}
 \label{subsec:DYN_adv_ubs}
 
@@ -527 +495 @@
 For example, the evaluation of $u_T^{ubs} $ is done as follows:
 \begin{equation}
-  \label{eq:dynadv_ubs}
+  \label{eq:DYN_adv_ubs}
   u_T^{ubs} =\overline u ^i-\;\frac{1}{6}
   \begin{cases}
@@ -542 +510 @@
 But the amplitudes of the false extrema are significantly reduced over those in the centred second order method.
 As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum
-(\ie\ \np{ln\_dynldf\_lap}\forcode{ = }\np{ln\_dynldf\_bilap}\forcode{ = .false.}),
+(\ie\ \np[=]{ln_dynldf_lap}{ln\_dynldf\_lap}\np[=.false.]{ln_dynldf_bilap}{ln\_dynldf\_bilap}),
 and it is recommended to do so.
 
 The UBS scheme is not used in all directions.
 In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie\ $u_{uw}^{ubs}$ and
-$u_{vw}^{ubs}$ in \autoref{eq:dynadv_cen2} are used.
+$u_{vw}^{ubs}$ in \autoref{eq:DYN_adv_cen2} are used.
 UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm  pursue the
 sentence:Since vertical mixing of momentum is a source term of the TKE equation...  }
 
-For stability reasons, the first term in (\autoref{eq:dynadv_ubs}),
+For stability reasons, the first term in (\autoref{eq:DYN_adv_ubs}),
 which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time),
 while the second term, which is the diffusion part of the scheme,
@@ -559 +527 @@
 Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by
 one coefficient.
-Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
+Replacing $1/6$ by $1/8$ in (\autoref{eq:DYN_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
 This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded.
 Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme.
@@ -566 +534 @@
 there is also the possibility of using a $4^{th}$ order evaluation of the advective velocity as in ROMS.
 This is an error and should be suppressed soon.
-%%%
 \gmcomment{action :  this have to be done}
-%%%
-
-% ================================================================
-%           Hydrostatic pressure gradient term
-% ================================================================
-\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]
-{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
+
+%% =================================================================================================
+\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
 \label{sec:DYN_hpg}
-%------------------------------------------nam_dynhpg---------------------------------------------------
-
-\nlst{namdyn_hpg}
-%-------------------------------------------------------------------------------------------------------------
-
-Options are defined through the \nam{dyn\_hpg} namelist variables.
+
+\begin{listing}
+  \nlst{namdyn_hpg}
+  \caption{\forcode{&namdyn_hpg}}
+  \label{lst:namdyn_hpg}
+\end{listing}
+
+Options are defined through the \nam{dyn_hpg}{dyn\_hpg} namelist variables.
 The key distinction between the different algorithms used for
 the hydrostatic pressure gradient is the vertical coordinate used,
@@ -593 +558 @@
 At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied.
 
-%--------------------------------------------------------------------------------------------------------------
-%           z-coordinate with full step
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco = .true.})]
-{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})}
+%% =================================================================================================
+\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco})]{Full step $Z$-coordinate (\protect\np{ln_dynhpg_zco}{ln\_dynhpg\_zco})}
 \label{subsec:DYN_hpg_zco}
 
@@ -607 +569 @@
 for $k=km$ (surface layer, $jk=1$ in the code)
 \begin{equation}
-  \label{eq:dynhpg_zco_surf}
+  \label{eq:DYN_hpg_zco_surf}
   \left\{
     \begin{aligned}
@@ -620 +582 @@
 for $1<k<km$ (interior layer)
 \begin{equation}
-  \label{eq:dynhpg_zco}
+  \label{eq:DYN_hpg_zco}
   \left\{
     \begin{aligned}
@@ -633 +595 @@
 \end{equation}
 
-Note that the $1/2$ factor in (\autoref{eq:dynhpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as
+Note that the $1/2$ factor in (\autoref{eq:DYN_hpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as
 the vertical derivative of the scale factor at the surface level ($z=0$).
 Note also that in case of variable volume level (\texttt{vvl?} defined),
-the surface pressure gradient is included in \autoref{eq:dynhpg_zco_surf} and
-\autoref{eq:dynhpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$.
-
-%--------------------------------------------------------------------------------------------------------------
-%           z-coordinate with partial step
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps = .true.})]
-{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})}
+the surface pressure gradient is included in \autoref{eq:DYN_hpg_zco_surf} and
+\autoref{eq:DYN_hpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$.
+
+%% =================================================================================================
+\subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps})]{Partial step $Z$-coordinate (\protect\np{ln_dynhpg_zps}{ln\_dynhpg\_zps})}
 \label{subsec:DYN_hpg_zps}
 
@@ -661 +620 @@
 module \mdl{zpsdhe} located in the TRA directory and described in \autoref{sec:TRA_zpshde}.
 
-%--------------------------------------------------------------------------------------------------------------
-%           s- and s-z-coordinates
-%--------------------------------------------------------------------------------------------------------------
+%% =================================================================================================
 \subsection{$S$- and $Z$-$S$-coordinates}
 \label{subsec:DYN_hpg_sco}
@@ -672 +629 @@
 density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation.
 
-$\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.})
-\begin{equation}
-  \label{eq:dynhpg_sco}
+$\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np[=.true.]{ln_dynhpg_sco}{ln\_dynhpg\_sco})
+\begin{equation}
+  \label{eq:DYN_hpg_sco}
   \left\{
     \begin{aligned}
@@ -686 +643 @@
 
 Where the first term is the pressure gradient along coordinates,
-computed as in \autoref{eq:dynhpg_zco_surf} - \autoref{eq:dynhpg_zco},
+computed as in \autoref{eq:DYN_hpg_zco_surf} - \autoref{eq:DYN_hpg_zco},
 and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point
 ($e_{3w}$).
 
-$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}\forcode{ = .true.}).
-This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}\forcode{ = .true.}).
-
-$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}\forcode{ = .true.})
+$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np[=.true.]{ln_dynhpg_isf}{ln\_dynhpg\_isf}).
+This scheme need the activation of ice shelf cavities (\np[=.true.]{ln_isfcav}{ln\_isfcav}).
+
+$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np[=.true.]{ln_dynhpg_prj}{ln\_dynhpg\_prj})
 
 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{shchepetkin.mcwilliams_OM05}
-(\np{ln\_dynhpg\_djc}\forcode{ = .true.}) (currently disabled; under development)
-
-Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is activated
+(\np[=.true.]{ln_dynhpg_djc}{ln\_dynhpg\_djc}) (currently disabled; under development)
+
+Note that expression \autoref{eq:DYN_hpg_sco} is commonly used when the variable volume formulation is activated
 (\texttt{vvl?}) because in that case, even with a flat bottom,
 the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}.
-The pressure jacobian scheme (\np{ln\_dynhpg\_prj}\forcode{ = .true.}) is available as
-an improved option to \np{ln\_dynhpg\_sco}\forcode{ = .true.} when \texttt{vvl?} is active.
+The pressure jacobian scheme (\np[=.true.]{ln_dynhpg_prj}{ln\_dynhpg\_prj}) is available as
+an improved option to \np[=.true.]{ln_dynhpg_sco}{ln\_dynhpg\_sco} when \texttt{vvl?} is active.
 The pressure Jacobian scheme uses a constrained cubic spline to
 reconstruct the density profile across the water column.
@@ -710 +667 @@
 This method can provide a more accurate calculation of the horizontal pressure gradient than the standard scheme.
 
+%% =================================================================================================
 \subsection{Ice shelf cavity}
 \label{subsec:DYN_hpg_isf}
+
 Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and
-the pressure gradient due to the ocean load (\np{ln\_dynhpg\_isf}\forcode{ = .true.}).\\
+the pressure gradient due to the ocean load (\np[=.true.]{ln_dynhpg_isf}{ln\_dynhpg\_isf}).\\
 
 The main hypothesis to compute the ice shelf load is that the ice shelf is in an isostatic equilibrium.
@@ -722 +681 @@
 A detailed description of this method is described in \citet{losch_JGR08}.\\
 
-The pressure gradient due to ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in
+The pressure gradient due to ocean load is computed using the expression \autoref{eq:DYN_hpg_sco} described in
 \autoref{subsec:DYN_hpg_sco}.
 
-%--------------------------------------------------------------------------------------------------------------
-%           Time-scheme
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Time-scheme (\forcode{ln_dynhpg_imp = .{true,false}.})]
-{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .\{true,false\}}.)}
+%% =================================================================================================
+\subsection[Time-scheme (\forcode{ln_dynhpg_imp})]{Time-scheme (\protect\np{ln_dynhpg_imp}{ln\_dynhpg\_imp})}
 \label{subsec:DYN_hpg_imp}
 
@@ -745 +701 @@
 rather than at the central time level $t$ only, as in the standard leapfrog scheme.
 
-$\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}):
-
-\begin{equation}
-  \label{eq:dynhpg_lf}
+$\bullet$ leapfrog scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}):
+
+\begin{equation}
+  \label{eq:DYN_hpg_lf}
   \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
   -\frac{1}{\rho_o \,e_{1u} }\delta_{i+1/2} \left[ {p_h^t } \right]
 \end{equation}
 
-$\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}):
-\begin{equation}
-  \label{eq:dynhpg_imp}
+$\bullet$ semi-implicit scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}):
+\begin{equation}
+  \label{eq:DYN_hpg_imp}
   \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
   -\frac{1}{4\,\rho_o \,e_{1u} } \delta_{i+1/2} \left[ p_h^{t+\rdt} +2\,p_h^t +p_h^{t-\rdt}  \right]
 \end{equation}
 
-The semi-implicit time scheme \autoref{eq:dynhpg_imp} is made possible without
+The semi-implicit time scheme \autoref{eq:DYN_hpg_imp} is made possible without
 significant additional computation since the density can be updated to time level $t+\rdt$ before
 computing the horizontal hydrostatic pressure gradient.
 It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using
-\autoref{eq:dynhpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:dynhpg_lf}.
-Note that \autoref{eq:dynhpg_imp} is equivalent to applying a time filter to the pressure gradient to
+\autoref{eq:DYN_hpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:DYN_hpg_lf}.
+Note that \autoref{eq:DYN_hpg_imp} is equivalent to applying a time filter to the pressure gradient to
 eliminate high frequency IGWs.
-Obviously, when using \autoref{eq:dynhpg_imp},
+Obviously, when using \autoref{eq:DYN_hpg_imp},
 the doubling of the time-step is achievable only if no other factors control the time-step,
 such as the stability limits associated with advection or diffusion.
 
-In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}\forcode{ = .true.}.
+In practice, the semi-implicit scheme is used when \np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}.
 In this case, we choose to apply the time filter to temperature and salinity used in the equation of state,
 instead of applying it to the hydrostatic pressure or to the density,
@@ -777 +733 @@
 The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows:
 \[
-  % \label{eq:rho_flt}
+  % \label{eq:DYN_rho_flt}
   \rho^t = \rho( \widetilde{T},\widetilde {S},z_t)
   \quad    \text{with}  \quad
@@ -785 +741 @@
 Note that in the semi-implicit case, it is necessary to save the filtered density,
 an extra three-dimensional field, in the restart file to restart the model with exact reproducibility.
-This option is controlled by  \np{nn\_dynhpg\_rst}, a namelist parameter.
-
-% ================================================================
-% Surface Pressure Gradient
-% ================================================================
-\section[Surface pressure gradient (\textit{dynspg.F90})]
-{Surface pressure gradient (\protect\mdl{dynspg})}
+This option is controlled by  \np{nn_dynhpg_rst}{nn\_dynhpg\_rst}, a namelist parameter.
+
+%% =================================================================================================
+\section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})}
 \label{sec:DYN_spg}
-%-----------------------------------------nam_dynspg----------------------------------------------------
-
-\nlst{namdyn_spg}
-%------------------------------------------------------------------------------------------------------------
-
-Options are defined through the \nam{dyn\_spg} namelist variables.
-The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:PE_hor_pg}).
+
+\begin{listing}
+  \nlst{namdyn_spg}
+  \caption{\forcode{&namdyn_spg}}
+  \label{lst:namdyn_spg}
+\end{listing}
+
+Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables.
+The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:MB_hor_pg}).
 The main distinction is between the fixed volume case (linear free surface) and
 the variable volume case (nonlinear free surface, \texttt{vvl?} is defined).
-In the linear free surface case (\autoref{subsec:PE_free_surface})
+In the linear free surface case (\autoref{subsec:MB_free_surface})
 the vertical scale factors $e_{3}$ are fixed in time,
-while they are time-dependent in the nonlinear case (\autoref{subsec:PE_free_surface}).
+while they are time-dependent in the nonlinear case (\autoref{subsec:MB_free_surface}).
 With both linear and nonlinear free surface, external gravity waves are allowed in the equations,
 which imposes a very small time step when an explicit time stepping is used.
 Two methods are proposed to allow a longer time step for the three-dimensional equations:
-the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt?}),
+the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:MB_flt?}),
 and the split-explicit free surface described below.
 The extra term introduced in the filtered method is calculated implicitly,
 so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
 
-
 The form of the surface pressure gradient term depends on how the user wants to
-handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:PE_hor_pg}).
+handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:MB_hor_pg}).
 Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx):
 an explicit formulation which requires a small time step;
@@ -825 +779 @@
 As a consequence the update of the $next$ velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
 
-
-%--------------------------------------------------------------------------------------------------------------
-% Explicit free surface formulation
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Explicit free surface (\texttt{ln\_dynspg\_exp}\forcode{ = .true.})]
-{Explicit free surface (\protect\np{ln\_dynspg\_exp}\forcode{ = .true.})}
+%% =================================================================================================
+\subsection[Explicit free surface (\forcode{ln_dynspg_exp})]{Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})}
 \label{subsec:DYN_spg_exp}
 
-In the explicit free surface formulation (\np{ln\_dynspg\_exp} set to true),
+In the explicit free surface formulation (\np{ln_dynspg_exp}{ln\_dynspg\_exp} set to true),
 the model time step is chosen to be small enough to resolve the external gravity waves
 (typically a few tens of seconds).
@@ -839 +789 @@
 is thus simply given by :
 \begin{equation}
-  \label{eq:dynspg_exp}
+  \label{eq:DYN_spg_exp}
   \left\{
     \begin{aligned}
@@ -853 +803 @@
 Thus, nothing is done in the \mdl{dynspg\_exp} module.
 
-%--------------------------------------------------------------------------------------------------------------
-% Split-explict free surface formulation
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Split-explicit free surface (\texttt{ln\_dynspg\_ts}\forcode{ = .true.})]
-{Split-explicit free surface (\protect\np{ln\_dynspg\_ts}\forcode{ = .true.})}
+%% =================================================================================================
+\subsection[Split-explicit free surface (\forcode{ln_dynspg_ts})]{Split-explicit free surface (\protect\np{ln_dynspg_ts}{ln\_dynspg\_ts})}
 \label{subsec:DYN_spg_ts}
-%------------------------------------------namsplit-----------------------------------------------------------
-%
 %\nlst{namsplit}
-%-------------------------------------------------------------------------------------------------------------
-
-The split-explicit free surface formulation used in \NEMO\ (\np{ln\_dynspg\_ts} set to true),
+
+The split-explicit free surface formulation used in \NEMO\ (\np{ln_dynspg_ts}{ln\_dynspg\_ts} set to true),
 also called the time-splitting formulation, follows the one proposed by \citet{shchepetkin.mcwilliams_OM05}.
 The general idea is to solve the free surface equation and the associated barotropic velocity equations with
 a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables
-(\autoref{fig:DYN_dynspg_ts}).
+(\autoref{fig:DYN_spg_ts}).
 The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through
-the \np{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$.
-This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}\forcode{ = .true.}) considering that
+the \np{nn_baro}{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$.
+This parameter can be optionally defined automatically (\np[=.true.]{ln_bt_nn_auto}{ln\_bt\_nn\_auto}) considering that
 the stability of the barotropic system is essentially controled by external waves propagation.
 Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry.
-Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn\_bt\_cmax}.
-
-%%%
+Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn_bt_cmax}{rn\_bt\_cmax}.
+
 The barotropic mode solves the following equations:
 % \begin{subequations}
-%  \label{eq:BT}
-\begin{equation}
-  \label{eq:BT_dyn}
+%  \label{eq:DYN_BT}
+\begin{equation}
+  \label{eq:DYN_BT_dyn}
   \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}=
   -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h}
@@ -887 +830 @@
 \end{equation}
 \[
-  % \label{eq:BT_ssh}
+  % \label{eq:DYN_BT_ssh}
   \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E
 \]
@@ -893 +836 @@
 where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing coupling term between modes,
 surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency.
-The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress
-(see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration.
+The third term on the right hand side of \autoref{eq:DYN_BT_dyn} represents the bottom stress
+(see section \autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration.
 Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm
 detailed in \citet{shchepetkin.mcwilliams_OM05}.
@@ -901 +844 @@
 (see their figure 12, lower left).
 
-%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
 \begin{figure}[!t]
-  \begin{center}
-    \includegraphics[width=\textwidth]{Fig_DYN_dynspg_ts}
-    \caption{
-      \protect\label{fig:DYN_dynspg_ts}
-      Schematic of the split-explicit time stepping scheme for the external and internal modes.
-      Time increases to the right. In this particular exemple,
-      a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_flt=1$) and $nn\_baro=5$.
-      Internal mode time steps (which are also the model time steps) are denoted by $t-\rdt$, $t$ and $t+\rdt$.
-      Variables with $k$ superscript refer to instantaneous barotropic variables,
-      $< >$ and $<< >>$ operator refer to time filtered variables using respectively primary (red vertical bars) and
-      secondary weights (blue vertical bars).
-      The former are used to obtain time filtered quantities at $t+\rdt$ while
-      the latter are used to obtain time averaged transports to advect tracers.
-      a) Forward time integration: \protect\np{ln\_bt\_fw}\forcode{ = .true.},
-      \protect\np{ln\_bt\_av}\forcode{ = .true.}.
-      b) Centred time integration: \protect\np{ln\_bt\_fw}\forcode{ = .false.},
-      \protect\np{ln\_bt\_av}\forcode{ = .true.}.
-      c) Forward time integration with no time filtering (POM-like scheme):
-      \protect\np{ln\_bt\_fw}\forcode{ = .true.}, \protect\np{ln\_bt\_av}\forcode{ = .false.}.
-    }
-  \end{center}
+  \centering
+  \includegraphics[width=0.66\textwidth]{Fig_DYN_dynspg_ts}
+  \caption[Split-explicit time stepping scheme for the external and internal modes]{
+    Schematic of the split-explicit time stepping scheme for the external and internal modes.
+    Time increases to the right.
+    In this particular exemple,
+    a boxcar averaging window over \np{nn_baro}{nn\_baro} barotropic time steps is used
+    (\np[=1]{nn_bt_flt}{nn\_bt\_flt}) and \np[=5]{nn_baro}{nn\_baro}.
+    Internal mode time steps (which are also the model time steps) are denoted by
+    $t-\rdt$, $t$ and $t+\rdt$.
+    Variables with $k$ superscript refer to instantaneous barotropic variables,
+    $< >$ and $<< >>$ operator refer to time filtered variables using respectively primary
+    (red vertical bars) and secondary weights (blue vertical bars).
+    The former are used to obtain time filtered quantities at $t+\rdt$ while
+    the latter are used to obtain time averaged transports to advect tracers.
+    a) Forward time integration:
+    \protect\np[=.true.]{ln_bt_fw}{ln\_bt\_fw},  \protect\np[=.true.]{ln_bt_av}{ln\_bt\_av}.
+    b) Centred time integration:
+    \protect\np[=.false.]{ln_bt_fw}{ln\_bt\_fw}, \protect\np[=.true.]{ln_bt_av}{ln\_bt\_av}.
+    c) Forward time integration with no time filtering (POM-like scheme):
+    \protect\np[=.true.]{ln_bt_fw}{ln\_bt\_fw},  \protect\np[=.false.]{ln_bt_av}{ln\_bt\_av}.}
+  \label{fig:DYN_spg_ts}
 \end{figure}
-%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
-
-In the default case (\np{ln\_bt\_fw}\forcode{ = .true.}),
+
+In the default case (\np[=.true.]{ln_bt_fw}{ln\_bt\_fw}),
 the external mode is integrated between \textit{now} and \textit{after} baroclinic time-steps
-(\autoref{fig:DYN_dynspg_ts}a).
+(\autoref{fig:DYN_spg_ts}a).
 To avoid aliasing of fast barotropic motions into three dimensional equations,
-time filtering is eventually applied on barotropic quantities (\np{ln\_bt\_av}\forcode{ = .true.}).
+time filtering is eventually applied on barotropic quantities (\np[=.true.]{ln_bt_av}{ln\_bt\_av}).
 In that case, the integration is extended slightly beyond \textit{after} time step to
 provide time filtered quantities.
@@ -938 +880 @@
 asselin filtering is not applied to barotropic quantities.\\
 Alternatively, one can choose to integrate barotropic equations starting from \textit{before} time step
-(\np{ln\_bt\_fw}\forcode{ = .false.}).
-Although more computationaly expensive ( \np{nn\_baro} additional iterations are indeed necessary),
+(\np[=.false.]{ln_bt_fw}{ln\_bt\_fw}).
+Although more computationaly expensive ( \np{nn_baro}{nn\_baro} additional iterations are indeed necessary),
 the baroclinic to barotropic forcing term given at \textit{now} time step become centred in
 the middle of the integration window.
@@ -946 +888 @@
 %references to Patrick Marsaleix' work here. Also work done by SHOM group.
 
-%%%
 
 As far as tracer conservation is concerned,
@@ -960 +901 @@
 obtain exact conservation.
 
-%%%
 
 One can eventually choose to feedback instantaneous values by not using any time filter
-(\np{ln\_bt\_av}\forcode{ = .false.}).
+(\np[=.false.]{ln_bt_av}{ln\_bt\_av}).
 In that case, external mode equations are continuous in time,
 \ie\ they are not re-initialized when starting a new sub-stepping sequence.
@@ -975 +915 @@
 it is still significant as shown by \citet{levier.treguier.ea_rpt07} in the case of an analytical barotropic Kelvin wave.
 
-%>>>>>===============
 \gmcomment{               %%% copy from griffies Book
 
@@ -1095 +1034 @@
 }            %%end gm comment (copy of griffies book)
 
-%>>>>>===============
-
-
-%--------------------------------------------------------------------------------------------------------------
-% Filtered free surface formulation
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Filtered free surface (\texttt{dynspg\_flt?})]
-{Filtered free surface (\protect\texttt{dynspg\_flt?})}
+%% =================================================================================================
+\subsection{Filtered free surface (\forcode{dynspg_flt?})}
 \label{subsec:DYN_spg_fltp}
 
 The filtered formulation follows the \citet{roullet.madec_JGR00} implementation.
-The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly.
+The extra term introduced in the equations (see \autoref{subsec:MB_free_surface}) is solved implicitly.
 The elliptic solvers available in the code are documented in \autoref{chap:MISC}.
 
@@ -1112 +1045 @@
 \gmcomment{               %%% copy from chap-model basics
   \[
-    % \label{eq:spg_flt}
+    % \label{eq:DYN_spg_flt}
     \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t}= {\mathrm {\mathbf M}}
     - g \nabla \left( \tilde{\rho} \ \eta \right)
@@ -1120 +1053 @@
   $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density,
   and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient,
-  non-linear and viscous terms in \autoref{eq:PE_dyn}.
+  non-linear and viscous terms in \autoref{eq:MB_dyn}.
 }   %end gmcomment
 
@@ -1127 +1060 @@
 It is computed once and for all and applies to all ocean time steps.
 
-% ================================================================
-% Lateral diffusion term
-% ================================================================
-\section[Lateral diffusion term and operators (\textit{dynldf.F90})]
-{Lateral diffusion term and operators (\protect\mdl{dynldf})}
+%% =================================================================================================
+\section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})}
 \label{sec:DYN_ldf}
-%------------------------------------------nam_dynldf----------------------------------------------------
-
-\nlst{namdyn_ldf}
-%-------------------------------------------------------------------------------------------------------------
-
-Options are defined through the \nam{dyn\_ldf} namelist variables.
+
+\begin{listing}
+  \nlst{namdyn_ldf}
+  \caption{\forcode{&namdyn_ldf}}
+  \label{lst:namdyn_ldf}
+\end{listing}
+
+Options are defined through the \nam{dyn_ldf}{dyn\_ldf} namelist variables.
 The options available for lateral diffusion are to use either laplacian (rotated or not) or biharmonic operators.
 The coefficients may be constant or spatially variable;
@@ -1145 +1077 @@
 \ie\ the velocity appearing in its expression is the \textit{before} velocity in time,
 except for the pure vertical component that appears when a tensor of rotation is used.
-This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}).
+This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}).
 
 At the lateral boundaries either free slip,
@@ -1164 +1096 @@
 }
 
-% ================================================================
-\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap = .true.})]
-{Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})}
+%% =================================================================================================
+\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap})]{Iso-level laplacian operator (\protect\np{ln_dynldf_lap}{ln\_dynldf\_lap})}
 \label{subsec:DYN_ldf_lap}
 
 For lateral iso-level diffusion, the discrete operator is:
 \begin{equation}
-  \label{eq:dynldf_lap}
+  \label{eq:DYN_ldf_lap}
   \left\{
     \begin{aligned}
@@ -1184 +1115 @@
 \end{equation}
 
-As explained in \autoref{subsec:PE_ldf},
+As explained in \autoref{subsec:MB_ldf},
 this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and
 ensures a complete separation between the vorticity and divergence parts of the momentum diffusion.
 
-%--------------------------------------------------------------------------------------------------------------
-%           Rotated laplacian operator
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Rotated laplacian (\forcode{ln_dynldf_iso = .true.})]
-{Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})}
+%% =================================================================================================
+\subsection[Rotated laplacian (\forcode{ln_dynldf_iso})]{Rotated laplacian operator (\protect\np{ln_dynldf_iso}{ln\_dynldf\_iso})}
 \label{subsec:DYN_ldf_iso}
 
 A rotation of the lateral momentum diffusion operator is needed in several cases:
-for iso-neutral diffusion in the $z$-coordinate (\np{ln\_dynldf\_iso}\forcode{ = .true.}) and
-for either iso-neutral (\np{ln\_dynldf\_iso}\forcode{ = .true.}) or
-geopotential (\np{ln\_dynldf\_hor}\forcode{ = .true.}) diffusion in the $s$-coordinate.
+for iso-neutral diffusion in the $z$-coordinate (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) and
+for either iso-neutral (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) or
+geopotential (\np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}) diffusion in the $s$-coordinate.
 In the partial step case, coordinates are horizontal except at the deepest level and
-no rotation is performed when \np{ln\_dynldf\_hor}\forcode{ = .true.}.
+no rotation is performed when \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}.
 The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on
 each momentum component.
@@ -1206 +1134 @@
 The resulting discrete representation is:
 \begin{equation}
-  \label{eq:dyn_ldf_iso}
+  \label{eq:DYN_ldf_iso}
   \begin{split}
     D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\
@@ -1247 +1175 @@
 The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}).
 
-%--------------------------------------------------------------------------------------------------------------
-%           Iso-level bilaplacian operator
-%--------------------------------------------------------------------------------------------------------------
-\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap = .true.})]
-{Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})}
+%% =================================================================================================
+\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap})]{Iso-level bilaplacian operator (\protect\np{ln_dynldf_bilap}{ln\_dynldf\_bilap})}
 \label{subsec:DYN_ldf_bilap}
 
-The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:dynldf_lap} twice.
+The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice.
 It requires an additional assumption on boundary conditions:
 the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen,
 while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}).
-%%%
 \gmcomment{add a remark on the the change in the position of the coefficient}
-%%%
-
-% ================================================================
-%           Vertical diffusion term
-% ================================================================
-\section[Vertical diffusion term (\textit{dynzdf.F90})]
-{Vertical diffusion term (\protect\mdl{dynzdf})}
+
+%% =================================================================================================
+\section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})}
 \label{sec:DYN_zdf}
-%----------------------------------------------namzdf------------------------------------------------------
-
-\nlst{namzdf}
-%-------------------------------------------------------------------------------------------------------------
-
-Options are defined through the \nam{zdf} namelist variables.
+
+Options are defined through the \nam{zdf}{zdf} namelist variables.
 The large vertical diffusion coefficient found in the surface mixed layer together with high vertical resolution implies that in the case of explicit time stepping there would be too restrictive a constraint on the time step.
 Two time stepping schemes can be used for the vertical diffusion term:
 $(a)$ a forward time differencing scheme
-(\np{ln\_zdfexp}\forcode{ = .true.}) using a time splitting technique (\np{nn\_zdfexp} $>$ 1) or
-$(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp}\forcode{ = .false.})
-(see \autoref{chap:STP}).
-Note that namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.
+(\np[=.true.]{ln_zdfexp}{ln\_zdfexp}) using a time splitting technique (\np{nn_zdfexp}{nn\_zdfexp} $>$ 1) or
+$(b)$ a backward (or implicit) time differencing scheme (\np[=.false.]{ln_zdfexp}{ln\_zdfexp})
+(see \autoref{chap:TD}).
+Note that namelist variables \np{ln_zdfexp}{ln\_zdfexp} and \np{nn_zdfexp}{nn\_zdfexp} apply to both tracers and dynamics.
 
 The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is.
-The vertical diffusion operators given by \autoref{eq:PE_zdf} take the following semi-discrete space form:
-\[
-  % \label{eq:dynzdf}
+The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form:
+\[
+  % \label{eq:DYN_zdf}
   \left\{
     \begin{aligned}
@@ -1303 +1219 @@
 the vertical turbulent momentum fluxes,
 \begin{equation}
-  \label{eq:dynzdf_sbc}
+  \label{eq:DYN_zdf_sbc}
   \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}
   = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v }
@@ -1316 +1232 @@
 
 The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation
-(see \autoref{sec:ZDF_bfr})
-
-% ================================================================
-% External Forcing
-% ================================================================
+(see \autoref{sec:ZDF_drg})
+
+%% =================================================================================================
 \section{External forcings}
 \label{sec:DYN_forcing}
@@ -1328 +1242 @@
 three other forcings may enter the dynamical equations by affecting the surface pressure gradient.
 
-(1) When \np{ln\_apr\_dyn}\forcode{ = .true.} (see \autoref{sec:SBC_apr}),
+(1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}),
 the atmospheric pressure is taken into account when computing the surface pressure gradient.
 
-(2) When \np{ln\_tide\_pot}\forcode{ = .true.} and \np{ln\_tide}\forcode{ = .true.} (see \autoref{sec:SBC_tide}),
+(2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see \autoref{sec:SBC_tide}),
 the tidal potential is taken into account when computing the surface pressure gradient.
 
-(3) When \np{nn\_ice\_embd}\forcode{ = 2} and LIM or CICE is used
+(3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and LIM or CICE is used
 (\ie\ when the sea-ice is embedded in the ocean),
 the snow-ice mass is taken into account when computing the surface pressure gradient.
 
-
 \gmcomment{ missing : the lateral boundary condition !!!   another external forcing
  }
 
-% ================================================================
-% Wetting and drying
-% ================================================================
+%% =================================================================================================
 \section{Wetting and drying }
 \label{sec:DYN_wetdry}
+
 There are two main options for wetting and drying code (wd):
 (a) an iterative limiter (il) and (b) a directional limiter (dl).
@@ -1356 +1268 @@
 by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$.
 
-\nlst{namwad}
+\begin{listing}
+  \nlst{namwad}
+  \caption{\forcode{&namwad}}
+  \label{lst:namwad}
+\end{listing}
 
 The following terminology is used. The depth of the topography (positive downwards)
@@ -1391 +1307 @@
 The final sub-section covers some additional considerations that are relevant to both schemes.
 
-
-%-----------------------------------------------------------------------------------------
 %   Iterative limiters
-%-----------------------------------------------------------------------------------------
-\subsection[Directional limiter (\textit{wet\_dry.F90})]
-{Directional limiter (\mdl{wet\_dry})}
+%% =================================================================================================
+\subsection[Directional limiter (\textit{wet\_dry.F90})]{Directional limiter (\mdl{wet\_dry})}
 \label{subsec:DYN_wd_directional_limiter}
+
 The principal idea of the directional limiter is that
-water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than \np{rn\_wdmin1}).
+water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than \np{rn_wdmin1}{rn\_wdmin1}).
 
 All the changes associated with this option are made to the barotropic solver for the non-linear
@@ -1409 +1323 @@
 
 The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj).
-If the user sets \np{ln\_wd\_dl\_ramp}\forcode{ = .false.} then zuwdmask is 1 when the
-flux is from a cell with water depth greater than \np{rn\_wdmin1} and 0 otherwise. If the user sets
-\np{ln\_wd\_dl\_ramp}\forcode{ = .true.} the flux across the face is ramped down as the water depth decreases
-from 2 * \np{rn\_wdmin1} to \np{rn\_wdmin1}. The use of this ramp reduced grid-scale noise in idealised test cases.
+If the user sets \np[=.false.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} then zuwdmask is 1 when the
+flux is from a cell with water depth greater than \np{rn_wdmin1}{rn\_wdmin1} and 0 otherwise. If the user sets
+\np[=.true.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} the flux across the face is ramped down as the water depth decreases
+from 2 * \np{rn_wdmin1}{rn\_wdmin1} to \np{rn_wdmin1}{rn\_wdmin1}. The use of this ramp reduced grid-scale noise in idealised test cases.
 
 At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen
@@ -1422 +1336 @@
 treatment in the calculation of the flux of mass across the cell face.
 
-
 \cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic
 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than
@@ -1428 +1341 @@
 fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because
 the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts
-to equal their mean value during the barotropic steps. If the user sets \np{ln\_wd\_dl\_bc}\forcode{ = .true.}, the
+to equal their mean value during the barotropic steps. If the user sets \np[=.true.]{ln_wd_dl_bc}{ln\_wd\_dl\_bc}, the
 baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask.
 
-%-----------------------------------------------------------------------------------------
 %   Iterative limiters
-%-----------------------------------------------------------------------------------------
-
-\subsection[Iterative limiter (\textit{wet\_dry.F90})]
-{Iterative limiter (\mdl{wet\_dry})}
+
+%% =================================================================================================
+\subsection[Iterative limiter (\textit{wet\_dry.F90})]{Iterative limiter (\mdl{wet\_dry})}
 \label{subsec:DYN_wd_iterative_limiter}
 
-\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]
-{Iterative flux limiter (\mdl{wet\_dry})}
-\label{subsubsec:DYN_wd_il_spg_limiter}
+%% =================================================================================================
+\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]{Iterative flux limiter (\mdl{wet\_dry})}
+\label{subsec:DYN_wd_il_spg_limiter}
 
 The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry''
@@ -1449 +1360 @@
 
 The continuity equation for the total water depth in a column
-\begin{equation} \label{dyn_wd_continuity}
- \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .
+\begin{equation}
+  \label{eq:DYN_wd_continuity}
+  \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .
 \end{equation}
 can be written in discrete form  as
 
-\begin{align} \label{dyn_wd_continuity_2}
-\frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) )
-&= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\
-&= \mathrm{zzflx}_{i,j} .
+\begin{align}
+  \label{eq:DYN_wd_continuity_2}
+  \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) )
+  &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\
+  &= \mathrm{zzflx}_{i,j} .
 \end{align}
 
@@ -1470 +1383 @@
 (zzflxp) and fluxes that are into the cell (zzflxn).  Clearly
 
-\begin{equation} \label{dyn_wd_zzflx_p_n_1}
-\mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .
+\begin{equation}
+  \label{eq:DYN_wd_zzflx_p_n_1}
+  \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .
 \end{equation}
 
@@ -1482 +1396 @@
 $\mathrm{zcoef}_{i,j}^{(m)}$ such that:
 
-\begin{equation} \label{dyn_wd_continuity_coef}
-\begin{split}
-\mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\
-\mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}
-\end{split}
+\begin{equation}
+  \label{eq:DYN_wd_continuity_coef}
+  \begin{split}
+    \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\
+    \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}
+  \end{split}
 \end{equation}
 
@@ -1494 +1409 @@
 The iteration is initialised by setting
 
-\begin{equation} \label{dyn_wd_zzflx_initial}
-\mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} .
+\begin{equation}
+  \label{eq:DYN_wd_zzflx_initial}
+  \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} .
 \end{equation}
 
 The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the
 cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell
-times the timestep divided by the cell area. Using (\ref{dyn_wd_continuity_2}) this
+times the timestep divided by the cell area. Using (\autoref{eq:DYN_wd_continuity_2}) this
 condition is
 
-\begin{equation} \label{dyn_wd_continuity_if}
-h_{i,j}(t_e)  - \mathrm{rn\_wdmin1} <  \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) .
-\end{equation}
-
-Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum
+\begin{equation}
+  \label{eq:DYN_wd_continuity_if}
+  h_{i,j}(t_e)  - \mathrm{rn\_wdmin1} <  \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) .
+\end{equation}
+
+Rearranging (\autoref{eq:DYN_wd_continuity_if}) we can obtain an expression for the maximum
 outward flux that can be allowed and still maintain the minimum wet depth:
 
-\begin{equation} \label{dyn_wd_max_flux}
-\begin{split}
-\mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
-\phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big]
-\end{split}
+\begin{equation}
+  \label{eq:DYN_wd_max_flux}
+  \begin{split}
+    \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
+    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big]
+  \end{split}
 \end{equation}
 
 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is
-this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an
+this necessary/desirable?]}. Substituting from (\autoref{eq:DYN_wd_continuity_coef}) gives an
 expression for the coefficient needed to multiply the outward flux at this cell in order
 to avoid drying.
 
-\begin{equation} \label{dyn_wd_continuity_nxtcoef}
-\begin{split}
-\mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
-\phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} }
-\end{split}
+\begin{equation}
+  \label{eq:DYN_wd_continuity_nxtcoef}
+  \begin{split}
+    \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
+    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} }
+  \end{split}
 \end{equation}
 
@@ -1541 +1460 @@
 directional limiter does.
 
-
-%----------------------------------------------------------------------------------------
 %      Surface pressure gradients
-%----------------------------------------------------------------------------------------
-\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]
-{Modification of surface pressure gradients (\mdl{dynhpg})}
-\label{subsubsec:DYN_wd_il_spg}
+%% =================================================================================================
+\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]{Modification of surface pressure gradients (\mdl{dynhpg})}
+\label{subsec:DYN_wd_il_spg}
 
 At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the
@@ -1560 +1476 @@
 neighbouring $(i+1,j)$ and $(i,j)$ tracer points.  zcpx is calculated using two logicals
 variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid
-column.  The three possible combinations are illustrated in figure \ref{Fig_WAD_dynhpg}.
-
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht] \begin{center}
-\includegraphics[width=\textwidth]{Fig_WAD_dynhpg}
-\caption{ \label{Fig_WAD_dynhpg}
-Illustrations of the three possible combinations of the logical variables controlling the
-limiting of the horizontal pressure gradient in wetting and drying regimes}
-\end{center}\end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+column.  The three possible combinations are illustrated in \autoref{fig:DYN_WAD_dynhpg}.
+
+\begin{figure}[!ht]
+  \centering
+  \includegraphics[width=0.66\textwidth]{Fig_WAD_dynhpg}
+  \caption[Combinations controlling the limiting of the horizontal pressure gradient in
+  wetting and drying regimes]{
+    Three possible combinations of the logical variables controlling the
+    limiting of the horizontal pressure gradient in wetting and drying regimes}
+  \label{fig:DYN_WAD_dynhpg}
+\end{figure}
 
 The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at
@@ -1576 +1493 @@
 of the topography at the two points:
 
-\begin{equation} \label{dyn_ll_tmp1}
-\begin{split}
-\mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\
+\begin{equation}
+  \label{eq:DYN_ll_tmp1}
+  \begin{split}
+    \mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\
                      & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\  .and.} \\
-& \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\
-& \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\
-& \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }
-\end{split}
+                     & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\
+                     & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\
+                     & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }
+  \end{split}
 \end{equation}
 
@@ -1590 +1508 @@
 at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$
 
-\begin{equation} \label{dyn_ll_tmp2}
-\begin{split}
-\mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\
-& \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\
-& \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .
-\end{split}
+\begin{equation}
+  \label{eq:DYN_ll_tmp2}
+  \begin{split}
+    \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\
+    & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\
+    & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .
+  \end{split}
 \end{equation}
 
@@ -1611 +1530 @@
 conditions.
 
-\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]
-{Additional considerations (\mdl{usrdef\_zgr})}
-\label{subsubsec:WAD_additional}
+%% =================================================================================================
+\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]{Additional considerations (\mdl{usrdef\_zgr})}
+\label{subsec:DYN_WAD_additional}
 
 In the very shallow water where wetting and drying occurs the parametrisation of
@@ -1623 +1542 @@
 in uncoupled integrations the net surface heat fluxes need to be appropriately limited.
 
-%----------------------------------------------------------------------------------------
 %      The WAD test cases
-%----------------------------------------------------------------------------------------
-\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]
-{The WAD test cases (\mdl{usrdef\_zgr})}
-\label{WAD_test_cases}
+%% =================================================================================================
+\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]{The WAD test cases (\mdl{usrdef\_zgr})}
+\label{subsec:DYN_WAD_test_cases}
 
 See the WAD tests MY\_DOC documention for details of the WAD test cases.
 
-
-
-% ================================================================
-% Time evolution term
-% ================================================================
-\section[Time evolution term (\textit{dynnxt.F90})]
-{Time evolution term (\protect\mdl{dynnxt})}
+%% =================================================================================================
+\section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})}
 \label{sec:DYN_nxt}
 
-%----------------------------------------------namdom----------------------------------------------------
-
-\nlst{namdom}
-%-------------------------------------------------------------------------------------------------------------
-
-Options are defined through the \nam{dom} namelist variables.
+Options are defined through the \nam{dom}{dom} namelist variables.
 The general framework for dynamics time stepping is a leap-frog scheme,
-\ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:STP}).
+\ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:TD}).
 The scheme is applied to the velocity, except when
 using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux})
 in the variable volume case (\texttt{vvl?} defined),
-where it has to be applied to the thickness weighted velocity (see \autoref{sec:A_momentum})
+where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum})
 
 $\bullet$ vector invariant form or linear free surface
-(\np{ln\_dynhpg\_vec}\forcode{ = .true.} ; \texttt{vvl?} not defined):
-\[
-  % \label{eq:dynnxt_vec}
+(\np[=.true.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} not defined):
+\[
+  % \label{eq:DYN_nxt_vec}
   \left\{
     \begin{aligned}
@@ -1667 +1574 @@
 
 $\bullet$ flux form and nonlinear free surface
-(\np{ln\_dynhpg\_vec}\forcode{ = .false.} ; \texttt{vvl?} defined):
-\[
-  % \label{eq:dynnxt_flux}
+(\np[=.false.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} defined):
+\[
+  % \label{eq:DYN_nxt_flux}
   \left\{
     \begin{aligned}
@@ -1680 +1587 @@
 where RHS is the right hand side of the momentum equation,
 the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient.
-$\gamma$ is initialized as \np{nn\_atfp} (namelist parameter).
-Its default value is \np{nn\_atfp}\forcode{ = 10.e-3}.
+$\gamma$ is initialized as \np{nn_atfp}{nn\_atfp} (namelist parameter).
+Its default value is \np[=10.e-3]{nn_atfp}{nn\_atfp}.
 In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for
 the momentum equations.
@@ -1689 +1596 @@
 and only array swapping and Asselin filtering is done in \mdl{dynnxt}.
 
-% ================================================================
-\biblio
-
-\pindex
+\onlyinsubfile{\input{../../global/epilogue}}
 
 \end{document}

@@ -1 @@
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
+*.ind
 *.ilg
-*.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*