Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex	(revision 11597)
@@ -8,4 +8,5 @@
 \vfill
 \begin{figure}[b]
+%% =================================================================================================
 \subsubsection*{Changes record}
 \begin{tabular}{m{0.08\linewidth}||m{0.32\linewidth}|m{0.6\linewidth}}
@@ -30,8 +31,8 @@
 of those described elsewhere in this manual.
 
+%% =================================================================================================
 \section{Choice of horizontal grid}
 \label{sec:DOMCFG_hor}
 
-%--------------------------------------------namdom-------------------------------------------------------
 
 \begin{listing}
@@ -40,5 +41,4 @@
   \label{lst:namdom_domcfg}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 The user has three options available in defining a horizontal grid, which involve the
@@ -94,7 +94,9 @@
 (and the number of grid points).
 
+%% =================================================================================================
 \section{Vertical grid}
 \label{sec:DOMCFG_vert}
 
+%% =================================================================================================
 \subsection{Vertical reference coordinate}
 \label{sec:DOMCFG_zref}
@@ -289,4 +291,5 @@
 %% %        Meter Bathymetry
 %% % -------------------------------------------------------------------------------------------------------------
+%% =================================================================================================
 \subsection{Model bathymetry}
 \label{subsec:DOMCFG_bathy}
@@ -317,4 +320,5 @@
 \end{description}
 
+%% =================================================================================================
 \subsection{Choice of vertical grid}
 \label{sec:DOMCFG_vgrd}
@@ -337,4 +341,5 @@
 %%%
 
+%% =================================================================================================
 \subsubsection[$Z$-coordinate with uniform thickness levels (\forcode{ln_zco})]{$Z$-coordinate with uniform thickness levels (\protect\np{ln_zco}{ln\_zco})}
 \label{subsec:DOMCFG_zco}
@@ -346,4 +351,5 @@
 rarely used in modern simulations but it can be useful for testing purposes.
 
+%% =================================================================================================
 \subsubsection[$Z$-coordinate with partial step (\forcode{ln_zps})]{$Z$-coordinate with partial step (\protect\np{ln_zps}{ln\_zps})}
 \label{subsec:DOMCFG_zps}
@@ -373,7 +379,7 @@
 the default thickness $e_{3t}(jk)$).
 
+%% =================================================================================================
 \subsubsection[$S$-coordinate (\forcode{ln_sco})]{$S$-coordinate (\protect\np{ln_sco}{ln\_sco})}
 \label{sec:DOMCFG_sco}
-%------------------------------------------nam_zgr_sco---------------------------------------------------
 %
 \begin{listing}
@@ -382,5 +388,4 @@
   \label{lst:namzgr_sco_domcfg}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 Options are defined in \nam{zgr_sco}{zgr\_sco} (\texttt{DOMAINcfg} only).
 In $s$-coordinate (\np[=.true.]{ln_sco}{ln\_sco}), the depth and thickness of the model levels are defined from
@@ -519,4 +524,5 @@
 and is output as part of the model mesh file at the start of the run.
 
+%% =================================================================================================
 \subsubsection[\zstar- or \sstar-coordinate (\forcode{ln_linssh})]{\zstar- or \sstar-coordinate (\protect\np{ln_linssh}{ln\_linssh})}
 \label{subsec:DOMCFG_zgr_star}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex	(revision 11597)
@@ -9,4 +9,5 @@
 This appendix some on going consideration on algorithms used or planned to be used in \NEMO.
 
+%% =================================================================================================
 \section{Upstream Biased Scheme (UBS) (\protect\np[=.true.]{ln_traadv_ubs}{ln\_traadv\_ubs})}
 \label{sec:ALGOS_tra_adv_ubs}
@@ -180,4 +181,5 @@
 which leads to ${A_u^{lT}} = \frac{1}{12} {e_{1u}}^3\ |u|$
 
+%% =================================================================================================
 \section{Leapfrog energetic}
 \label{sec:ALGOS_LF}
@@ -233,6 +235,8 @@
 In time this boundary condition is not physical and \textbf{add something here!!!}
 
+%% =================================================================================================
 \section{Lateral diffusion operator}
 
+%% =================================================================================================
 \subsection{Griffies iso-neutral diffusion operator}
 
@@ -426,4 +430,5 @@
 \end{description}
 
+%% =================================================================================================
 \subsection{Eddy induced velocity and skew flux formulation}
 
@@ -579,4 +584,5 @@
 
 $\ $\newpage      %force an empty line
+%% =================================================================================================
 \subsection{Discrete invariants of the iso-neutral diffrusion}
 \label{subsec:ALGOS_Gf_operator}
@@ -741,4 +747,5 @@
 There is no need to develop a specific to obtain it.
 
+%% =================================================================================================
 \subsection{Discrete invariants of the skew flux formulation}
 \label{subsec:ALGOS_eiv_skew}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_diff_opers.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_diff_opers.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_diff_opers.tex	(revision 11597)
@@ -7,7 +7,9 @@
 \chaptertoc
 
+%% =================================================================================================
 \section{Horizontal/Vertical $2^{nd}$ order tracer diffusive operators}
 \label{sec:DIFFOPERS_1}
 
+%% =================================================================================================
 \subsubsection*{In z-coordinates}
 
@@ -22,4 +24,5 @@
 \end{align}
 
+%% =================================================================================================
 \subsubsection*{In generalized vertical coordinates}
 
@@ -148,7 +151,9 @@
 %\addtocounter{equation}{-2}
 
+%% =================================================================================================
 \section{Iso/Diapycnal $2^{nd}$ order tracer diffusive operators}
 \label{sec:DIFFOPERS_2}
 
+%% =================================================================================================
 \subsubsection*{In z-coordinates}
 
@@ -252,4 +257,5 @@
 the property becomes obvious.
 
+%% =================================================================================================
 \subsubsection*{In generalized vertical coordinates}
 
@@ -309,4 +315,5 @@
 \autoref{sec:DIFFOPERS_1} onto $s$-coordinates is exact, however steep the $s$-surfaces.
 
+%% =================================================================================================
 \section{Lateral/Vertical momentum diffusive operators}
 \label{sec:DIFFOPERS_3}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex	(revision 11597)
@@ -12,4 +12,5 @@
 %\gmcomment{
 
+%% =================================================================================================
 \section{Introduction / Notations}
 \label{sec:INVARIANTS_0}
@@ -85,4 +86,5 @@
 \end{flalign}
 
+%% =================================================================================================
 \section{Continuous conservation}
 \label{sec:INVARIANTS_1}
@@ -311,7 +313,9 @@
 %
 
+%% =================================================================================================
 \section{Discrete total energy conservation: vector invariant form}
 \label{sec:INVARIANTS_2}
 
+%% =================================================================================================
 \subsection{Total energy conservation}
 \label{subsec:INVARIANTS_KE+PE_vect}
@@ -337,4 +341,5 @@
 leads to the discrete equivalent of the four equations \autoref{eq:INVARIANTS_E_tot_flux}.
 
+%% =================================================================================================
 \subsection{Vorticity term (coriolis + vorticity part of the advection)}
 \label{subsec:INVARIANTS_vor}
@@ -343,4 +348,5 @@
 or the planetary ($q=f/e_{3f}$), or the total potential vorticity ($q=(\zeta +f) /e_{3f}$).
 Two discretisation of the vorticity term (ENE and EEN) allows the conservation of the kinetic energy.
+%% =================================================================================================
 \subsubsection{Vorticity term with ENE scheme (\protect\np[=.true.]{ln_dynvor_ene}{ln\_dynvor\_ene})}
 \label{subsec:INVARIANTS_vorENE}
@@ -380,4 +386,5 @@
 In other words, the domain averaged kinetic energy does not change due to the vorticity term.
 
+%% =================================================================================================
 \subsubsection{Vorticity term with EEN scheme (\protect\np[=.true.]{ln_dynvor_een}{ln\_dynvor\_een})}
 \label{subsec:INVARIANTS_vorEEN_vect}
@@ -449,4 +456,5 @@
 \end{flalign*}
 
+%% =================================================================================================
 \subsubsection{Gradient of kinetic energy / Vertical advection}
 \label{subsec:INVARIANTS_zad}
@@ -556,4 +564,5 @@
 Blah blah required on the the step representation of bottom topography.....
 
+%% =================================================================================================
 \subsection{Pressure gradient term}
 \label{subsec:INVARIANTS_2.6}
@@ -698,7 +707,9 @@
 Nevertheless, it is almost never satisfied since a linear equation of state is rarely used.
 
+%% =================================================================================================
 \section{Discrete total energy conservation: flux form}
 \label{sec:INVARIANTS_3}
 
+%% =================================================================================================
 \subsection{Total energy conservation}
 \label{subsec:INVARIANTS_KE+PE_flux}
@@ -721,7 +732,9 @@
 vector invariant or in flux form, leads to the discrete equivalent of the ????
 
+%% =================================================================================================
 \subsection{Coriolis and advection terms: flux form}
 \label{subsec:INVARIANTS_3.2}
 
+%% =================================================================================================
 \subsubsection{Coriolis plus ``metric'' term}
 \label{subsec:INVARIANTS_3.3}
@@ -742,4 +755,5 @@
 The derivation is the same as for the vorticity term in the vector invariant form (\autoref{subsec:INVARIANTS_vor}).
 
+%% =================================================================================================
 \subsubsection{Flux form advection}
 \label{subsec:INVARIANTS_3.4}
@@ -820,7 +834,9 @@
 The horizontal kinetic energy is not conserved, but forced to decay (\ie\ the scheme is diffusive).
 
+%% =================================================================================================
 \section{Discrete enstrophy conservation}
 \label{sec:INVARIANTS_4}
 
+%% =================================================================================================
 \subsubsection{Vorticity term with ENS scheme  (\protect\np[=.true.]{ln_dynvor_ens}{ln\_dynvor\_ens})}
 \label{subsec:INVARIANTS_vorENS}
@@ -889,4 +905,5 @@
 The later equality is obtain only when the flow is horizontally non-divergent, \ie\ $\chi$=$0$.
 
+%% =================================================================================================
 \subsubsection{Vorticity Term with EEN scheme (\protect\np[=.true.]{ln_dynvor_een}{ln\_dynvor\_een})}
 \label{subsec:INVARIANTS_vorEEN}
@@ -959,4 +976,5 @@
 \end{flalign*}
 
+%% =================================================================================================
 \section{Conservation properties on tracers}
 \label{sec:INVARIANTS_5}
@@ -972,4 +990,5 @@
 as the equation of state is non linear with respect to $T$ and $S$.
 In practice, the mass is conserved to a very high accuracy.
+%% =================================================================================================
 \subsection{Advection term}
 \label{subsec:INVARIANTS_5.1}
@@ -1035,4 +1054,5 @@
 which is the discrete form of $ \frac{1}{2} \int_D {  T^2 \frac{1}{e_3} \frac{\partial  e_3 }{\partial t} \;dv }$.
 
+%% =================================================================================================
 \section{Conservation properties on lateral momentum physics}
 \label{sec:INVARIANTS_dynldf_properties}
@@ -1053,4 +1073,5 @@
 the term associated with the horizontal gradient of the divergence is locally zero.
 
+%% =================================================================================================
 \subsection{Conservation of potential vorticity}
 \label{subsec:INVARIANTS_6.1}
@@ -1084,4 +1105,5 @@
 \end{flalign*}
 
+%% =================================================================================================
 \subsection{Dissipation of horizontal kinetic energy}
 \label{subsec:INVARIANTS_6.2}
@@ -1133,4 +1155,5 @@
 \]
 
+%% =================================================================================================
 \subsection{Dissipation of enstrophy}
 \label{subsec:INVARIANTS_6.3}
@@ -1154,4 +1177,5 @@
 \end{flalign*}
 
+%% =================================================================================================
 \subsection{Conservation of horizontal divergence}
 \label{subsec:INVARIANTS_6.4}
@@ -1178,4 +1202,5 @@
 \end{flalign*}
 
+%% =================================================================================================
 \subsection{Dissipation of horizontal divergence variance}
 \label{subsec:INVARIANTS_6.5}
@@ -1201,4 +1226,5 @@
 \end{flalign*}
 
+%% =================================================================================================
 \section{Conservation properties on vertical momentum physics}
 \label{sec:INVARIANTS_7}
@@ -1369,4 +1395,5 @@
 \end{flalign*}
 
+%% =================================================================================================
 \section{Conservation properties on tracer physics}
 \label{sec:INVARIANTS_8}
@@ -1378,4 +1405,5 @@
 As for the advection term, there is conservation of mass only if the Equation Of Seawater is linear.
 
+%% =================================================================================================
 \subsection{Conservation of tracers}
 \label{subsec:INVARIANTS_8.1}
@@ -1408,4 +1436,5 @@
 In fact, this property simply results from the flux form of the operator.
 
+%% =================================================================================================
 \subsection{Dissipation of tracer variance}
 \label{subsec:INVARIANTS_8.2}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_s_coord.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_s_coord.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_s_coord.tex	(revision 11597)
@@ -10,4 +10,5 @@
 \vfill
 \begin{figure}[b]
+%% =================================================================================================
 \subsubsection*{Changes record}
 \begin{tabular}{l||l|m{0.65\linewidth}}
@@ -18,4 +19,5 @@
 \end{figure}
 
+%% =================================================================================================
 \section{Chain rule for $s-$coordinates}
 \label{sec:SCOORD_chain}
@@ -112,4 +114,5 @@
 \end{equation}
 
+%% =================================================================================================
 \section{Continuity equation in $s-$coordinates}
 \label{sec:SCOORD_continuity}
@@ -227,4 +230,5 @@
 the contribution of the time variation of the vertical coordinate to the volume budget.
 
+%% =================================================================================================
 \section{Momentum equation in $s-$coordinate}
 \label{sec:SCOORD_momentum}
@@ -554,4 +558,5 @@
 \ie\ the volume flux across the moving $s$-surfaces per unit horizontal area.
 
+%% =================================================================================================
 \section{Tracer equation}
 \label{sec:SCOORD_tracer}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex	(revision 11597)
@@ -17,8 +17,7 @@
 \chaptertoc
 
+%% =================================================================================================
 \section[Choice of \forcode{namtra\_ldf} namelist parameters]{Choice of \protect\nam{tra_ldf}{tra\_ldf} namelist parameters}
-%-----------------------------------------nam_traldf------------------------------------------------------
-
-%---------------------------------------------------------------------------------------------------------
+
 
 Two scheme are available to perform the iso-neutral diffusion.
@@ -63,4 +62,5 @@
 \end{description}
 
+%% =================================================================================================
 \section{Triad formulation of iso-neutral diffusion}
 \label{sec:TRIADS_iso}
@@ -69,4 +69,5 @@
 but formulated within the \NEMO\ framework, using scale factors rather than grid-sizes.
 
+%% =================================================================================================
 \subsection{Iso-neutral diffusion operator}
 
@@ -156,4 +157,5 @@
 is evaluated at $w$-points but involves horizontal gradients defined at $u$-points.
 
+%% =================================================================================================
 \subsection{Standard discretization}
 
@@ -187,4 +189,5 @@
 (\ie\ they enter the computation of density), but it does not work for a passive tracer.
 
+%% =================================================================================================
 \subsection{Expression of the skew-flux in terms of triad slopes}
 
@@ -282,4 +285,5 @@
 and in \autoref{eq:TRIADS_i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.
 
+%% =================================================================================================
 \subsection{Full triad fluxes}
 
@@ -365,4 +369,5 @@
 \end{flalign}
 
+%% =================================================================================================
 \subsection{Ensuring the scheme does not increase tracer variance}
 \label{subsec:TRIADS_variance}
@@ -462,4 +467,5 @@
 \]
 
+%% =================================================================================================
 \subsection{Triad volumes in Griffes's scheme and in \NEMO}
 
@@ -490,4 +496,5 @@
 we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above.
 
+%% =================================================================================================
 \subsection{Summary of the scheme}
 
@@ -625,4 +632,5 @@
 \end{description}
 
+%% =================================================================================================
 \subsection{Treatment of the triads at the boundaries}
 \label{sec:TRIADS_iso_bdry}
@@ -673,4 +681,5 @@
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
+%% =================================================================================================
 \subsection{ Limiting of the slopes within the interior}
 \label{sec:TRIADS_limit}
@@ -701,4 +710,5 @@
 and so acts to reduce gravitational potential energy.
 
+%% =================================================================================================
 \subsection{Tapering within the surface mixed layer}
 \label{sec:TRIADS_taper}
@@ -707,4 +717,5 @@
 When the Griffies triads are used, we offer two options for this.
 
+%% =================================================================================================
 \subsubsection{Linear slope tapering within the surface mixed layer}
 \label{sec:TRIADS_lintaper}
@@ -822,4 +833,5 @@
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
+%% =================================================================================================
 \subsubsection{Additional truncation of skew iso-neutral flux components}
 \label{subsec:TRIADS_Gerdes-taper}
@@ -864,7 +876,9 @@
 % This may give strange looking results,
 % particularly where the mixed-layer depth varies strongly laterally.
+%% =================================================================================================
 \section{Eddy induced advection formulated as a skew flux}
 \label{sec:TRIADS_skew-flux}
 
+%% =================================================================================================
 \subsection{Continuous skew flux formulation}
 \label{sec:TRIADS_continuous-skew-flux}
@@ -975,4 +989,5 @@
 Since it has the same divergence as the advective form it also preserves the tracer variance.
 
+%% =================================================================================================
 \subsection{Discrete skew flux formulation}
 
@@ -1017,4 +1032,5 @@
 It also ensures the following two key properties.
 
+%% =================================================================================================
 \subsubsection{No change in tracer variance}
 
@@ -1039,4 +1055,5 @@
 Hence the two fluxes associated with each triad make no net contribution to the variance budget.
 
+%% =================================================================================================
 \subsubsection{Reduction in gravitational PE}
 
@@ -1093,4 +1110,5 @@
 \beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative.
 
+%% =================================================================================================
 \subsection{Treatment of the triads at the boundaries}
 \label{sec:TRIADS_skew_bdry}
@@ -1104,4 +1122,5 @@
 The namelist parameter \np{ln_botmix_triad}{ln\_botmix\_triad} has no effect on the eddy-induced skew-fluxes.
 
+%% =================================================================================================
 \subsection{Limiting of the slopes within the interior}
 \label{sec:TRIADS_limitskew}
@@ -1111,4 +1130,5 @@
 Each individual triad \rtriadt{R} is so limited.
 
+%% =================================================================================================
 \subsection{Tapering within the surface mixed layer}
 \label{sec:TRIADS_taperskew}
@@ -1131,4 +1151,5 @@
 (the horizontal flux convergence is relatively insignificant within the mixed-layer).
 
+%% =================================================================================================
 \subsection{Streamfunction diagnostics}
 \label{sec:TRIADS_sfdiag}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex	(revision 11597)
@@ -9,4 +9,5 @@
 \vfill
 \begin{figure}[b]
+%% =================================================================================================
 \subsubsection*{Changes record}
 \begin{tabular}{l||l|m{0.65\linewidth}}
@@ -27,4 +28,5 @@
 %===============================================================
 
+%% =================================================================================================
 \section{Direct initialization}
 \label{sec:ASM_DI}
@@ -34,4 +36,5 @@
 DI is used when \np{ln_asmdin}{ln\_asmdin} is set to true.
 
+%% =================================================================================================
 \section{Incremental analysis updates}
 \label{sec:ASM_IAU}
@@ -92,4 +95,5 @@
 %==========================================================================
 % Divergence damping description %%%
+%% =================================================================================================
 \section{Divergence damping initialisation}
 \label{sec:ASM_div_dmp}
@@ -135,4 +139,5 @@
 %==========================================================================
 
+%% =================================================================================================
 \section{Implementation details}
 \label{sec:ASM_details}
@@ -141,5 +146,4 @@
 the ORCA2 grid.
 
-%------------------------------------------nam_asminc-----------------------------------------------------
 %
 \begin{listing}
@@ -148,5 +152,4 @@
   \label{lst:nam_asminc}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 The header of an assimilation increments file produced using the NetCDF tool
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex	(revision 11597)
@@ -9,4 +9,5 @@
 \vfill
 \begin{figure}[b]
+%% =================================================================================================
 \subsubsection*{Changes record}
 \begin{tabular}{l||l|m{0.65\linewidth}}
@@ -20,4 +21,5 @@
 \end{figure}
 
+%% =================================================================================================
 \section{Model output}
 \label{sec:DIA_io_old}
@@ -45,4 +47,5 @@
 %\gmcomment{                    % start of gmcomment
 
+%% =================================================================================================
 \section{Standard model output (IOM)}
 \label{sec:DIA_iom}
@@ -103,4 +106,5 @@
 even without a parallel-enabled NetCDF4 library, simply by employing only one dedicated I/O server.
 
+%% =================================================================================================
 \subsection{XIOS: Reading and writing restart file}
 
@@ -142,6 +146,8 @@
 An older versions of XIOS do not support reading functionality. It's recommended to use at least XIOS2@1451.
 
+%% =================================================================================================
 \subsection{XIOS: XML Inputs-Outputs Server}
 
+%% =================================================================================================
 \subsubsection{Attached or detached mode?}
 
@@ -178,4 +184,5 @@
 The following subsection provides a typical example but the syntax will vary in different MPP environments.
 
+%% =================================================================================================
 \subsubsection{Number of cpu used by XIOS in detached mode}
 
@@ -192,4 +199,5 @@
 \cmd|mpirun -np 62 ./nemo.exe : -np 2 ./xios_server.exe|
 
+%% =================================================================================================
 \subsubsection{Control of XIOS: the context in iodef.xml}
 
@@ -237,6 +245,8 @@
 \end{table}
 
+%% =================================================================================================
 \subsection{Practical issues}
 
+%% =================================================================================================
 \subsubsection{Installation}
 
@@ -247,4 +257,5 @@
 {Extract and install XIOS} guide provides an example illustration of how this can be achieved.
 
+%% =================================================================================================
 \subsubsection{Add your own outputs}
 
@@ -303,6 +314,8 @@
 \end{enumerate}
 
+%% =================================================================================================
 \subsection{XML fundamentals}
 
+%% =================================================================================================
 \subsubsection{ XML basic rules}
 
@@ -314,4 +327,5 @@
 See \href{http://www.xmlnews.org/docs/xml-basics.html}{here} for more details.
 
+%% =================================================================================================
 \subsubsection{Structure of the XML file used in \NEMO}
 
@@ -419,4 +433,5 @@
 \end{table}
 
+%% =================================================================================================
 \subsubsection{Nesting XML files}
 
@@ -435,4 +450,5 @@
 \end{xmllines}
 
+%% =================================================================================================
 \subsubsection{Use of inheritance}
 
@@ -475,4 +491,5 @@
 Inherit (and overwrite, if needed) the attributes of a tag you are refering to.
 
+%% =================================================================================================
 \subsubsection{Use of groups}
 
@@ -516,4 +533,5 @@
 \end{xmllines}
 
+%% =================================================================================================
 \subsection{Detailed functionalities}
 
@@ -521,4 +539,5 @@
 the new functionalities offered by the XML interface of XIOS.
 
+%% =================================================================================================
 \subsubsection{Define horizontal subdomains}
 
@@ -562,4 +581,5 @@
 We are therefore advising to use the ''one\_file'' type in this case.
 
+%% =================================================================================================
 \subsubsection{Define vertical zooms}
 
@@ -585,4 +605,5 @@
 \end{xmllines}
 
+%% =================================================================================================
 \subsubsection{Control of the output file names}
 
@@ -657,4 +678,5 @@
 \noindent will give the following file name radical: \ifile{myfile\_ORCA2\_19891231\_freq1d}
 
+%% =================================================================================================
 \subsubsection{Other controls of the XML attributes from \NEMO}
 
@@ -710,6 +732,8 @@
 \end{table}
 
+%% =================================================================================================
 \subsubsection{Advanced use of XIOS functionalities}
 
+%% =================================================================================================
 \subsection{XML reference tables}
 \label{subsec:DIA_IOM_xmlref}
@@ -858,4 +882,5 @@
 field\_ref="sst" means that attributes not explicitely defined, are inherited from sst field.
 
+%% =================================================================================================
 \subsubsection{Tag list per family}
 
@@ -1051,4 +1076,5 @@
 \end{table}
 
+%% =================================================================================================
 \subsubsection{Attributes list per family}
 
@@ -1286,4 +1312,5 @@
 \end{table}
 
+%% =================================================================================================
 \subsection{CF metadata standard compliance}
 
@@ -1298,4 +1325,5 @@
 This must be set to true if these metadata are to be included in the output files.
 
+%% =================================================================================================
 \section[NetCDF4 support (\texttt{\textbf{key\_netcdf4}})]{NetCDF4 support (\protect\key{netcdf4})}
 \label{sec:DIA_nc4}
@@ -1315,5 +1343,4 @@
 setting the \np{ln_nc4zip}{ln\_nc4zip} logical to false in the \nam{nc4}{nc4} namelist:
 
-%------------------------------------------namnc4----------------------------------------------------
 
 \begin{listing}
@@ -1322,5 +1349,4 @@
   \label{lst:namnc4}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 If \key{netcdf4} has not been defined, these namelist parameters are not read.
@@ -1367,5 +1393,4 @@
 each processing region.
 
-%------------------------------------------TABLE----------------------------------------------------
 \begin{table}
   \centering
@@ -1403,5 +1428,4 @@
   \label{tab:DIA_NC4}
 \end{table}
-%----------------------------------------------------------------------------------------------------
 
 When \key{iomput} is activated with \key{netcdf4} chunking and compression parameters for fields produced via
@@ -1414,8 +1438,8 @@
 the invidual processing regions and different chunking choices may be desired.
 
+%% =================================================================================================
 \section[Tracer/Dynamics trends (\forcode{&namtrd})]{Tracer/Dynamics trends (\protect\nam{trd}{trd})}
 \label{sec:DIA_trd}
 
-%------------------------------------------namtrd----------------------------------------------------
 
 \begin{listing}
@@ -1424,5 +1448,4 @@
   \label{lst:namtrd}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 Each trend of the dynamics and/or temperature and salinity time evolution equations can be send to
@@ -1462,7 +1485,7 @@
 and none of the options have been tested with variable volume (\ie\ \np[=.true.]{ln_linssh}{ln\_linssh}).
 
+%% =================================================================================================
 \section[FLO: On-Line Floats trajectories (\texttt{\textbf{key\_floats}})]{FLO: On-Line Floats trajectories (\protect\key{floats})}
 \label{sec:DIA_FLO}
-%--------------------------------------------namflo-------------------------------------------------------
 
 \begin{listing}
@@ -1471,5 +1494,4 @@
   \label{lst:namflo}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 The on-line computation of floats advected either by the three dimensional velocity field or constraint to
@@ -1481,4 +1503,5 @@
 are consistent with the numeric of the code, so that the trajectories never intercept the bathymetry.
 
+%% =================================================================================================
 \subsubsection{Input data: initial coordinates}
 
@@ -1538,4 +1561,5 @@
 \np{jpnflnewflo}{jpnflnewflo} can be added in the initialization file.
 
+%% =================================================================================================
 \subsubsection{Output data}
 
@@ -1564,8 +1588,8 @@
 \end{xmllines}
 
+%% =================================================================================================
 \section[Harmonic analysis of tidal constituents (\texttt{\textbf{key\_diaharm}})]{Harmonic analysis of tidal constituents (\protect\key{diaharm})}
 \label{sec:DIA_diag_harm}
 
-%------------------------------------------nam_diaharm----------------------------------------------------
 %
 \begin{listing}
@@ -1574,5 +1598,4 @@
   \label{lst:nam_diaharm}
 \end{listing}
-%----------------------------------------------------------------------------------------------------------
 
 A module is available to compute the amplitude and phase of tidal waves.
@@ -1612,8 +1635,8 @@
 We obtain in output $C_{j}$ and $S_{j}$ for each tidal wave.
 
+%% =================================================================================================
 \section[Transports across sections (\texttt{\textbf{key\_diadct}})]{Transports across sections (\protect\key{diadct})}
 \label{sec:DIA_diag_dct}
 
-%------------------------------------------nam_diadct----------------------------------------------------
 
 \begin{listing}
@@ -1622,5 +1645,4 @@
   \label{lst:nam_diadct}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 A module is available to compute the transport of volume, heat and salt through sections.
@@ -1646,4 +1668,5 @@
 \np{nn_debug}{nn\_debug} : debugging of the section
 
+%% =================================================================================================
 \subsubsection{Creating a binary file containing the pathway of each section}
 
@@ -1717,4 +1740,5 @@
  }
 
+%% =================================================================================================
 \subsubsection{To read the output files}
 
@@ -1755,4 +1779,5 @@
 \end{table}
 
+%% =================================================================================================
 \section{Diagnosing the steric effect in sea surface height}
 \label{sec:DIA_steric}
@@ -1931,4 +1956,5 @@
 Both steric and thermosteric sea level are computed in \mdl{diaar5}.
 
+%% =================================================================================================
 \section{Other diagnostics}
 \label{sec:DIA_diag_others}
@@ -1937,4 +1963,5 @@
 The available ready-to-add diagnostics modules can be found in directory DIA.
 
+%% =================================================================================================
 \subsection[Depth of various quantities (\textit{diahth.F90})]{Depth of various quantities (\protect\mdl{diahth})}
 
@@ -1969,4 +1996,5 @@
 %       CMIP specific diagnostics
 % -----------------------------------------------------------
+%% =================================================================================================
 \subsection[CMIP specific diagnostics (\textit{diaar5.F90}, \textit{diaptr.F90})]{CMIP specific diagnostics (\protect\mdl{diaar5})}
 
@@ -1987,5 +2015,4 @@
 the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:DIA_mask_subasins}).
 
-%------------------------------------------namptr-----------------------------------------
 
 \begin{listing}
@@ -1994,12 +2021,11 @@
   \label{lst:namptr}
 \end{listing}
-%-----------------------------------------------------------------------------------------
 
 % -----------------------------------------------------------
 %       25 hour mean and hourly Surface, Mid and Bed
 % -----------------------------------------------------------
+%% =================================================================================================
 \subsection{25 hour mean output for tidal models}
 
-%------------------------------------------nam_dia25h-------------------------------------
 
 \begin{listing}
@@ -2008,5 +2034,4 @@
   \label{lst:nam_dia25h}
 \end{listing}
-%-----------------------------------------------------------------------------------------
 
 A module is available to compute a crudely detided M2 signal by obtaining a 25 hour mean.
@@ -2018,7 +2043,7 @@
 %     Top Middle and Bed hourly output
 % -----------------------------------------------------------
+%% =================================================================================================
 \subsection{Top middle and bed hourly output}
 
-%------------------------------------------nam_diatmb-----------------------------------------------------
 
 \begin{listing}
@@ -2027,5 +2052,4 @@
   \label{lst:nam_diatmb}
 \end{listing}
-%----------------------------------------------------------------------------------------------------------
 
 A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables.
@@ -2038,4 +2062,5 @@
 %     Courant numbers
 % -----------------------------------------------------------
+%% =================================================================================================
 \subsection{Courant numbers}
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex	(revision 11597)
@@ -52,4 +52,5 @@
 
 %===============================================================
+%% =================================================================================================
 \section{Warm layer model}
 \label{sec:DIU_warm_layer_sec}
@@ -109,4 +110,5 @@
 %===============================================================
 
+%% =================================================================================================
 \section{Cool skin model}
 \label{sec:DIU_cool_skin_sec}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex	(revision 11597)
@@ -38,7 +38,9 @@
 and other relevant information about the DOM (DOMain) source code modules.
 
+%% =================================================================================================
 \section{Fundamentals of the discretisation}
 \label{sec:DOM_basics}
 
+%% =================================================================================================
 \subsection{Arrangement of variables}
 \label{subsec:DOM_cell}
@@ -145,4 +147,5 @@
 \end{figure}
 
+%% =================================================================================================
 \subsection{Discrete operators}
 \label{subsec:DOM_operators}
@@ -235,4 +238,5 @@
 demonstrate integral conservative properties of the discrete formulation chosen.
 
+%% =================================================================================================
 \subsection{Numerical indexing}
 \label{subsec:DOM_Num_Index}
@@ -258,4 +262,5 @@
 %        Horizontal Indexing
 % -----------------------------------
+%% =================================================================================================
 \subsubsection{Horizontal indexing}
 \label{subsec:DOM_Num_Index_hor}
@@ -270,4 +275,5 @@
 %        Vertical indexing
 % -----------------------------------
+%% =================================================================================================
 \subsubsection{Vertical indexing}
 \label{subsec:DOM_Num_Index_vertical}
@@ -301,4 +307,5 @@
 \end{figure}
 
+%% =================================================================================================
 \section{Spatial domain configuration}
 \label{subsec:DOM_config}
@@ -338,4 +345,5 @@
 %        Domain Size
 % -----------------------------------
+%% =================================================================================================
 \subsection{Domain size}
 \label{subsec:DOM_size}
@@ -355,7 +363,9 @@
 See \autoref{sec:LBC_jperio} for details on the available options and the corresponding values for \jp{jperio}.
 
+%% =================================================================================================
 \subsection[Horizontal grid mesh (\textit{domhgr.F90}]{Horizontal grid mesh (\protect\mdl{domhgr})}
 \label{subsec:DOM_hgr}
 
+%% =================================================================================================
 \subsubsection{Required fields}
 \label{sec:DOM_hgr_fields}
@@ -380,4 +390,5 @@
 evaluated for the same arguments as $\lambda$ and $\varphi$.
 
+%% =================================================================================================
 \subsubsection{Optional fields}
 
@@ -417,7 +428,7 @@
 thus no specific arrays are defined at $w$ points.
 
+%% =================================================================================================
 \subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})}
 \label{subsec:DOM_zgr}
-%-----------------------------------------namdom-------------------------------------------
 \begin{listing}
   \nlst{namdom}
@@ -425,5 +436,4 @@
   \label{lst:namdom}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 In the vertical, the model mesh is determined by four things:
@@ -506,4 +516,5 @@
 their reference counterpart and remain fixed.
 
+%% =================================================================================================
 \subsubsection{Required fields}
 \label{sec:DOM_zgr_fields}
@@ -541,4 +552,5 @@
 With ice cavities, \jp{top\_level} determines the first wet point below the overlying ice shelf.
 
+%% =================================================================================================
 \subsubsection{Level bathymetry and mask}
 \label{subsec:DOM_msk}
@@ -572,7 +584,6 @@
 %% (see \autoref{fig:LBC_jperio}).
 
-%-------------------------------------------------------------------------------------------------
 %        Closed seas
-%-------------------------------------------------------------------------------------------------
+%% =================================================================================================
 \subsection{Closed seas}
 \label{subsec:DOM_closea}
@@ -595,4 +606,5 @@
 \end{clines}
 
+%% =================================================================================================
 \subsection{Output grid files}
 \label{subsec:DOM_meshmask}
@@ -613,7 +625,7 @@
 This file contains additional fields that can be useful for post-processing applications.
 
+%% =================================================================================================
 \section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})}
 \label{sec:DOM_DTA_tsd}
-%-----------------------------------------namtsd-------------------------------------------
 \begin{listing}
   \nlst{namtsd}
@@ -621,5 +633,4 @@
   \label{lst:namtsd}
 \end{listing}
-%------------------------------------------------------------------------------------------
 
 Basic initial state options are defined in \nam{tsd}{tsd}.
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex	(revision 11597)
@@ -53,10 +53,9 @@
 MISC correspond to "extracting tendency terms" or "vorticity balance"?}
 
+%% =================================================================================================
 \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})}
 \label{subsec:DYN_divcur}
@@ -92,7 +91,5 @@
 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})}
 \label{subsec:DYN_sshwzv}
@@ -152,7 +149,8 @@
 (see \autoref{subsec:DOM_Num_Index_vertical}).
 
+
+%% =================================================================================================
 \section{Coriolis and advection: vector invariant form}
 \label{sec:DYN_adv_cor_vect}
-%-----------------------------------------nam_dynadv----------------------------------------------------
 
 \begin{listing}
@@ -161,5 +159,4 @@
   \label{lst:namdyn_adv}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 The vector invariant form of the momentum equations is the one most often used in
@@ -172,7 +169,7 @@
 \autoref{chap:LBC}.
 
+%% =================================================================================================
 \subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})}
 \label{subsec:DYN_vor}
-%------------------------------------------nam_dynvor----------------------------------------------------
 
 \begin{listing}
@@ -181,5 +178,4 @@
   \label{lst:namdyn_vor}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{dyn_vor}{dyn\_vor} namelist variables.
@@ -195,7 +191,6 @@
 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})]{Enstrophy conserving scheme (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens})}
 \label{subsec:DYN_vor_ens}
@@ -218,7 +213,6 @@
 \end{equation}
 
-%-------------------------------------------------------------
 %                 energy conserving scheme
-%-------------------------------------------------------------
+%% =================================================================================================
 \subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene})]{Energy conserving scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene})}
 \label{subsec:DYN_vor_ene}
@@ -238,7 +232,6 @@
 \end{equation}
 
-%-------------------------------------------------------------
 %                 mix energy/enstrophy conserving scheme
-%-------------------------------------------------------------
+%% =================================================================================================
 \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}
@@ -263,7 +256,6 @@
 \]
 
-%-------------------------------------------------------------
 %                 energy and enstrophy conserving scheme
-%-------------------------------------------------------------
+%% =================================================================================================
 \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}
@@ -353,7 +345,5 @@
 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})}
 \label{subsec:DYN_keg}
@@ -373,7 +363,5 @@
 \]
 
-%--------------------------------------------------------------------------------------------------------------
-%           Vertical advection term
-%--------------------------------------------------------------------------------------------------------------
+%% =================================================================================================
 \subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})}
 \label{subsec:DYN_zad}
@@ -400,9 +388,8 @@
 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----------------------------------------------------
-
-%-------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables.
@@ -413,7 +400,5 @@
 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})}
 \label{subsec:DYN_cor_flux}
@@ -434,7 +419,5 @@
 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})}
 \label{subsec:DYN_adv_flux}
@@ -467,7 +450,6 @@
 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})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln_dynadv_cen2}{ln\_dynadv\_cen2})}
 \label{subsec:DYN_adv_cen2}
@@ -490,7 +472,6 @@
 so $u$ and $v$ are the \emph{now} velocities.
 
-%-------------------------------------------------------------
 %                 UBS scheme
-%-------------------------------------------------------------
+%% =================================================================================================
 \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}
@@ -543,7 +524,7 @@
 %%%
 
+%% =================================================================================================
 \section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
 \label{sec:DYN_hpg}
-%------------------------------------------nam_dynhpg---------------------------------------------------
 
 \begin{listing}
@@ -552,5 +533,4 @@
   \label{lst:namdyn_hpg}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{dyn_hpg}{dyn\_hpg} namelist variables.
@@ -566,7 +546,5 @@
 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})]{Full step $Z$-coordinate (\protect\np{ln_dynhpg_zco}{ln\_dynhpg\_zco})}
 \label{subsec:DYN_hpg_zco}
@@ -611,7 +589,5 @@
 \autoref{eq:DYN_hpg_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})]{Partial step $Z$-coordinate (\protect\np{ln_dynhpg_zps}{ln\_dynhpg\_zps})}
 \label{subsec:DYN_hpg_zps}
@@ -632,7 +608,5 @@
 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}
@@ -681,4 +655,5 @@
 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}
@@ -697,7 +672,5 @@
 \autoref{subsec:DYN_hpg_sco}.
 
-%--------------------------------------------------------------------------------------------------------------
-%           Time-scheme
-%--------------------------------------------------------------------------------------------------------------
+%% =================================================================================================
 \subsection[Time-scheme (\forcode{ln_dynhpg_imp})]{Time-scheme (\protect\np{ln_dynhpg_imp}{ln\_dynhpg\_imp})}
 \label{subsec:DYN_hpg_imp}
@@ -758,7 +731,7 @@
 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----------------------------------------------------
 
 \begin{listing}
@@ -767,5 +740,4 @@
   \label{lst:namdyn_spg}
 \end{listing}
-%------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables.
@@ -795,7 +767,5 @@
 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 (\forcode{ln_dynspg_exp})]{Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})}
 \label{subsec:DYN_spg_exp}
@@ -821,13 +791,9 @@
 Thus, nothing is done in the \mdl{dynspg\_exp} module.
 
-%--------------------------------------------------------------------------------------------------------------
-% Split-explict free surface formulation
-%--------------------------------------------------------------------------------------------------------------
+%% =================================================================================================
 \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}{ln\_dynspg\_ts} set to true),
@@ -1065,7 +1031,5 @@
 %>>>>>===============
 
-%--------------------------------------------------------------------------------------------------------------
-% Filtered free surface formulation
-%--------------------------------------------------------------------------------------------------------------
+%% =================================================================================================
 \subsection{Filtered free surface (\forcode{dynspg_flt?})}
 \label{subsec:DYN_spg_fltp}
@@ -1093,7 +1057,7 @@
 It is computed once and for all and applies to all ocean time steps.
 
+%% =================================================================================================
 \section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})}
 \label{sec:DYN_ldf}
-%------------------------------------------nam_dynldf----------------------------------------------------
 
 \begin{listing}
@@ -1102,5 +1066,4 @@
   \label{lst:namdyn_ldf}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{dyn_ldf}{dyn\_ldf} namelist variables.
@@ -1130,6 +1093,516 @@
 }
 
-%           Vertical diffusion term
-% External Forcing
-% Wetting and drying
-% Time evolution term
+%% =================================================================================================
+\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:DYN_ldf_lap}
+  \left\{
+    \begin{aligned}
+      D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}
+          \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[
+        {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\
+      D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}
+          \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[
+        {A_f^{lm} \;e_{3f} \zeta } \right]
+    \end{aligned}
+  \right.
+\end{equation}
+
+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.
+
+%% =================================================================================================
+\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[=.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[=.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.
+It must be emphasized that this formulation ignores constraints on the stress tensor such as symmetry.
+The resulting discrete representation is:
+\begin{equation}
+  \label{eq:DYN_ldf_iso}
+  \begin{split}
+    D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} }	\\
+    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_T^{lm}  \left(
+              {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta_{i}[u]
+                -e_{2t} \; r_{1t} \,\overline{\overline {\delta_{k+1/2}[u]}}^{\,i,\,k}}
+            \right)} \right]} 	\right. \\
+    & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f}
+            }\,\delta_{j+1/2} [u] - e_{1f}\, r_{2f}
+            \,\overline{\overline {\delta_{k+1/2} [u]}} ^{\,j+1/2,\,k}}
+        \right)} \right] \\
+    &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline
+              {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} }
+        \right.} \right. \\
+    &  \ \qquad \qquad \qquad \quad\
+    - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} \\
+    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
+                +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2}
+                \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} \\ \\
+    D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\
+    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_f^{lm}  \left(
+              {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta_{i+1/2}[v]
+                -e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}}
+            \right)} \right]} 	\right. \\
+    & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t}
+            }\,\delta_{j} [v] - e_{1t}\, r_{2t}
+            \,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}}
+        \right)} \right] \\
+    & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline
+              {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\
+    &  \ \qquad \qquad \qquad \quad\
+    - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\
+    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
+                +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2}
+                \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\}
+  \end{split}
+\end{equation}
+where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and
+the surface of computation ($z$- or $s$-surfaces).
+The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}).
+
+%% =================================================================================================
+\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: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}
+%%%
+
+%% =================================================================================================
+\section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})}
+\label{sec:DYN_zdf}
+
+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[=.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:MB_zdf} take the following semi-discrete space form:
+\[
+  % \label{eq:DYN_zdf}
+  \left\{
+    \begin{aligned}
+      D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta_k \left[ \frac{A_{uw}^{vm} }{e_{3uw} }
+        \ \delta_{k+1/2} [\,u\,]         \right]     \\
+      \\
+      D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta_k \left[ \frac{A_{vw}^{vm} }{e_{3vw} }
+        \ \delta_{k+1/2} [\,v\,]         \right]
+    \end{aligned}
+  \right.
+\]
+where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and diffusivity coefficients.
+The way these coefficients are evaluated depends on the vertical physics used (see \autoref{chap:ZDF}).
+
+The surface boundary condition on momentum is the stress exerted by the wind.
+At the surface, the momentum fluxes are prescribed as the boundary condition on
+the vertical turbulent momentum fluxes,
+\begin{equation}
+  \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 }
+\end{equation}
+where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in
+the (\textbf{i},\textbf{j}) coordinate system.
+The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in
+the vertical over the mixed layer depth.
+If the vertical mixing coefficient is small (when no mixed layer scheme is used)
+the surface stress enters only the top model level, as a body force.
+The surface wind stress is calculated in the surface module routines (SBC, see \autoref{chap:SBC}).
+
+The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation
+(see \autoref{sec:ZDF_drg})
+
+%% =================================================================================================
+\section{External forcings}
+\label{sec:DYN_forcing}
+
+Besides the surface and bottom stresses (see the above section)
+which are introduced as boundary conditions on the vertical mixing,
+three other forcings may enter the dynamical equations by affecting the surface pressure gradient.
+
+(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[=.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[=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
+ }
+
+%% =================================================================================================
+\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).
+The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO
+MS
+which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative
+limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$
+and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated
+by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$.
+
+\begin{listing}
+  \nlst{namwad}
+  \caption{\forcode{&namwad}}
+  \label{lst:namwad}
+\end{listing}
+
+The following terminology is used. The depth of the topography (positive downwards)
+at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the \NEMO\ code.
+The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign
+conventions used, the water depth, $h$, is the height of the free surface plus the depth of the
+topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$).
+
+Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be
+covered by water. They require the topography specified with a model
+configuration to have negative depths at points where the land is higher than the
+topography's reference sea-level. The vertical grid in \NEMO\ is normally computed relative to an
+initial state with zero sea surface height elevation.
+The user can choose to compute the vertical grid and heights in the model relative to
+a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths
+(i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays).
+
+Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''.
+$\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies
+with very steep slopes require larger values for normal choices of time-step. Surface fluxes
+are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers.
+The fluxes are tappered down using a $\mathrm{tanh}$ weighting function
+to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow.
+The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$.
+The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$.
+
+Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration
+and in ``realistic'' configurations covering parts of the north-west European shelf.
+All these configurations have used pure sigma coordinates. It is expected that
+the wetting and drying code will work in domains with more general s-coordinates provided
+the coordinates are pure sigma in the region where wetting and drying actually occurs.
+
+The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter.
+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})}
+\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}{rn\_wdmin1}).
+
+All the changes associated with this option are made to the barotropic solver for the non-linear
+free surface code within dynspg\_ts.
+On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells
+and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells
+whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells
+to be at least the minimum depth and hence is able to conserve mass / volume.
+
+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[=.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
+also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have
+chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells.
+The rationale for setting the velocity to zero is that it is the momentum equations that are being solved
+and the total momentum of the upstream cell (treating it as a finite volume) should be considered
+to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its
+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
+or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer
+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[=.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})}
+\label{subsec:DYN_wd_iterative_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''
+or may become dry within the next time-step using an iterative method.
+
+The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows:
+
+The continuity equation for the total water depth in a column
+\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{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}
+
+In the above $h$ is the depth of the water in the column at point $(i,j)$,
+$\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and
+$\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is
+the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t =
+t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and
+$\mathrm{zzflx}$ is the sum of the fluxes through all the faces.
+
+The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell
+(zzflxp) and fluxes that are into the cell (zzflxn).  Clearly
+
+\begin{equation}
+  \label{eq:DYN_wd_zzflx_p_n_1}
+  \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .
+\end{equation}
+
+The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until
+none of the cells will ``dry out''. To be precise the fluxes are limited until none of the
+cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$.
+
+Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and
+$\mathrm{flxv}^{(m)}$.  Then the adjustment is achieved by seeking a set of coefficients,
+$\mathrm{zcoef}_{i,j}^{(m)}$ such that:
+
+\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}
+
+where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around
+cells that would otherwise dry.
+
+The iteration is initialised by setting
+
+\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 (\autoref{eq:DYN_wd_continuity_2}) this
+condition is
+
+\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{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 (\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{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}
+
+Only the outward flux components are altered but, of course, outward fluxes from one cell
+are inward fluxes to adjacent cells and the balance in these cells may need subsequent
+adjustment; hence the iterative nature of this scheme.  Note, for example, that the flux
+across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration
+if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then
+the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by
+the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across
+the faces of the cells do not ``compete'' (they do not over-write each other) and one
+would expect the scheme to converge relatively quickly. The scheme is flux based so
+conserves mass. It also conserves constant tracers for the same reason that the
+directional limiter does.
+
+
+%      Surface pressure gradients
+%% =================================================================================================
+\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
+topography is sloping at these points the sea-surface will have a similar slope and there
+will hence be very large horizontal pressure gradients at these points. The WAD modifies
+the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such
+points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie
+between $0$ and $1$.
+
+We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated
+at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the
+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 \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
+both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and
+the minimum height of the sea surface at the two points is greater than the maximum height
+of the topography at the two points:
+
+\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}
+\end{equation}
+
+The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height
+of the sea surface at the two points is greater than the maximum height of the topography
+at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$
+
+\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}
+
+If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$
+point is unmodified. If both logicals are false zhpi is set to zero.
+
+If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure
+gradient is multiplied through by zcpx which is the absolute value of the difference in
+the water depths at the two points divided by the difference in the surface heights at the
+two points. Thus the sign of the sea surface height gradient is retained but the magnitude
+of the pressure force is determined by the difference in water depths rather than the
+difference in surface height between the two points. Note that dividing by the difference
+between the sea surface heights can be problematic if the heights approach parity. An
+additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such
+conditions.
+
+%% =================================================================================================
+\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
+bottom drag is clearly very important. In order to promote stability
+it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach.
+
+Suitable specifcation of the surface heat flux in wetting and drying domains in forced and
+coupled simulations needs further consideration. In order to prevent freezing or boiling
+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{subsec:DYN_WAD_test_cases}
+
+See the WAD tests MY\_DOC documention for details of the WAD test cases.
+
+
+
+%% =================================================================================================
+\section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})}
+\label{sec:DYN_nxt}
+
+
+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: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:SCOORD_momentum})
+
+$\bullet$ vector invariant form or linear free surface
+(\np[=.true.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} not defined):
+\[
+  % \label{eq:DYN_nxt_vec}
+  \left\{
+    \begin{aligned}
+      &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt  \ \text{RHS}_u^t  	\\
+      &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right]
+    \end{aligned}
+  \right.
+\]
+
+$\bullet$ flux form and nonlinear free surface
+(\np[=.false.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} defined):
+\[
+  % \label{eq:DYN_nxt_flux}
+  \left\{
+    \begin{aligned}
+      &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t  	\\
+      &\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t
+      +\gamma \,\left[ {\left(e_{3u}\,u\right)_f^{t-\rdt} -2\left(e_{3u}\,u\right)^t+\left(e_{3u}\,u\right)^{t+\rdt}} \right]
+    \end{aligned}
+  \right.
+\]
+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}{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.
+
+Note that with the filtered free surface,
+the update of the \textit{after} velocities is done in the \mdl{dynsp\_flt} module,
+and only array swapping and Asselin filtering is done in \mdl{dynnxt}.
+
+\onlyinsubfile{\input{../../global/epilogue}}
+
+\end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex	(revision 11597)
@@ -9,7 +9,7 @@
 %gm% add here introduction to this chapter
 
+%% =================================================================================================
 \section[Boundary condition at the coast (\forcode{rn_shlat})]{Boundary condition at the coast (\protect\np{rn_shlat}{rn\_shlat})}
 \label{sec:LBC_coast}
-%--------------------------------------------namlbc-------------------------------------------------------
 
 \begin{listing}
@@ -18,5 +18,4 @@
   \label{lst:namlbc}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 %The lateral ocean boundary conditions contiguous to coastlines are Neumann conditions for heat and salt
@@ -142,4 +141,5 @@
 it is only applied next to the coast where the minimum water depth can be quite shallow.
 
+%% =================================================================================================
 \section[Model domain boundary condition (\forcode{jperio})]{Model domain boundary condition (\protect\jp{jperio})}
 \label{sec:LBC_jperio}
@@ -150,4 +150,5 @@
 The north-fold boundary condition is associated with the 3-pole ORCA mesh.
 
+%% =================================================================================================
 \subsection[Closed, cyclic (\forcode{=0,1,2,7})]{Closed, cyclic (\protect\jp{jperio}\forcode{=0,1,2,7})}
 \label{subsec:LBC_jperio012}
@@ -191,4 +192,5 @@
 \end{figure}
 
+%% =================================================================================================
 \subsection[North-fold (\forcode{=3,6})]{North-fold (\protect\jp{jperio}\forcode{=3,6})}
 \label{subsec:LBC_north_fold}
@@ -210,8 +212,8 @@
 \end{figure}
 
+%% =================================================================================================
 \section[Exchange with neighbouring processors (\textit{lbclnk.F90}, \textit{lib\_mpp.F90})]{Exchange with neighbouring processors (\protect\mdl{lbclnk}, \protect\mdl{lib\_mpp})}
 \label{sec:LBC_mpp}
 
-%-----------------------------------------nammpp--------------------------------------------
 
 \begin{listing}
@@ -220,5 +222,4 @@
   \label{lst:nammpp}
 \end{listing}
-%-----------------------------------------------------------------------------------------------
 
 For massively parallel processing (mpp), a domain decomposition method is used.
@@ -321,8 +322,8 @@
 \end{figure}
 
+%% =================================================================================================
 \section{Unstructured open boundary conditions (BDY)}
 \label{sec:LBC_bdy}
 
-%-----------------------------------------nambdy--------------------------------------------
 
 \begin{listing}
@@ -331,6 +332,4 @@
   \label{lst:nambdy}
 \end{listing}
-%-----------------------------------------------------------------------------------------------
-%-----------------------------------------nambdy_dta--------------------------------------------
 
 \begin{listing}
@@ -339,5 +338,4 @@
   \label{lst:nambdy_dta}
 \end{listing}
-%-----------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{bdy}{bdy} and \nam{bdy_dta}{bdy\_dta} namelist variables.
@@ -352,5 +350,5 @@
 See the section on the Input Boundary Data Files for details.
 
-%----------------------------------------------
+%% =================================================================================================
 \subsection{Namelists}
 \label{subsec:LBC_bdy_namelist}
@@ -419,5 +417,5 @@
 FRS conditions are applied on temperature and salinity and climatological data is read from initial condition files.
 
-%----------------------------------------------
+%% =================================================================================================
 \subsection{Flow relaxation scheme}
 \label{subsec:LBC_bdy_FRS_scheme}
@@ -458,5 +456,5 @@
 This is typically set to a value between 8 and 10.
 
-%----------------------------------------------
+%% =================================================================================================
 \subsection{Flather radiation scheme}
 \label{subsec:LBC_bdy_flather_scheme}
@@ -480,5 +478,5 @@
 $U$ and $U_{e}$ are defined on the $U$ or $V$ points with $nbr=1$, \ie\ between the two $T$ grid points.
 
-%----------------------------------------------
+%% =================================================================================================
 \subsection{Orlanski radiation scheme}
 \label{subsec:LBC_bdy_orlanski_scheme}
@@ -528,5 +526,5 @@
 (\autoref{eq:LBC_wave_continuous}) - (\autoref{eq:LBC_tau_in}) correspond to equations (13) - (15) and (2) - (3) in MMS.\\
 
-%----------------------------------------------
+%% =================================================================================================
 \subsection{Relaxation at the boundary}
 \label{subsec:LBC_bdy_relaxation}
@@ -544,5 +542,5 @@
 The same scaling is applied in the Orlanski damping.
 
-%----------------------------------------------
+%% =================================================================================================
 \subsection{Boundary geometry}
 \label{subsec:LBC_bdy_geometry}
@@ -590,5 +588,5 @@
 \end{figure}
 
-%----------------------------------------------
+%% =================================================================================================
 \subsection{Input boundary data files}
 \label{subsec:LBC_bdy_data}
@@ -626,5 +624,5 @@
 \end{figure}
 
-%----------------------------------------------
+%% =================================================================================================
 \subsection{Volume correction}
 \label{subsec:LBC_bdy_vol_corr}
@@ -643,9 +641,8 @@
 applied to all boundaries at once.
 
-%----------------------------------------------
+%% =================================================================================================
 \subsection{Tidal harmonic forcing}
 \label{subsec:LBC_bdy_tides}
 
-%-----------------------------------------nambdy_tide--------------------------------------------
 
 \begin{listing}
@@ -654,5 +651,4 @@
   \label{lst:nambdy_tide}
 \end{listing}
-%-----------------------------------------------------------------------------------------------
 
 Tidal forcing at open boundaries requires the activation of surface
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex	(revision 11597)
@@ -22,12 +22,12 @@
 is described in \autoref{apdx:TRIADS}
 
-%-----------------------------------namtra_ldf - namdyn_ldf--------------------------------------------
-
-%--------------------------------------------------------------------------------------------------------------
-
+
+
+%% =================================================================================================
 \section[Lateral mixing operators]{Lateral mixing operators}
 \label{sec:LDF_op}
 We remind here the different lateral mixing operators that can be used. Further details can be found in \autoref{subsec:TRA_ldf_op} and  \autoref{sec:DYN_ldf}.
 
+%% =================================================================================================
 \subsection[No lateral mixing (\forcode{ln_traldf_OFF} \& \forcode{ln_dynldf_OFF})]{No lateral mixing (\protect\np{ln_traldf_OFF}{ln\_traldf\_OFF} \& \protect\np{ln_dynldf_OFF}{ln\_dynldf\_OFF})}
 
@@ -37,4 +37,5 @@
 see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes.
 
+%% =================================================================================================
 \subsection[Laplacian mixing (\forcode{ln_traldf_lap} \& \forcode{ln_dynldf_lap})]{Laplacian mixing (\protect\np{ln_traldf_lap}{ln\_traldf\_lap} \& \protect\np{ln_dynldf_lap}{ln\_dynldf\_lap})}
 Setting \protect\np[=.true.]{ln_traldf_lap}{ln\_traldf\_lap} and/or \protect\np[=.true.]{ln_dynldf_lap}{ln\_dynldf\_lap} enables
@@ -42,4 +43,5 @@
 Laplacian and Bilaplacian operators for the same variable.
 
+%% =================================================================================================
 \subsection[Bilaplacian mixing (\forcode{ln_traldf_blp} \& \forcode{ln_dynldf_blp})]{Bilaplacian mixing (\protect\np{ln_traldf_blp}{ln\_traldf\_blp} \& \protect\np{ln_dynldf_blp}{ln\_dynldf\_blp})}
 Setting \protect\np[=.true.]{ln_traldf_blp}{ln\_traldf\_blp} and/or \protect\np[=.true.]{ln_dynldf_blp}{ln\_dynldf\_blp} enables
@@ -47,4 +49,5 @@
 We stress again that from \NEMO\ 4, the simultaneous use Laplacian and Bilaplacian operators is not allowed.
 
+%% =================================================================================================
 \section[Direction of lateral mixing (\textit{ldfslp.F90})]{Direction of lateral mixing (\protect\mdl{ldfslp})}
 \label{sec:LDF_slp}
@@ -69,4 +72,5 @@
 %gm% add here afigure of the slope in i-direction
 
+%% =================================================================================================
 \subsection{Slopes for tracer geopotential mixing in the $s$-coordinate}
 
@@ -99,4 +103,5 @@
 and either \np[=.true.]{ln_traldf_hor}{ln\_traldf\_hor} or \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}.
 
+%% =================================================================================================
 \subsection{Slopes for tracer iso-neutral mixing}
 \label{subsec:LDF_slp_iso}
@@ -273,4 +278,5 @@
 \colorbox{yellow}{add here a discussion about the flattening of the slopes, vs tapering the coefficient.}
 
+%% =================================================================================================
 \subsection{Slopes for momentum iso-neutral mixing}
 
@@ -299,4 +305,5 @@
 (see \autoref{sec:LBC_coast}).
 
+%% =================================================================================================
 \section[Lateral mixing coefficient (\forcode{nn_aht_ijk_t} \& \forcode{nn_ahm_ijk_t})]{Lateral mixing coefficient (\protect\np{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 \label{sec:LDF_coef}
@@ -305,4 +312,5 @@
 The way the mixing coefficients are set in the reference version can be described as follows:
 
+%% =================================================================================================
 \subsection[Mixing coefficients read from file (\forcode{=-20, -30})]{ Mixing coefficients read from file (\protect\np[=-20, -30]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=-20, -30]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
@@ -313,5 +321,4 @@
 The provided fields can either be 2d (\np[=-20]{nn_aht_ijk_t}{nn\_aht\_ijk\_t}, \np[=-20]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}) or 3d (\np[=-30]{nn_aht_ijk_t}{nn\_aht\_ijk\_t},  \np[=-30]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}). They must be given at U, V points for tracers and T, F points for momentum (see \autoref{tab:LDF_files}).
 
-%-------------------------------------------------TABLE---------------------------------------------------
 \begin{table}[htb]
   \centering
@@ -327,6 +334,6 @@
   \label{tab:LDF_files}
 \end{table}
-%--------------------------------------------------------------------------------------------------------------
-
+
+%% =================================================================================================
 \subsection[Constant mixing coefficients (\forcode{=0})]{ Constant mixing coefficients (\protect\np[=0]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=0]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
@@ -345,4 +352,5 @@
  $U_{scl}$ and $L_{scl}$ are given by the namelist parameters \np{rn_Ud}{rn\_Ud}, \np{rn_Uv}{rn\_Uv}, \np{rn_Ld}{rn\_Ld} and \np{rn_Lv}{rn\_Lv}.
 
+%% =================================================================================================
 \subsection[Vertically varying mixing coefficients (\forcode{=10})]{Vertically varying mixing coefficients (\protect\np[=10]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=10]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
@@ -352,4 +360,5 @@
 This profile is hard coded in module \mdl{ldfc1d\_c2d}, but can be easily modified by users.
 
+%% =================================================================================================
 \subsection[Mesh size dependent mixing coefficients (\forcode{=20})]{Mesh size dependent mixing coefficients (\protect\np[=20]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=20]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
@@ -377,4 +386,5 @@
 \colorbox{yellow}{CASE \np{nn_aht_ijk_t}{nn\_aht\_ijk\_t} = 21 to be added}
 
+%% =================================================================================================
 \subsection[Mesh size and depth dependent mixing coefficients (\forcode{=30})]{Mesh size and depth dependent mixing coefficients (\protect\np[=30]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=30]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
@@ -383,4 +393,5 @@
 the magnitude of the coefficient.
 
+%% =================================================================================================
 \subsection[Velocity dependent mixing coefficients (\forcode{=31})]{Flow dependent mixing coefficients (\protect\np[=31]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=31]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 In that case, the eddy coefficient is proportional to the local velocity magnitude so that the Reynolds number $Re =  \lvert U \rvert  e / A_l$ is constant (and here hardcoded to $12$):
@@ -397,4 +408,5 @@
 \end{equation}
 
+%% =================================================================================================
 \subsection[Deformation rate dependent viscosities (\forcode{nn_ahm_ijk_t=32})]{Deformation rate dependent viscosities (\protect\np[=32]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
@@ -433,4 +445,5 @@
 where $C_{min}$ and $C_{max}$ are adimensional namelist parameters given by \np{rn_minfac}{rn\_minfac} and \np{rn_maxfac}{rn\_maxfac} respectively.
 
+%% =================================================================================================
 \subsection{About space and time varying mixing coefficients}
 
@@ -446,9 +459,9 @@
 (\autoref{sec:INVARIANTS_dynldf_properties}).
 
+%% =================================================================================================
 \section[Eddy induced velocity (\forcode{ln_ldfeiv})]{Eddy induced velocity (\protect\np{ln_ldfeiv}{ln\_ldfeiv})}
 
 \label{sec:LDF_eiv}
 
-%--------------------------------------------namtra_eiv---------------------------------------------------
 
 \begin{listing}
@@ -458,5 +471,4 @@
 \end{listing}
 
-%--------------------------------------------------------------------------------------------------------------
 
 %%gm  from Triad appendix  : to be incorporated....
@@ -513,8 +525,8 @@
 In case of setting \np[=.true.]{ln_traldf_triad}{ln\_traldf\_triad}, a skew form of the eddy induced advective fluxes is used, which is described in \autoref{apdx:TRIADS}.
 
+%% =================================================================================================
 \section[Mixed layer eddies (\forcode{ln_mle})]{Mixed layer eddies (\protect\np{ln_mle}{ln\_mle})}
 \label{sec:LDF_mle}
 
-%--------------------------------------------namtra_eiv---------------------------------------------------
 
 \begin{listing}
@@ -524,5 +536,4 @@
 \end{listing}
 
-%--------------------------------------------------------------------------------------------------------------
 
 If  \np[=.true.]{ln_mle}{ln\_mle} in \nam{tra_mle}{tra\_mle} namelist, a parameterization of the mixing due to unresolved mixed layer instabilities is activated (\citet{foxkemper.ferrari_JPO08}). Additional transport is computed in \rou{ldf\_mle\_trp} and added to the eulerian transport in \rou{tra\_adv} as done for eddy induced advection.
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex	(revision 11597)
@@ -9,4 +9,5 @@
 \vfill
 \begin{figure}[b]
+%% =================================================================================================
 \subsubsection*{Changes record}
 \begin{tabular}{l||l|m{0.65\linewidth}}
@@ -55,4 +56,5 @@
 In \autoref{sec:OBS_obsutils} we describe some utilities to help work with the files produced by the OBS code.
 
+%% =================================================================================================
 \section{Running the observation operator code example}
 \label{sec:OBS_example}
@@ -103,4 +105,5 @@
 \autoref{sec:OBS_obsutils}.
 
+%% =================================================================================================
 \section{Technical details (feedback type observation file headers)}
 \label{sec:OBS_details}
@@ -109,5 +112,4 @@
 the observation files that may be used with the observation operator.
 
-%------------------------------------------namobs--------------------------------------------------------
 
 \begin{listing}
@@ -116,5 +118,4 @@
   \label{lst:namobs}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 The observation operator code uses the feedback observation file format for all data types.
@@ -123,4 +124,5 @@
 sea surface temperature are in the following subsections.
 
+%% =================================================================================================
 \subsection{Profile feedback file}
 
@@ -279,4 +281,5 @@
 \end{clines}
 
+%% =================================================================================================
 \subsection{Sea level anomaly feedback file}
 
@@ -425,4 +428,5 @@
 \end{clines}
 
+%% =================================================================================================
 \subsection{Sea surface temperature feedback file}
 
@@ -542,7 +546,9 @@
 \end{clines}
 
+%% =================================================================================================
 \section{Theoretical details}
 \label{sec:OBS_theory}
 
+%% =================================================================================================
 \subsection{Horizontal interpolation and averaging methods}
 
@@ -579,4 +585,5 @@
 Below is some more detail on the various options for interpolation and averaging available in \NEMO.
 
+%% =================================================================================================
 \subsubsection{Horizontal interpolation}
 
@@ -667,4 +674,5 @@
 \end{enumerate}
 
+%% =================================================================================================
 \subsubsection{Horizontal averaging}
 
@@ -708,4 +716,5 @@
 \end{figure}
 
+%% =================================================================================================
 \subsection{Grid search}
 
@@ -758,4 +767,5 @@
 the $i$ and $j$ ranges of this point searched to determine the precise four points surrounding the observation.
 
+%% =================================================================================================
 \subsection{Parallel aspects of horizontal interpolation}
 \label{subsec:OBS_parallel}
@@ -772,4 +782,5 @@
 and 2) round-robin.
 
+%% =================================================================================================
 \subsubsection{Geographical distribution of observations among processors}
 
@@ -798,4 +809,5 @@
 this could lead to load imbalance.
 
+%% =================================================================================================
 \subsubsection{Round-robin distribution of observations among processors}
 
@@ -819,4 +831,5 @@
 a subroutine has been developed that retrieves any grid points in the global space.
 
+%% =================================================================================================
 \subsection{Vertical interpolation operator}
 
@@ -830,7 +843,9 @@
 %\usepackage{framed}
 
+%% =================================================================================================
 \section{Standalone observation operator}
 \label{sec:OBS_sao}
 
+%% =================================================================================================
 \subsection{Concept}
 
@@ -849,10 +864,10 @@
 By forecast, we mean any method which produces an estimate of physical reality which is not an observed value.
 
-%--------------------------------------------------------------------------------------------------------
 % sao.exe
-%--------------------------------------------------------------------------------------------------------
-
+
+%% =================================================================================================
 \subsection{Using the standalone observation operator}
 
+%% =================================================================================================
 \subsubsection{Building}
 
@@ -862,7 +877,6 @@
 Note this a similar approach to that taken by the standalone surface scheme \emph{SAS\_SRC} and the offline TOP model \emph{OFF\_SRC}.
 
-%--------------------------------------------------------------------------------------------------------
 % Running
-%--------------------------------------------------------------------------------------------------------
+%% =================================================================================================
 \subsubsection{Running}
 
@@ -870,11 +884,11 @@
 a full \NEMO\ namelist and then to run the executable as if it were nemo.exe.
 
-%--------------------------------------------------------------------------------------------------------
 % Configuration section
-%--------------------------------------------------------------------------------------------------------
+%% =================================================================================================
 \subsection{Configuring the standalone observation operator}
 The observation files and settings understood by \nam{obs}{obs} have been outlined in the online observation operator section.
 In addition is a further namelist \nam{sao}{sao} which used to set the input model fields for the SAO
 
+%% =================================================================================================
 \subsubsection{Single field}
 
@@ -901,4 +915,5 @@
 \end{forlines}
 
+%% =================================================================================================
 \subsubsection{Multiple fields per run}
 
@@ -936,4 +951,5 @@
 This approach is referred to as \emph{Class 4} since it is the fourth metric defined by the GODAE intercomparison project. This requires multiple runs of the SAO and running an additional utility (not currently in the \NEMO\ repository) to combine the feedback files into one class 4 file.
 
+%% =================================================================================================
 \section{Observation utilities}
 \label{sec:OBS_obsutils}
@@ -948,4 +964,5 @@
 OBSTOOLS and dataplot are described in more detail below.
 
+%% =================================================================================================
 \subsection{Obstools}
 
@@ -953,4 +970,5 @@
 This are helpful in handling observation files and the feedback file output from the observation operator. A brief description of some of the utilities follows
 
+%% =================================================================================================
 \subsubsection{corio2fb}
 
@@ -962,4 +980,5 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsubsection{enact2fb}
 
@@ -971,4 +990,5 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsubsection{fbcomb}
 
@@ -981,4 +1001,5 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsubsection{fbmatchup}
 
@@ -990,4 +1011,5 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsubsection{fbprint}
 
@@ -1018,4 +1040,5 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsubsection{fbsel}
 
@@ -1027,4 +1050,5 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsubsection{fbstat}
 
@@ -1036,4 +1060,5 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsubsection{fbthin}
 
@@ -1046,4 +1071,5 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsubsection{sla2fb}
 
@@ -1058,4 +1084,5 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsubsection{vel2fb}
 
@@ -1067,8 +1094,10 @@
 \end{cmds}
 
+%% =================================================================================================
 \subsection{Building the obstools}
 
 To build the obstools use in the tools directory use ./maketools -n OBSTOOLS -m [ARCH].
 
+%% =================================================================================================
 \subsection{Dataplot}
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 11597)
@@ -8,5 +8,4 @@
 \chaptertoc
 
-%---------------------------------------namsbc--------------------------------------------------
 
 \begin{listing}
@@ -15,5 +14,4 @@
   \label{lst:namsbc}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 The ocean needs seven fields as surface boundary condition:
@@ -87,4 +85,5 @@
 which provides additional sources of fresh water.
 
+%% =================================================================================================
 \section{Surface boundary condition for the ocean}
 \label{sec:SBC_ocean}
@@ -147,5 +146,4 @@
 these averaged fields are used to compute the surface fluxes at the frequency of \np{nn_fsbc}{nn\_fsbc} time-steps.
 
-%-------------------------------------------------TABLE---------------------------------------------------
 \begin{table}[tb]
   \centering
@@ -167,8 +165,8 @@
   \label{tab:SBC_ssm}
 \end{table}
-%--------------------------------------------------------------------------------------------------------------
 
 %\colorbox{yellow}{Penser a} mettre dans le restant l'info nn\_fsbc ET nn\_fsbc*rdt de sorte de reinitialiser la moyenne si on change la frequence ou le pdt
 
+%% =================================================================================================
 \section{Input data generic interface}
 \label{sec:SBC_input}
@@ -203,4 +201,5 @@
 By default, the data are assumed to be in the same directory as the executable, so that cn\_dir='./'.
 
+%% =================================================================================================
 \subsection[Input data specification (\textit{fldread.F90})]{Input data specification (\protect\mdl{fldread})}
 \label{subsec:SBC_fldread}
@@ -220,5 +219,4 @@
   whether it is a climatological file or not, and to the open/close frequency (see below for definition).
 
-%--------------------------------------------------TABLE--------------------------------------------------
   \begin{table}[htbp]
     \centering
@@ -247,5 +245,4 @@
     \label{tab:SBC_fldread}
   \end{table}
-%--------------------------------------------------------------------------------------------------------------
 
 \item [Record frequency]:
@@ -325,4 +322,5 @@
 a useful feature for user considering that it is too heavy to manipulate the complete file for year Y-1.
 
+%% =================================================================================================
 \subsection{Interpolation on-the-fly}
 \label{subsec:SBC_iof}
@@ -347,4 +345,5 @@
 Note that nn\_lsm=0 forces the code to not apply the procedure, even if a land/sea mask file is supplied.
 
+%% =================================================================================================
 \subsubsection{Bilinear interpolation}
 \label{subsec:SBC_iof_bilinear}
@@ -368,4 +367,5 @@
 and wgt(1) corresponds to variable "wgt01" for example.
 
+%% =================================================================================================
 \subsubsection{Bicubic interpolation}
 \label{subsec:SBC_iof_bicubic}
@@ -386,4 +386,5 @@
 the spatial dependency has been included into the weights.
 
+%% =================================================================================================
 \subsubsection{Implementation}
 \label{subsec:SBC_iof_imp}
@@ -421,4 +422,5 @@
 or is a copy of one from the first few columns on the opposite side of the grid (cyclical case).
 
+%% =================================================================================================
 \subsubsection{Limitations}
 \label{subsec:SBC_iof_lim}
@@ -435,4 +437,5 @@
 \end{enumerate}
 
+%% =================================================================================================
 \subsubsection{Utilities}
 \label{subsec:SBC_iof_util}
@@ -442,8 +445,8 @@
 (see the directory NEMOGCM/TOOLS/WEIGHTS).
 
+%% =================================================================================================
 \subsection{Standalone surface boundary condition scheme (SAS)}
 \label{subsec:SBC_SAS}
 
-%---------------------------------------namsbc_sas--------------------------------------------------
 
 \begin{listing}
@@ -452,5 +455,4 @@
   \label{lst:namsbc_sas}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 In some circumstances, it may be useful to avoid calculating the 3D temperature,
@@ -513,7 +515,7 @@
  (\np[=.true.]{ln_flx}{ln\_flx}) and to provide 3D oceanic velocities instead of 2D ones (\np{ln_flx}{ln\_flx}\forcode{=.true.}). In that last case, only the 1st level will be read in.
 
+%% =================================================================================================
 \section[Flux formulation (\textit{sbcflx.F90})]{Flux formulation (\protect\mdl{sbcflx})}
 \label{sec:SBC_flx}
-%------------------------------------------namsbc_flx----------------------------------------------------
 
 \begin{listing}
@@ -522,5 +524,4 @@
   \label{lst:namsbc_flx}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 In the flux formulation (\np[=.true.]{ln_flx}{ln\_flx}),
@@ -534,7 +535,7 @@
 See \autoref{subsec:SBC_ssr} for its specification.
 
+%% =================================================================================================
 \section[Bulk formulation (\textit{sbcblk.F90})]{Bulk formulation (\protect\mdl{sbcblk})}
 \label{sec:SBC_blk}
-%---------------------------------------namsbc_blk--------------------------------------------------
 
 \begin{listing}
@@ -543,5 +544,4 @@
   \label{lst:namsbc_blk}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 In the bulk formulation, the surface boundary condition fields are computed with bulk formulae using atmospheric fields
@@ -562,5 +562,4 @@
 The required 9 input fields are:
 
-%--------------------------------------------------TABLE--------------------------------------------------
 \begin{table}[htbp]
   \centering
@@ -590,5 +589,4 @@
   \label{tab:SBC_BULK}
 \end{table}
-%--------------------------------------------------------------------------------------------------------------
 
 Note that the air velocity is provided at a tracer ocean point, not at a velocity ocean point ($u$- and $v$-points).
@@ -615,4 +613,5 @@
 the namsbc\_blk namelist (see \autoref{subsec:SBC_fldread}).
 
+%% =================================================================================================
 \subsection[Ocean-Atmosphere Bulk formulae (\textit{sbcblk\_algo\_coare.F90, sbcblk\_algo\_coare3p5.F90, sbcblk\_algo\_ecmwf.F90, sbcblk\_algo\_ncar.F90})]{Ocean-Atmosphere Bulk formulae (\mdl{sbcblk\_algo\_coare}, \mdl{sbcblk\_algo\_coare3p5}, \mdl{sbcblk\_algo\_ecmwf}, \mdl{sbcblk\_algo\_ncar})}
 \label{subsec:SBC_blk_ocean}
@@ -640,4 +639,5 @@
 \end{itemize}
 
+%% =================================================================================================
 \subsection{Ice-Atmosphere Bulk formulae}
 \label{subsec:SBC_blk_ice}
@@ -665,7 +665,7 @@
 \end{itemize}
 
+%% =================================================================================================
 \section[Coupled formulation (\textit{sbccpl.F90})]{Coupled formulation (\protect\mdl{sbccpl})}
 \label{sec:SBC_cpl}
-%------------------------------------------namsbc_cpl----------------------------------------------------
 
 \begin{listing}
@@ -674,5 +674,4 @@
   \label{lst:namsbc_cpl}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 In the coupled formulation of the surface boundary condition,
@@ -702,7 +701,7 @@
 In cases where this is definitely not possible, the model should abort with an error message.
 
+%% =================================================================================================
 \section[Atmospheric pressure (\textit{sbcapr.F90})]{Atmospheric pressure (\protect\mdl{sbcapr})}
 \label{sec:SBC_apr}
-%------------------------------------------namsbc_apr----------------------------------------------------
 
 \begin{listing}
@@ -711,5 +710,4 @@
   \label{lst:namsbc_apr}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 The optional atmospheric pressure can be used to force ocean and ice dynamics
@@ -739,8 +737,8 @@
 \np{ln_apr_obc}{ln\_apr\_obc}  might be set to true.
 
+%% =================================================================================================
 \section[Surface tides (\textit{sbctide.F90})]{Surface tides (\protect\mdl{sbctide})}
 \label{sec:SBC_tide}
 
-%------------------------------------------nam_tide---------------------------------------
 
 \begin{listing}
@@ -749,5 +747,4 @@
   \label{lst:nam_tide}
 \end{listing}
-%-----------------------------------------------------------------------------------------
 
 The tidal forcing, generated by the gravity forces of the Earth-Moon and Earth-Sun sytems,
@@ -791,7 +788,7 @@
 \forcode{.false.} removes the SAL contribution.
 
+%% =================================================================================================
 \section[River runoffs (\textit{sbcrnf.F90})]{River runoffs (\protect\mdl{sbcrnf})}
 \label{sec:SBC_rnf}
-%------------------------------------------namsbc_rnf----------------------------------------------------
 
 \begin{listing}
@@ -800,5 +797,4 @@
   \label{lst:namsbc_rnf}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 %River runoff generally enters the ocean at a nonzero depth rather than through the surface.
@@ -913,7 +909,7 @@
 %To do this we need to treat evaporation/precipitation fluxes and river runoff differently in the tra_sbc module.  We decided to separate them throughout the code, so that the variable emp represented solely evaporation minus precipitation fluxes, and a new 2d variable rnf was added which represents the volume flux of river runoff (in kg/m2s to remain consistent with emp).  This meant many uses of emp and emps needed to be changed, a list of all modules which use emp or emps and the changes made are below:
 
+%% =================================================================================================
 \section[Ice shelf melting (\textit{sbcisf.F90})]{Ice shelf melting (\protect\mdl{sbcisf})}
 \label{sec:SBC_isf}
-%------------------------------------------namsbc_isf----------------------------------------------------
 
 \begin{listing}
@@ -922,5 +918,4 @@
   \label{lst:namsbc_isf}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------
 
 The namelist variable in \nam{sbc}{sbc}, \np{nn_isf}{nn\_isf}, controls the ice shelf representation.
@@ -1029,7 +1024,7 @@
 \end{figure}
 
+%% =================================================================================================
 \section{Ice sheet coupling}
 \label{sec:SBC_iscpl}
-%------------------------------------------namsbc_iscpl----------------------------------------------------
 
 \begin{listing}
@@ -1038,5 +1033,4 @@
   \label{lst:namsbc_iscpl}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------
 
 Ice sheet/ocean coupling is done through file exchange at the restart step.
@@ -1093,7 +1087,7 @@
 The corrective increment is apply into the cell itself (if it is a wet cell), the neigbouring cells or the closest wet cell (if the cell is now dry).
 
+%% =================================================================================================
 \section{Handling of icebergs (ICB)}
 \label{sec:SBC_ICB_icebergs}
-%------------------------------------------namberg----------------------------------------------------
 
 \begin{listing}
@@ -1102,5 +1096,4 @@
   \label{lst:namberg}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 Icebergs are modelled as lagrangian particles in \NEMO\ \citep{marsh.ivchenko.ea_GMD15}.
@@ -1162,7 +1155,7 @@
 since its trajectory data may be spread across multiple files.
 
+%% =================================================================================================
 \section[Interactions with waves (\textit{sbcwave.F90}, \forcode{ln_wave})]{Interactions with waves (\protect\mdl{sbcwave}, \protect\np{ln_wave}{ln\_wave})}
 \label{sec:SBC_wave}
-%------------------------------------------namsbc_wave--------------------------------------------------------
 
 \begin{listing}
@@ -1171,5 +1164,4 @@
   \label{lst:namsbc_wave}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 Ocean waves represent the interface between the ocean and the atmosphere, so \NEMO\ is extended to incorporate
@@ -1196,4 +1188,5 @@
 
 % ----------------------------------------------------------------
+%% =================================================================================================
 \subsection[Neutral drag coefficient from wave model (\forcode{ln_cdgw})]{Neutral drag coefficient from wave model (\protect\np{ln_cdgw}{ln\_cdgw})}
 \label{subsec:SBC_wave_cdgw}
@@ -1208,4 +1201,5 @@
 % 3D Stokes Drift (ln_sdw, nn_sdrift)
 % ----------------------------------------------------------------
+%% =================================================================================================
 \subsection[3D Stokes Drift (\forcode{ln_sdw} \& \forcode{nn_sdrift})]{3D Stokes Drift (\protect\np{ln_sdw}{ln\_sdw} \& \np{nn_sdrift}{nn\_sdrift})}
 \label{subsec:SBC_wave_sdw}
@@ -1303,4 +1297,5 @@
 % Stokes-Coriolis term (ln_stcor)
 % ----------------------------------------------------------------
+%% =================================================================================================
 \subsection[Stokes-Coriolis term (\forcode{ln_stcor})]{Stokes-Coriolis term (\protect\np{ln_stcor}{ln\_stcor})}
 \label{subsec:SBC_wave_stcor}
@@ -1316,4 +1311,5 @@
 % Waves modified stress (ln_tauwoc, ln_tauw)
 % ----------------------------------------------------------------
+%% =================================================================================================
 \subsection[Wave modified stress (\forcode{ln_tauwoc} \& \forcode{ln_tauw})]{Wave modified sress (\protect\np{ln_tauwoc}{ln\_tauwoc} \& \np{ln_tauw}{ln\_tauw})}
 \label{subsec:SBC_wave_tauw}
@@ -1353,13 +1349,13 @@
 meridional stress components by setting \np[=.true.]{ln_tauw}{ln\_tauw}.
 
+%% =================================================================================================
 \section{Miscellaneous options}
 \label{sec:SBC_misc}
 
+%% =================================================================================================
 \subsection[Diurnal cycle (\textit{sbcdcy.F90})]{Diurnal cycle (\protect\mdl{sbcdcy})}
 \label{subsec:SBC_dcy}
-%------------------------------------------namsbc-------------------------------------------------------------
 %
 
-%-------------------------------------------------------------------------------------------------------------
 
 \begin{figure}[!t]
@@ -1411,4 +1407,5 @@
 an inconsistency between the scale of the vertical resolution and the forcing acting on that scale.
 
+%% =================================================================================================
 \subsection{Rotation of vector pairs onto the model grid directions}
 \label{subsec:SBC_rotation}
@@ -1427,7 +1424,7 @@
 The rot\_rep routine from the \mdl{geo2ocean} module is used to perform the rotation.
 
+%% =================================================================================================
 \subsection[Surface restoring to observed SST and/or SSS (\textit{sbcssr.F90})]{Surface restoring to observed SST and/or SSS (\protect\mdl{sbcssr})}
 \label{subsec:SBC_ssr}
-%------------------------------------------namsbc_ssr----------------------------------------------------
 
 \begin{listing}
@@ -1436,5 +1433,4 @@
   \label{lst:namsbc_ssr}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{sbc_ssr}{sbc\_ssr} namelist variables.
@@ -1471,4 +1467,5 @@
 reduce the uncertainties we have on the observed freshwater budget.
 
+%% =================================================================================================
 \subsection{Handling of ice-covered area  (\textit{sbcice\_...})}
 \label{subsec:SBC_ice-cover}
@@ -1508,4 +1505,5 @@
 %GS: ocean-ice (SI3) interface is not located in SBC directory anymore, so it should be included in SI3 doc
 
+%% =================================================================================================
 \subsection[Interface to CICE (\textit{sbcice\_cice.F90})]{Interface to CICE (\protect\mdl{sbcice\_cice})}
 \label{subsec:SBC_cice}
@@ -1538,4 +1536,5 @@
 there is no sea ice.
 
+%% =================================================================================================
 \subsection[Freshwater budget control (\textit{sbcfwb.F90})]{Freshwater budget control (\protect\mdl{sbcfwb})}
 \label{subsec:SBC_fwb}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex	(revision 11597)
@@ -9,4 +9,5 @@
 % \vfill
 % \begin{figure}[b]
+%% =================================================================================================
 % \subsubsection*{Changes record}
 % \begin{tabular}{l||l|m{0.65\linewidth}}
@@ -39,4 +40,5 @@
 $\mathbf{\xi}$ are uncorrelated over the horizontal and fully correlated along the vertical.
 
+%% =================================================================================================
 \section{Stochastic processes}
 \label{sec:STO_the_details}
@@ -114,4 +116,5 @@
 for any other configuration or resolution of the model.
 
+%% =================================================================================================
 \section{Implementation details}
 \label{sec:STO_thech_details}
@@ -165,5 +168,4 @@
 The set of parameters is available in \nam{sto}{sto} namelist
 (only the subset for equation of state stochastic parametrisation is listed below):
-%---------------------------------------namsto--------------------------------------------------
 
 \begin{listing}
@@ -172,5 +174,4 @@
   \label{lst:namsto}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 The variables of stochastic paramtetrisation itself (based on the global 2° experiments as in \cite{brankart_OM13} are:
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 11597)
@@ -54,7 +54,7 @@
 (\np{ln_tra_trd}{ln\_tra\_trd} or \np[=.true.]{ln_tra_mxl}{ln\_tra\_mxl}), as described in \autoref{chap:DIA}.
 
+%% =================================================================================================
 \section[Tracer advection (\textit{traadv.F90})]{Tracer advection (\protect\mdl{traadv})}
 \label{sec:TRA_adv}
-%------------------------------------------namtra_adv-----------------------------------------------------
 
 \begin{listing}
@@ -63,5 +63,4 @@
   \label{lst:namtra_adv}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 When considered (\ie\ when \np{ln_traadv_OFF}{ln\_traadv\_OFF} is not set to \forcode{.true.}),
@@ -171,4 +170,5 @@
 their results.
 
+%% =================================================================================================
 \subsection[CEN: Centred scheme (\forcode{ln_traadv_cen})]{CEN: Centred scheme (\protect\np{ln_traadv_cen}{ln\_traadv\_cen})}
 \label{subsec:TRA_adv_cen}
@@ -235,4 +235,5 @@
 these near boundary grid points.
 
+%% =================================================================================================
 \subsection[FCT: Flux Corrected Transport scheme (\forcode{ln_traadv_fct})]{FCT: Flux Corrected Transport scheme (\protect\np{ln_traadv_fct}{ln\_traadv\_fct})}
 \label{subsec:TRA_adv_tvd}
@@ -274,4 +275,5 @@
 while a forward scheme is used for the diffusive part.
 
+%% =================================================================================================
 \subsection[MUSCL: Monotone Upstream Scheme for Conservative Laws (\forcode{ln_traadv_mus})]{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln_traadv_mus}{ln\_traadv\_mus})}
 \label{subsec:TRA_adv_mus}
@@ -307,4 +309,5 @@
 (\np[=.true.]{ln_mus_ups}{ln\_mus\_ups}).
 
+%% =================================================================================================
 \subsection[UBS a.k.a. UP3: Upstream-Biased Scheme (\forcode{ln_traadv_ubs})]{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln_traadv_ubs}{ln\_traadv\_ubs})}
 \label{subsec:TRA_adv_ubs}
@@ -376,4 +379,5 @@
 Note the current version of \NEMO\ uses the computationally more efficient formulation \autoref{eq:TRA_adv_ubs}.
 
+%% =================================================================================================
 \subsection[QCK: QuiCKest scheme (\forcode{ln_traadv_qck})]{QCK: QuiCKest scheme (\protect\np{ln_traadv_qck}{ln\_traadv\_qck})}
 \label{subsec:TRA_adv_qck}
@@ -396,7 +400,7 @@
 %%%gmcomment   :  Cross term are missing in the current implementation....
 
+%% =================================================================================================
 \section[Tracer lateral diffusion (\textit{traldf.F90})]{Tracer lateral diffusion (\protect\mdl{traldf})}
 \label{sec:TRA_ldf}
-%-----------------------------------------nam_traldf------------------------------------------------------
 
 \begin{listing}
@@ -405,5 +409,4 @@
   \label{lst:namtra_ldf}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{tra_ldf}{tra\_ldf} namelist variables.
@@ -425,4 +428,5 @@
 the pure vertical component is split into an explicit and an implicit part \citep{lemarie.debreu.ea_OM12}.
 
+%% =================================================================================================
 \subsection[Type of operator (\forcode{ln_traldf_}\{\forcode{OFF,lap,blp}\})]{Type of operator (\protect\np{ln_traldf_OFF}{ln\_traldf\_OFF}, \protect\np{ln_traldf_lap}{ln\_traldf\_lap}, or \protect\np{ln_traldf_blp}{ln\_traldf\_blp})}
 \label{subsec:TRA_ldf_op}
@@ -456,4 +460,5 @@
 whereas the laplacian damping time scales only like $\lambda^{-2}$.
 
+%% =================================================================================================
 \subsection[Action direction (\forcode{ln_traldf_}\{\forcode{lev,hor,iso,triad}\})]{Direction of action (\protect\np{ln_traldf_lev}{ln\_traldf\_lev}, \protect\np{ln_traldf_hor}{ln\_traldf\_hor}, \protect\np{ln_traldf_iso}{ln\_traldf\_iso}, or \protect\np{ln_traldf_triad}{ln\_traldf\_triad})}
 \label{subsec:TRA_ldf_dir}
@@ -479,4 +484,5 @@
 the next two sub-sections.
 
+%% =================================================================================================
 \subsection[Iso-level (bi-)laplacian operator (\forcode{ln_traldf_iso})]{Iso-level (bi-)laplacian operator ( \protect\np{ln_traldf_iso}{ln\_traldf\_iso})}
 \label{subsec:TRA_ldf_lev}
@@ -507,4 +513,5 @@
 They are calculated in the \mdl{zpshde} module, described in \autoref{sec:TRA_zpshde}.
 
+%% =================================================================================================
 \subsection{Standard and triad (bi-)laplacian operator}
 \label{subsec:TRA_ldf_iso_triad}
@@ -512,4 +519,5 @@
 %&&    Standard rotated (bi-)laplacian operator
 %&& ----------------------------------------------
+%% =================================================================================================
 \subsubsection[Standard rotated (bi-)laplacian operator (\textit{traldf\_iso.F90})]{Standard rotated (bi-)laplacian operator (\protect\mdl{traldf\_iso})}
 \label{subsec:TRA_ldf_iso}
@@ -555,4 +563,5 @@
 %&&     Triad rotated (bi-)laplacian operator
 %&&  -------------------------------------------
+%% =================================================================================================
 \subsubsection[Triad rotated (bi-)laplacian operator (\forcode{ln_traldf_triad})]{Triad rotated (bi-)laplacian operator (\protect\np{ln_traldf_triad}{ln\_traldf\_triad})}
 \label{subsec:TRA_ldf_triad}
@@ -573,4 +582,5 @@
 %&&    Option for the rotated operators
 %&& ----------------------------------------------
+%% =================================================================================================
 \subsubsection{Option for the rotated operators}
 \label{subsec:TRA_ldf_options}
@@ -584,9 +594,8 @@
 \end{itemize}
 
+%% =================================================================================================
 \section[Tracer vertical diffusion (\textit{trazdf.F90})]{Tracer vertical diffusion (\protect\mdl{trazdf})}
 \label{sec:TRA_zdf}
-%--------------------------------------------namzdf---------------------------------------------------------
-
-%--------------------------------------------------------------------------------------------------------------
+
 
 Options are defined through the \nam{zdf}{zdf} namelist variables.
@@ -618,7 +627,9 @@
 it overcomes the stability constraint.
 
+%% =================================================================================================
 \section{External forcing}
 \label{sec:TRA_sbc_qsr_bbc}
 
+%% =================================================================================================
 \subsection[Surface boundary condition (\textit{trasbc.F90})]{Surface boundary condition (\protect\mdl{trasbc})}
 \label{subsec:TRA_sbc}
@@ -686,7 +697,7 @@
 This is the reason why the modified filter is not applied in the linear free surface case (see \autoref{chap:TD}).
 
+%% =================================================================================================
 \subsection[Solar radiation penetration (\textit{traqsr.F90})]{Solar radiation penetration (\protect\mdl{traqsr})}
 \label{subsec:TRA_qsr}
-%--------------------------------------------namqsr--------------------------------------------------------
 
 \begin{listing}
@@ -695,5 +706,4 @@
   \label{lst:namtra_qsr}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{tra_qsr}{tra\_qsr} namelist variables.
@@ -804,7 +814,7 @@
 \end{figure}
 
+%% =================================================================================================
 \subsection[Bottom boundary condition (\textit{trabbc.F90}) - \forcode{ln_trabbc})]{Bottom boundary condition (\protect\mdl{trabbc} - \protect\np{ln_trabbc}{ln\_trabbc})}
 \label{subsec:TRA_bbc}
-%--------------------------------------------nambbc--------------------------------------------------------
 
 \begin{listing}
@@ -813,5 +823,4 @@
   \label{lst:nambbc}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 \begin{figure}[!t]
   \centering
@@ -839,7 +848,7 @@
 the \ifile{geothermal\_heating} NetCDF file (\autoref{fig:TRA_geothermal}) \citep{emile-geay.madec_OS09}.
 
+%% =================================================================================================
 \section[Bottom boundary layer (\textit{trabbl.F90} - \forcode{ln_trabbl})]{Bottom boundary layer (\protect\mdl{trabbl} - \protect\np{ln_trabbl}{ln\_trabbl})}
 \label{sec:TRA_bbl}
-%--------------------------------------------nambbl---------------------------------------------------------
 
 \begin{listing}
@@ -848,5 +857,4 @@
   \label{lst:nambbl}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{bbl}{bbl} namelist variables.
@@ -872,4 +880,5 @@
 \citet{campin.goosse_T99}.
 
+%% =================================================================================================
 \subsection[Diffusive bottom boundary layer (\forcode{nn_bbl_ldf=1})]{Diffusive bottom boundary layer (\protect\np[=1]{nn_bbl_ldf}{nn\_bbl\_ldf})}
 \label{subsec:TRA_bbl_diff}
@@ -908,4 +917,5 @@
 $\overline H^\sigma$, the along bottom mean temperature, salinity and depth, respectively.
 
+%% =================================================================================================
 \subsection[Advective bottom boundary layer (\forcode{nn_bbl_adv=1,2})]{Advective bottom boundary layer (\protect\np[=1,2]{nn_bbl_adv}{nn\_bbl\_adv})}
 \label{subsec:TRA_bbl_adv}
@@ -994,7 +1004,7 @@
 It has to be used to compute the effective velocity as well as the effective overturning circulation.
 
+%% =================================================================================================
 \section[Tracer damping (\textit{tradmp.F90})]{Tracer damping (\protect\mdl{tradmp})}
 \label{sec:TRA_dmp}
-%--------------------------------------------namtra_dmp-------------------------------------------------
 
 \begin{listing}
@@ -1003,5 +1013,4 @@
   \label{lst:namtra_dmp}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 In some applications it can be useful to add a Newtonian damping term into the temperature and salinity equations:
@@ -1050,8 +1059,7 @@
 \path{./tools/DMP_TOOLS}.
 
+%% =================================================================================================
 \section[Tracer time evolution (\textit{tranxt.F90})]{Tracer time evolution (\protect\mdl{tranxt})}
 \label{sec:TRA_nxt}
-%--------------------------------------------namdom-----------------------------------------------------
-%--------------------------------------------------------------------------------------------------------------
 
 Options are defined through the \nam{dom}{dom} namelist variables.
@@ -1082,7 +1090,7 @@
 $T^{t - \rdt} = T^t$ and $T^t = T_f$.
 
+%% =================================================================================================
 \section[Equation of state (\textit{eosbn2.F90})]{Equation of state (\protect\mdl{eosbn2})}
 \label{sec:TRA_eosbn2}
-%--------------------------------------------nameos-----------------------------------------------------
 
 \begin{listing}
@@ -1091,6 +1099,6 @@
   \label{lst:nameos}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
-
+
+%% =================================================================================================
 \subsection[Equation of seawater (\forcode{ln_}\{\forcode{teos10,eos80,seos}\})]{Equation of seawater (\protect\np{ln_teos10}{ln\_teos10}, \protect\np{ln_teos80}{ln\_teos80}, or \protect\np{ln_seos}{ln\_seos})}
 \label{subsec:TRA_eos}
@@ -1212,4 +1220,5 @@
 \end{table}
 
+%% =================================================================================================
 \subsection[Brunt-V\"{a}is\"{a}l\"{a} frequency]{Brunt-V\"{a}is\"{a}l\"{a} frequency}
 \label{subsec:TRA_bn2}
@@ -1232,4 +1241,5 @@
 They are computed through \textit{eos\_rab}, a \fortran\ function that can be found in \mdl{eosbn2}.
 
+%% =================================================================================================
 \subsection{Freezing point of seawater}
 \label{subsec:TRA_fzp}
@@ -1251,4 +1261,5 @@
 a \fortran\ function that can be found in \mdl{eosbn2}.
 
+%% =================================================================================================
 %\subsection{Potential Energy anomalies}
 %\label{subsec:TRA_bn2}
@@ -1257,4 +1268,5 @@
 %
 
+%% =================================================================================================
 \section[Horizontal derivative in \textit{zps}-coordinate (\textit{zpshde.F90})]{Horizontal derivative in \textit{zps}-coordinate (\protect\mdl{zpshde})}
 \label{sec:TRA_zpshde}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 11597)
@@ -12,4 +12,5 @@
 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN.
 
+%% =================================================================================================
 \section{Vertical mixing}
 \label{sec:ZDF}
@@ -38,5 +39,4 @@
 %and thus of the formulation used (see \autoref{chap:TD}).
 
-%--------------------------------------------namzdf--------------------------------------------------------
 
 \begin{listing}
@@ -45,6 +45,6 @@
   \label{lst:namzdf}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
-
+
+%% =================================================================================================
 \subsection[Constant (\forcode{ln_zdfcst})]{Constant (\protect\np{ln_zdfcst}{ln\_zdfcst})}
 \label{subsec:ZDF_cst}
@@ -66,8 +66,8 @@
 $\sim10^{-9}~m^2.s^{-1}$ for salinity.
 
+%% =================================================================================================
 \subsection[Richardson number dependent (\forcode{ln_zdfric})]{Richardson number dependent (\protect\np{ln_zdfric}{ln\_zdfric})}
 \label{subsec:ZDF_ric}
 
-%--------------------------------------------namric---------------------------------------------------------
 
 \begin{listing}
@@ -76,5 +76,4 @@
   \label{lst:namzdf_ric}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 When \np[=.true.]{ln_zdfric}{ln\_zdfric}, a local Richardson number dependent formulation for the vertical momentum and
@@ -124,7 +123,7 @@
 the empirical values \np{rn_wtmix}{rn\_wtmix} and \np{rn_wvmix}{rn\_wvmix} \citep{lermusiaux_JMS01}.
 
+%% =================================================================================================
 \subsection[TKE turbulent closure scheme (\forcode{ln_zdftke})]{TKE turbulent closure scheme (\protect\np{ln_zdftke}{ln\_zdftke})}
 \label{subsec:ZDF_tke}
-%--------------------------------------------namzdf_tke--------------------------------------------------
 
 \begin{listing}
@@ -133,5 +132,4 @@
   \label{lst:namzdf_tke}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 The vertical eddy viscosity and diffusivity coefficients are computed from a TKE turbulent closure model based on
@@ -196,4 +194,5 @@
 \np{rn_avt0}{rn\_avt0} (\nam{zdf}{zdf} namelist, see \autoref{subsec:ZDF_cst}).
 
+%% =================================================================================================
 \subsubsection{Turbulent length scale}
 
@@ -266,6 +265,6 @@
 $\bar{e}$ reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
 
+%% =================================================================================================
 \subsubsection{Surface wave breaking parameterization}
-%-----------------------------------------------------------------------%
 
 Following \citet{mellor.blumberg_JPO04}, the TKE turbulence closure model has been modified to
@@ -300,6 +299,6 @@
 surface $\bar{e}$ value.
 
+%% =================================================================================================
 \subsubsection{Langmuir cells}
-%--------------------------------------%
 
 Langmuir circulations (LC) can be described as ordered large-scale vertical motions in
@@ -354,6 +353,6 @@
 \]
 
+%% =================================================================================================
 \subsubsection{Mixing just below the mixed layer}
-%--------------------------------------------------------------%
 
 Vertical mixing parameterizations commonly used in ocean general circulation models tend to
@@ -402,8 +401,8 @@
 % (\eg\ Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
 
+%% =================================================================================================
 \subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls})]{GLS: Generic Length Scale (\protect\np{ln_zdfgls}{ln\_zdfgls})}
 \label{subsec:ZDF_gls}
 
-%--------------------------------------------namzdf_gls---------------------------------------------------------
 
 \begin{listing}
@@ -412,5 +411,4 @@
   \label{lst:namzdf_gls}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on two prognostic equations:
@@ -463,5 +461,4 @@
 They are made available through the \np{nn_clo}{nn\_clo} namelist parameter.
 
-%--------------------------------------------------TABLE--------------------------------------------------
 \begin{table}[htbp]
   \centering
@@ -490,5 +487,4 @@
   \label{tab:ZDF_GLS}
 \end{table}
-%--------------------------------------------------------------------------------------------------------------
 
 In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of
@@ -522,7 +518,7 @@
  in \citet{reffray.guillaume.ea_GMD15} for the \NEMO\ model.
 
+%% =================================================================================================
 \subsection[OSM: OSMosis boundary layer scheme (\forcode{ln_zdfosm})]{OSM: OSMosis boundary layer scheme (\protect\np{ln_zdfosm}{ln\_zdfosm})}
 \label{subsec:ZDF_osm}
-%--------------------------------------------namzdf_osm---------------------------------------------------------
 
 \begin{listing}
@@ -531,8 +527,8 @@
   \label{lst:namzdf_osm}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 The OSMOSIS turbulent closure scheme is based on......   TBC
 
+%% =================================================================================================
 \subsection[ Discrete energy conservation for TKE and GLS schemes]{Discrete energy conservation for TKE and GLS schemes}
 \label{subsec:ZDF_tke_ene}
@@ -635,4 +631,5 @@
 %For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.
 
+%% =================================================================================================
 \section{Convection}
 \label{sec:ZDF_conv}
@@ -645,4 +642,5 @@
 or/and the use of a turbulent closure scheme.
 
+%% =================================================================================================
 \subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc})]{Non-penetrative convective adjustment (\protect\np{ln_tranpc}{ln\_tranpc})}
 \label{subsec:ZDF_npc}
@@ -705,4 +703,5 @@
 having to recompute the expansion coefficients at each mixing iteration.
 
+%% =================================================================================================
 \subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd})]{Enhanced vertical diffusion (\protect\np{ln_zdfevd}{ln\_zdfevd})}
 \label{subsec:ZDF_evd}
@@ -728,4 +727,5 @@
 a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:TD_mLF}).
 
+%% =================================================================================================
 \subsection[Handling convection with turbulent closure schemes (\forcode{ln_zdf_}\{\forcode{tke,gls,osm}\})]{Handling convection with turbulent closure schemes (\forcode{ln_zdf{tke,gls,osm}})}
 \label{subsec:ZDF_tcs}
@@ -752,11 +752,10 @@
 % gm%  + one word on non local flux with KPP scheme trakpp.F90 module...
 
+%% =================================================================================================
 \section[Double diffusion mixing (\forcode{ln_zdfddm})]{Double diffusion mixing (\protect\np{ln_zdfddm}{ln\_zdfddm})}
 \label{subsec:ZDF_ddm}
 
-%-------------------------------------------namzdf_ddm-------------------------------------------------
 %
 %\nlst{namzdf_ddm}
-%--------------------------------------------------------------------------------------------------------------
 
 This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter
@@ -840,8 +839,8 @@
 This avoids duplication in the computation of $\alpha$ and $\beta$ (which is usually quite expensive).
 
+%% =================================================================================================
 \section[Bottom and top friction (\textit{zdfdrg.F90})]{Bottom and top friction (\protect\mdl{zdfdrg})}
 \label{sec:ZDF_drg}
 
-%--------------------------------------------namdrg--------------------------------------------------------
 %
 \begin{listing}
@@ -861,5 +860,4 @@
 \end{listing}
 
-%--------------------------------------------------------------------------------------------------------------
 
 Options to define the top and bottom friction are defined through the \nam{drg}{drg} namelist variables.
@@ -916,4 +914,5 @@
 Note than from \NEMO\ 4.0, drag coefficients are only computed at cell centers (\ie\ at T-points) and refer to as $c_b^T$ in the following. These are then linearly interpolated in space to get $c_b^\textbf{U}$ at velocity points.
 
+%% =================================================================================================
 \subsection[Linear top/bottom friction (\forcode{ln_lin})]{Linear top/bottom friction (\protect\np{ln_lin}{ln\_lin})}
 \label{subsec:ZDF_drg_linear}
@@ -952,4 +951,5 @@
 $mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}.
 
+%% =================================================================================================
 \subsection[Non-linear top/bottom friction (\forcode{ln_non_lin})]{Non-linear top/bottom friction (\protect\np{ln_non_lin}{ln\_non\_lin})}
 \label{subsec:ZDF_drg_nonlinear}
@@ -984,4 +984,5 @@
 $mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}.
 
+%% =================================================================================================
 \subsection[Log-layer top/bottom friction (\forcode{ln_loglayer})]{Log-layer top/bottom friction (\protect\np{ln_loglayer}{ln\_loglayer})}
 \label{subsec:ZDF_drg_loglayer}
@@ -1007,4 +1008,5 @@
 %In this case, the relevant namelist parameters are \np{rn_tfrz0}{rn\_tfrz0}, \np{rn_tfri2}{rn\_tfri2} and \np{rn_tfri2_max}{rn\_tfri2\_max}.
 
+%% =================================================================================================
 \subsection[Explicit top/bottom friction (\forcode{ln_drgimp=.false.})]{Explicit top/bottom friction (\protect\np[=.false.]{ln_drgimp}{ln\_drgimp})}
 \label{subsec:ZDF_drg_stability}
@@ -1065,4 +1067,5 @@
 The number of potential breaches of the explicit stability criterion are still reported for information purposes.
 
+%% =================================================================================================
 \subsection[Implicit top/bottom friction (\forcode{ln_drgimp=.true.})]{Implicit top/bottom friction (\protect\np[=.true.]{ln_drgimp}{ln\_drgimp})}
 \label{subsec:ZDF_drg_imp}
@@ -1092,4 +1095,5 @@
 Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit.
 
+%% =================================================================================================
 \subsection[Bottom friction with split-explicit free surface]{Bottom friction with split-explicit free surface}
 \label{subsec:ZDF_drg_ts}
@@ -1105,8 +1109,8 @@
 Note that other strategies are possible, like considering vertical diffusion step in advance, \ie\ prior barotropic integration.
 
+%% =================================================================================================
 \section[Internal wave-driven mixing (\forcode{ln_zdfiwm})]{Internal wave-driven mixing (\protect\np{ln_zdfiwm}{ln\_zdfiwm})}
 \label{subsec:ZDF_tmx_new}
 
-%--------------------------------------------namzdf_iwm------------------------------------------
 %
 \begin{listing}
@@ -1115,5 +1119,4 @@
   \label{lst:namzdf_iwm}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 The parameterization of mixing induced by breaking internal waves is a generalization of
@@ -1166,4 +1169,5 @@
 % Jc: input files names ?
 
+%% =================================================================================================
 \section[Surface wave-induced mixing (\forcode{ln_zdfswm})]{Surface wave-induced mixing (\protect\np{ln_zdfswm}{ln\_zdfswm})}
 \label{subsec:ZDF_swm}
@@ -1196,4 +1200,5 @@
 (for more information on wave parameters and settings see \autoref{sec:SBC_wave})
 
+%% =================================================================================================
 \section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp})]{Adaptive-implicit vertical advection(\protect\np{ln_zad_Aimp}{ln\_zad\_Aimp})}
 \label{subsec:ZDF_aimp}
@@ -1319,4 +1324,5 @@
 \end{figure}
 
+%% =================================================================================================
 \subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case}
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_cfgs.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_cfgs.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_cfgs.tex	(revision 11597)
@@ -7,4 +7,5 @@
 \chaptertoc
 
+%% =================================================================================================
 \section{Introduction}
 \label{sec:CFGS_intro}
@@ -18,5 +19,4 @@
 Configuration is defined manually through the \nam{cfg}{cfg} namelist variables.
 
-%------------------------------------------namcfg----------------------------------------------------
 
 \begin{listing}
@@ -25,6 +25,6 @@
   \label{lst:namcfg}
 \end{listing}
-%-------------------------------------------------------------------------------------------------------------
-
+
+%% =================================================================================================
 \section[C1D: 1D Water column model (\texttt{\textbf{key\_c1d}})]{C1D: 1D Water column model (\protect\key{c1d})}
 \label{sec:CFGS_c1d}
@@ -64,4 +64,5 @@
 % to be added:  a test case on the yearlong Ocean Weather Station (OWS) Papa dataset of Martin (1985)
 
+%% =================================================================================================
 \section{ORCA family: global ocean with tripolar grid}
 \label{sec:CFGS_orca}
@@ -90,4 +91,5 @@
 \end{figure}
 
+%% =================================================================================================
 \subsection{ORCA tripolar grid}
 \label{subsec:CFGS_orca_grid}
@@ -128,4 +130,5 @@
 while the ratio of anisotropy remains close to one except near the Victoria Island in the Canadian Archipelago.
 
+%% =================================================================================================
 \subsection{ORCA pre-defined resolution}
 \label{subsec:CFGS_orca_resolution}
@@ -138,5 +141,4 @@
 (\autoref{tab:CFGS_ORCA}).
 
-%--------------------------------------------------TABLE--------------------------------------------------
 \begin{table}[!t]
   \centering
@@ -156,5 +158,4 @@
   \label{tab:CFGS_ORCA}
 \end{table}
-%--------------------------------------------------------------------------------------------------------------
 
 The ORCA\_R2 configuration has the following specificity: starting from a 2\deg\ ORCA mesh,
@@ -196,4 +197,5 @@
 sponge layers at open boundaries.
 
+%% =================================================================================================
 \section{GYRE family: double gyre basin}
 \label{sec:CFGS_gyre}
@@ -255,4 +257,5 @@
 \end{figure}
 
+%% =================================================================================================
 \section{AMM: atlantic margin configuration}
 \label{sec:CFGS_config_AMM}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex	(revision 11597)
@@ -39,4 +39,5 @@
 \citep{Marti1992?, Levy1996?, Levy1998?}.
 
+%% =================================================================================================
 \section{Conservation properties on ocean dynamics}
 \label{sec:CONS_Invariant_dyn}
@@ -152,4 +153,5 @@
 otherwise there is no guarantee that the surface pressure force work vanishes.
 
+%% =================================================================================================
 \section{Conservation properties on ocean thermodynamics}
 \label{sec:CONS_Invariant_tra}
@@ -170,4 +172,5 @@
 In practice, the mass is conserved with a very good accuracy.
 
+%% =================================================================================================
 \subsection{Conservation properties on momentum physics}
 \label{subsec:CONS_Invariant_dyn_physics}
@@ -273,4 +276,5 @@
 \ie\ the vertical momentum physics conserve momentum, potential vorticity, and horizontal divergence.
 
+%% =================================================================================================
 \subsection{Conservation properties on tracer physics}
 \label{subsec:CONS_Invariant_tra_physics}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_misc.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_misc.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_misc.tex	(revision 11597)
@@ -7,4 +7,5 @@
 \chaptertoc
 
+%% =================================================================================================
 \section{Representation of unresolved straits}
 \label{sec:MISC_strait}
@@ -28,4 +29,5 @@
 lateral friction.
 
+%% =================================================================================================
 \subsection{Hand made geometry changes}
 \label{subsec:MISC_strait_hand}
@@ -103,4 +105,5 @@
 \end{figure}
 
+%% =================================================================================================
 \section[Closed seas (\textit{closea.F90})]{Closed seas (\protect\mdl{closea})}
 \label{sec:MISC_closea}
@@ -167,7 +170,9 @@
 them to the domain configuration file in the utils/tools/DOMAINcfg directory.
 
+%% =================================================================================================
 \section{Sub-domain functionality}
 \label{sec:MISC_zoom}
 
+%% =================================================================================================
 \subsection{Simple subsetting of input files via NetCDF attributes}
 
@@ -230,7 +235,9 @@
 conditions. Experimenting with this remains an exercise for the user.
 
+%% =================================================================================================
 \section[Accuracy and reproducibility (\textit{lib\_fortran.F90})]{Accuracy and reproducibility (\protect\mdl{lib\_fortran})}
 \label{sec:MISC_fortran}
 
+%% =================================================================================================
 \subsection[Issues with intrinsinc SIGN function (\texttt{\textbf{key\_nosignedzero}})]{Issues with intrinsinc SIGN function (\protect\key{nosignedzero})}
 \label{subsec:MISC_sign}
@@ -258,4 +265,5 @@
 some computers/compilers.
 
+%% =================================================================================================
 \subsection{MPP reproducibility}
 \label{subsec:MISC_glosum}
@@ -287,4 +295,5 @@
 Note also that this implementation may be sensitive to the optimization level.
 
+%% =================================================================================================
 \subsection{MPP scalability}
 \label{subsec:MISC_mppsca}
@@ -311,7 +320,7 @@
 non-reference configuration.
 
+%% =================================================================================================
 \section{Model optimisation, control print and benchmark}
 \label{sec:MISC_opt}
-%--------------------------------------------namctl-------------------------------------------------------
 
 \begin{listing}
@@ -320,8 +329,8 @@
   \label{lst:namctl}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 Options are defined through the  \nam{ctl}{ctl} namelist variables.
 
+%% =================================================================================================
 \subsection{Vector optimisation}
 
@@ -334,4 +343,5 @@
 % Add also one word on NEC specific optimisation (Novercheck option for example)
 
+%% =================================================================================================
 \subsection{Control print}
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex	(revision 11597)
@@ -1,3 +1,2 @@
-
 \documentclass[../main/NEMO_manual]{subfiles}
 
@@ -10,7 +9,9 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \section{Primitive equations}
 \label{sec:MB_PE}
 
+%% =================================================================================================
 %% =================================================================================================
 \subsection{Vector invariant formulation}
@@ -91,4 +92,5 @@
 Their nature and formulation are discussed in \autoref{sec:MB_zdf_ldf} and \autoref{subsec:MB_boundary_condition}.
 
+%% =================================================================================================
 %% =================================================================================================
 \subsection{Boundary conditions}
@@ -164,7 +166,9 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \section{Horizontal pressure gradient}
 \label{sec:MB_hor_pg}
 
+%% =================================================================================================
 %% =================================================================================================
 \subsection{Pressure formulation}
@@ -200,4 +204,5 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \subsection{Free surface formulation}
 \label{subsec:MB_free_surface}
@@ -250,7 +255,9 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \section{Curvilinear \textit{z-}coordinate system}
 \label{sec:MB_zco}
 
+%% =================================================================================================
 %% =================================================================================================
 \subsection{Tensorial formalism}
@@ -337,4 +344,5 @@
 where $q$ is a scalar quantity and $\vect A = (a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinates system.
 
+%% =================================================================================================
 %% =================================================================================================
 \subsection{Continuous model equations}
@@ -522,4 +530,5 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \section{Curvilinear generalised vertical coordinate system}
 \label{sec:MB_gco}
@@ -602,4 +611,5 @@
 %}
 
+%% =================================================================================================
 %% =================================================================================================
 \subsection{\textit{S}-coordinate formulation}
@@ -691,4 +701,5 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \subsection{Curvilinear \zstar-coordinate system}
 \label{subsec:MB_zco_star}
@@ -777,7 +788,9 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \subsection{Curvilinear terrain-following \textit{s}--coordinate}
 \label{subsec:MB_sco}
 
+%% =================================================================================================
 %% =================================================================================================
 \subsubsection{Introduction}
@@ -862,4 +875,5 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \subsection{\texorpdfstring{Curvilinear \ztilde-coordinate}{}}
 \label{subsec:MB_zco_tilde}
@@ -870,4 +884,5 @@
 Its use is therefore not recommended.
 
+%% =================================================================================================
 %% =================================================================================================
 \section{Subgrid scale physics}
@@ -894,4 +909,5 @@
 The formulation of these terms and their underlying physics are briefly discussed in the next two subsections.
 
+%% =================================================================================================
 %% =================================================================================================
 \subsection{Vertical subgrid scale physics}
@@ -927,4 +943,5 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \subsection{Formulation of the lateral diffusive and viscous operators}
 \label{subsec:MB_ldf}
@@ -981,4 +998,5 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \subsubsection{Lateral laplacian tracer diffusive operator}
 
@@ -1022,4 +1040,5 @@
 while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$.
 
+%% =================================================================================================
 %% =================================================================================================
 \subsubsection{Eddy induced velocity}
@@ -1060,4 +1079,5 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \subsubsection{Lateral bilaplacian tracer diffusive operator}
 
@@ -1071,4 +1091,5 @@
 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one.
 
+%% =================================================================================================
 %% =================================================================================================
 \subsubsection{Lateral Laplacian momentum diffusive operator}
@@ -1105,4 +1126,5 @@
 
 %% =================================================================================================
+%% =================================================================================================
 \subsubsection{Lateral bilaplacian momentum diffusive operator}
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex	(revision 11597)
@@ -3,4 +3,5 @@
 \begin{document}
 \chapter{ essai \zstar \sstar}
+%% =================================================================================================
 \section{Curvilinear \zstar- or \sstar coordinate system}
 
@@ -60,10 +61,9 @@
 %%%
 
+%% =================================================================================================
 \section[Surface pressure gradient and sea surface heigth (\textit{dynspg.F90})]{Surface pressure gradient and sea surface heigth (\protect\mdl{dynspg})}
 \label{sec:MBZ_dyn_hpg_spg}
-%-----------------------------------------nam_dynspg----------------------------------------------------
 
 %\nlst{nam_dynspg}
-%------------------------------------------------------------------------------------------------------------
 Options are defined through the \nam{_dynspg}{\_dynspg} namelist variables.
 The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:MB_hor_pg}).
@@ -81,7 +81,6 @@
 so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
 
-%-------------------------------------------------------------
 % Explicit
-%-------------------------------------------------------------
+%% =================================================================================================
 \subsubsection[Explicit (\texttt{\textbf{key\_dynspg\_exp}})]{Explicit (\protect\key{dynspg\_exp})}
 \label{subsec:MBZ_dyn_spg_exp}
@@ -117,12 +116,9 @@
 (\autoref{eq:DYN_spg_exp}).
 
-%-------------------------------------------------------------
 % Split-explicit time-stepping
-%-------------------------------------------------------------
+%% =================================================================================================
 \subsubsection[Split-explicit time-stepping (\texttt{\textbf{key\_dynspg\_ts}})]{Split-explicit time-stepping (\protect\key{dynspg\_ts})}
 \label{subsec:MBZ_dyn_spg_ts}
-%--------------------------------------------namdom----------------------------------------------------
-
-%--------------------------------------------------------------------------------------------------------------
+
 The split-explicit free surface formulation used in OPA follows the one proposed by \citet{Griffies2004?}.
 The general idea is to solve the free surface equation with a small time step,
@@ -274,7 +270,6 @@
 be more conservative, and so is recommended.
 
-%-------------------------------------------------------------
 % Filtered formulation
-%-------------------------------------------------------------
+%% =================================================================================================
 \subsubsection[Filtered formulation (\texttt{\textbf{key\_dynspg\_flt}})]{Filtered formulation (\protect\key{dynspg\_flt})}
 \label{subsec:MBZ_dyn_spg_flt}
@@ -288,7 +283,6 @@
 \colorbox{red}{\np[=1]{rnu}{rnu} to be suppressed from namelist !}
 
-%-------------------------------------------------------------
 % Non-linear free surface formulation
-%-------------------------------------------------------------
+%% =================================================================================================
 \subsection[Non-linear free surface formulation (\texttt{\textbf{key\_vvl}})]{Non-linear free surface formulation (\protect\key{vvl})}
 \label{subsec:MBZ_dyn_spg_vvl}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex	(revision 11596)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex	(revision 11597)
@@ -13,4 +13,5 @@
   would help  ==> to be added}
 %%%%
+
 
 Having defined the continuous equations in \autoref{chap:MB}, we need now to choose a time discretization,
@@ -20,4 +21,5 @@
 the consequences for the order in which the equations are solved.
 
+%% =================================================================================================
 \section{Time stepping environment}
 \label{sec:TD_environment}
@@ -47,4 +49,5 @@
 The time stepping itself is performed once at each time step where implicit vertical diffusion is computed, \ie\ in the \mdl{trazdf} and \mdl{dynzdf} modules.
 
+%% =================================================================================================
 \section{Non-diffusive part --- Leapfrog scheme}
 \label{sec:TD_leap_frog}
@@ -87,4 +90,5 @@
 filter parameter and the viscosity and diffusion coefficients.
 
+%% =================================================================================================
 \section{Diffusive part --- Forward or backward scheme}
 \label{sec:TD_forward_imp}
@@ -151,4 +155,5 @@
 (see for example \citet{richtmyer.morton_bk67}).
 
+%% =================================================================================================
 \section{Surface pressure gradient}
 \label{sec:TD_spg_ts}
@@ -181,4 +186,5 @@
 %}
 
+%% =================================================================================================
 \section{Modified Leapfrog -- Asselin filter scheme}
 \label{sec:TD_mLF}
@@ -239,8 +245,8 @@
 \end{figure}
 
+%% =================================================================================================
 \section{Start/Restart strategy}
 \label{sec:TD_rst}
 
-%--------------------------------------------namrun-------------------------------------------
 \begin{listing}
   \nlst{namrun}
@@ -248,5 +254,4 @@
   \label{lst:namrun}
 \end{listing}
-%--------------------------------------------------------------------------------------------------------------
 
 The first time step of this three level scheme when starting from initial conditions is a forward step
@@ -286,4 +291,75 @@
 \gmcomment{       % add a subsection here
 
-%-------------------------------------------------------------------------------------------------------------
 %        Time Domain
+% ------------------------------------------------------------------------------------------------------------
+%% =================================================================================================
+\subsection{Time domain}
+\label{subsec:TD_time}
+
+Options are defined through the  \nam{dom}{dom} namelist variables.
+ \colorbox{yellow}{add here a few word on nit000 and nitend}
+
+ \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj}
+
+add a description of daymod, and the model calandar (leap-year and co)
+
+}     	%% end add
+
+
+
+%%
+\gmcomment{       % add implicit in vvl case  and Crant-Nicholson scheme
+
+Implicit time stepping in case of variable volume thickness.
+
+Tracer case (NB for momentum in vector invariant form take care!)
+
+\begin{flalign*}
+  &\frac{\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}}{2\rdt}
+  \equiv \text{RHS}+ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]}
+  \rt]      \\
+  &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}
+  \equiv {2\rdt} \ \text{RHS}+ {2\rdt} \ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]}
+  \rt]      \\
+  &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}
+  \equiv 2\rdt \ \text{RHS}
+  + 2\rdt \ \lt\{ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} [ T_{k +1}^{t+1} - T_k      ^{t+1} ]
+    - \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} [ T_k       ^{t+1} - T_{k -1}^{t+1} ]  \rt\}     \\
+  &\\
+  &\lt( e_{3t}\,T \rt)_k^{t+1}
+  -  {2\rdt} \           \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}                  T_{k +1}^{t+1}
+  + {2\rdt} \ \lt\{  \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}
+    +  \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}     \rt\}   T_{k    }^{t+1}
+  -  {2\rdt} \           \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}                  T_{k -1}^{t+1}      \\
+  &\equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}    \\
+  %
+\end{flalign*}
+\begin{flalign*}
+  \allowdisplaybreaks
+  \intertext{ Tracer case }
+  %
+  &  \qquad \qquad  \quad   -  {2\rdt}                  \ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}
+  \qquad \qquad \qquad  \qquad  T_{k +1}^{t+1}   \\
+  &+ {2\rdt} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}
+  +   \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\
+  & \qquad \qquad  \qquad \qquad \qquad \quad \ \ -  {2\rdt} \                          \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}                          \quad \ \ T_{k -1}^{t+1}
+  \ \equiv \ \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}  \\
+  %
+\end{flalign*}
+\begin{flalign*}
+  \allowdisplaybreaks
+  \intertext{ Tracer content case }
+  %
+  & -  {2\rdt} \              & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_{k +1}^{t+1}}  && \  \lt( e_{3t}\,T \rt)_{k +1}^{t+1}   &\\
+  & + {2\rdt} \ \lt[ 1  \rt.+ & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_k^{t+1}}
+  + & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_k^{t+1}}  \lt.  \rt]  & \lt( e_{3t}\,T \rt)_{k   }^{t+1}  &\\
+  & -  {2\rdt} \               & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_{k -1}^{t+1}}     &\  \lt( e_{3t}\,T \rt)_{k -1}^{t+1}
+  \equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}  &
+\end{flalign*}
+
+%%
+}
+
+\onlyinsubfile{\input{../../global/epilogue}}
+
+\end{document}
