Changeset 11543

NEMO/trunk/doc/latex/NEMO/main/appendices.tex

-                      r11330
+                      r11543
 \subfile{../subfiles/annex_A}             %% Generalised vertical coordinate
 \subfile{../subfiles/annex_B}             %% Diffusive operator
 \subfile{../subfiles/annex_C}             %% Discrete invariants of the eqs.
 \subfile{../subfiles/annex_iso}            %% Isoneutral diffusion using triads
 \subfile{../subfiles/annex_DOMAINcfg}     %% Brief notes on DOMAINcfg
+\subfile{../subfiles/apdx_s_coord}      %% Generalised vertical coordinate
+\subfile{../subfiles/apdx_diff_opers}   %% Diffusive operators
+\subfile{../subfiles/apdx_invariants}   %% Discrete invariants of the eqs.
+\subfile{../subfiles/apdx_triads}       %% Isoneutral diffusion using triads
+\subfile{../subfiles/apdx_DOMAINcfg}    %% Brief notes on DOMAINcfg
 %% Not included
 …
 %\subfile{../subfiles/chap_DIU}
 %\subfile{../subfiles/chap_conservation}
+%\subfile{../subfiles/annex_E}            %% Notes on some on going staff
+%\subfile{../subfiles/apdx_algos}   %% Notes on some on going staff

NEMO/trunk/doc/latex/NEMO/main/chapters.tex

r11522	r11543
13	13	\subfile{../subfiles/chap_STO} %% Stochastic param.
14	14	\subfile{../subfiles/chap_misc} %% Miscellaneous topics
15		\subfile{../subfiles/chap_~~CONFIG}~~ %% Predefined configurations
	15	\subfile{../subfiles/chap_cfgs} %% Predefined configurations
16	16
17	17	%% Not included

NEMO/trunk/doc/latex/NEMO/main/introduction.tex

-                      r11522
+                      r11543
 \begin{description}
 \item [\nameref{chap:PE}] presents the equations and their assumptions, the vertical coordinates used,
+\item [\nameref{chap:MB}] presents the equations and their assumptions, the vertical coordinates used,
 and the subgrid scale physics.
 The equations are written in a curvilinear coordinate system, with a choice of vertical coordinates
 …
 Dimensional units in the meter, kilogram, second (MKS) international system are used throughout.
 The following chapters deal with the discrete equations.
 \item [\nameref{chap:STP}] presents the model time stepping environment.
+\item [\nameref{chap:TD}] presents the model time stepping environment.
 it is a three level scheme in which the tendency terms of the equations are evaluated either
 centered in time, or forward, or backward depending of the nature of the term.
 …
 \item [\nameref{chap:ASM}] describes how increments produced by
 data \textbf{A}s\textbf{S}i\textbf{M}ilation may be applied to the model equations.
+\item [\nameref{chap:STO}]
 \item [\nameref{chap:MISC}] (including solvers)
 \item [\nameref{chap:CFG}] provides finally a brief introduction to
+\item [\nameref{chap:CFGS}] provides finally a brief introduction to
 the pre-defined model configurations
 (water column model \texttt{C1D}, ORCA and GYRE families of configurations).
 …
 \begin{description}
 \item [\nameref{apdx:s_coord}]
 \item [\nameref{apdx:diff_oper}]
 \item [\nameref{apdx:invariants}]
 \item [\nameref{apdx:triads}]
 \item [\nameref{apdx:DOMAINcfg}]
 \item [\nameref{apdx:coding}]
+\item [\nameref{apdx:SCOORD}]
+\item [\nameref{apdx:DIFFOPERS}]
+\item [\nameref{apdx:INVARIANTS}]
+\item [\nameref{apdx:TRIADS}]
+\item [\nameref{apdx:DOMCFG}]
+\item [\nameref{apdx:CODING}]
 \end{description}

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex

-                      r11529
+                      r11543
 % ================================================================
 \chapter{A brief guide to the DOMAINcfg tool}
 \label{apdx:DOMAINcfg}
+\label{apdx:DOMCFG}
 \chaptertoc
 …
 The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and
 $gdepw_0$ for $t$- and $w$-points, respectively. See \autoref{sec:DOMCFG_sco} for the
 S-coordinate options.  As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of
+S-coordinate options.  As indicated on \autoref{fig:DOM_index_vert} \jp{jpk} is the number of
 $w$-levels.  $gdepw_0(1)$ is the ocean surface.  There are at most \jp{jpk}-1 $t$-points
 inside the ocean, the additional $t$-point at $jk = jpk$ is below the sea floor and is not
 …
 The depth field $h$ is not necessary the ocean depth,
 since a mixed step-like and bottom-following representation of the topography can be used
 (\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}).
+(\autoref{fig:DOM_z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:DOM_z_zps_s_sps}).
 The namelist parameter \np{rn\_rmax} determines the slope at which
 the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate.
 …
 \[
   z = s_{min} + C (s) (H - s_{min})
   % \label{eq:SH94_1}
+  % \label{eq:DOMCFG_SH94_1}
 \]
 …
          + b       \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] -   \tanh \lt( \frac{\theta}{2} \rt)}
                         {                                                  2 \tanh \lt( \frac{\theta}{2} \rt)}
  \label{eq:SH94_2}
+ \label{eq:DOMCFG_SH94_2}
 \]
 …
     \includegraphics[width=\textwidth]{Fig_sco_function}
     \caption{
       \protect\label{fig:sco_function}
+      \protect\label{fig:DOMCFG_sco_function}
       Examples of the stretching function applied to a seamount;
       from left to right: surface, surface and bottom, and bottom intensified resolutions
 …
 bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$.
 $b$ has been designed to allow surface and/or bottom increase of the vertical resolution
 (\autoref{fig:sco_function}).
+(\autoref{fig:DOMCFG_sco_function}).
 Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in
 …
 \begin{equation}
   z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1
   % \label{eq:z}
+  % \label{eq:DOMCFG_z}
 \end{equation}
 …
     For clarity every third coordinate surface is shown.
+  }
   \label{fig:fig_compare_coordinates_surface}
+  \label{fig:DOMCFG_fig_compare_coordinates_surface}
 \end{figure}
  % >>>>>>>>>>>>>>>>>>>>>>>>>>>>

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex

-                      r11529
+                      r11543
 % ================================================================
 \chapter{Note on some algorithms}
 \label{apdx:E}
+\label{apdx:ALGOS}
 \chaptertoc
 …
 \newpage
 This appendix some on going consideration on algorithms used or planned to be used in \NEMO.
+This appendix some on going consideration on algorithms used or planned to be used in \NEMO.
 % -------------------------------------------------------------------------------------------------------------
 %        UBS scheme
+%        UBS scheme
 % -------------------------------------------------------------------------------------------------------------
 \section{Upstream Biased Scheme (UBS) (\protect\np{ln\_traadv\_ubs}\forcode{ = .true.})}
 …
 a constant i-grid spacing ($\Delta i=1$).
 Alternative choice: introduce the scale factors:
+Alternative choice: introduce the scale factors:
 $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta_i \left[ \frac{e_{2u} e_{3u} }{e_{1u} }\delta_{i+1/2}[\tau] \right]$.
 …
 It is not a \emph{positive} scheme meaning false extrema are permitted but
 the amplitude of such are significantly reduced over the centred second order method.
 Nevertheless it is not recommended to apply it to a passive tracer that requires positivity.
+Nevertheless it is not recommended to apply it to a passive tracer that requires positivity.
 The intrinsic diffusion of UBS makes its use risky in the vertical direction where
 …
 \np{ln\_traadv\_ubs}\forcode{ = .true.}.
 For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds to
+For stability reasons, in \autoref{eq:TRA_adv_ubs}, the first term which corresponds to
 a second order centred scheme is evaluated using the \textit{now} velocity (centred in time) while
 the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity
 …
 This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme.
 UBS and QUICK schemes only differ by one coefficient.
 Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
+Substituting 1/6 with 1/8 in (\autoref{eq:TRA_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
 This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded.
 Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme.
 …
 Computer time can be saved by using a time-splitting technique on vertical advection.
 This possibility have been implemented and validated in ORCA05-L301.
 It is not currently offered in the current reference version.
+It is not currently offered in the current reference version.
 NB 2: In a forthcoming release four options will be proposed for the vertical component used in the UBS scheme.
 …
 The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme.
 NB 3: It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows:
+NB 3: It is straight forward to rewrite \autoref{eq:TRA_adv_ubs} as follows:
 \begin{equation}
   \label{eq:tra_adv_ubs2}
 …
   \right.
 \end{equation}
 or equivalently
+or equivalently
 \begin{equation}
   \label{eq:tra_adv_ubs2}
 …
   \end{split}
 \end{equation}
 \autoref{eq:tra_adv_ubs2} has several advantages.
+\autoref{eq:TRA_adv_ubs2} has several advantages.
 First it clearly evidences that the UBS scheme is based on the fourth order scheme to which
 is added an upstream biased diffusive term.
 Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step,
 not only the $2^{th}$ order part as stated above using \autoref{eq:tra_adv_ubs}.
+not only the $2^{th}$ order part as stated above using \autoref{eq:TRA_adv_ubs}.
 Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient which
 is simply proportional to the velocity.
 …
   \end{split}
 \end{equation}
 with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$,
+with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$,
 \ie\ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$
 it comes:
 …
 % -------------------------------------------------------------------------------------------------------------
 %        Leap-Frog energetic
+%        Leap-Frog energetic
 % -------------------------------------------------------------------------------------------------------------
 \section{Leapfrog energetic}
 …
   \equiv \frac{1}{\rdt} \overline{ \delta_{t+\rdt/2}[q]}^{\,t}
   =         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt}
 \]
+\]
 Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$,
 not $2\rdt$ as it can be found sometimes in literature.
 …
 \]
 is satisfied in discrete form.
 Indeed,
+Indeed,
 \[
   \begin{split}
 …
 \]
 NB here pb of boundary condition when applying the adjoint!
 In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition
+In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition
 (equivalently of the boundary value of the integration by part).
 In time this boundary condition is not physical and \textbf{add something here!!!}
 % ================================================================
 % Iso-neutral diffusion :
+% Iso-neutral diffusion :
 % ================================================================
 …
 % ================================================================
 % Griffies' iso-neutral diffusion operator :
+% Griffies' iso-neutral diffusion operator :
 % ================================================================
 \subsection{Griffies iso-neutral diffusion operator}
 …
 but is formulated within the \NEMO\ framework
 (\ie\ using scale factors rather than grid-size and having a position of $T$-points that
 is not necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}).
 In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO,
+is not necessary in the middle of vertical velocity points, see \autoref{fig:DOM_zgr_e3}).
+In the formulation \autoref{eq:TRA_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO,
 the off-diagonal terms of the small angle diffusion tensor contain several double spatial averages of a gradient,
 for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$.
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 The four iso-neutral fluxes associated with the triads are defined at $T$-point.
+The four iso-neutral fluxes associated with the triads are defined at $T$-point.
 They take the following expression:
 \begin{flalign*}
 …
 The resulting iso-neutral fluxes at $u$- and $w$-points are then given by
 the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}):
+the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:TRIADS_ISO_triad}):
 \begin{flalign}
   \label{eq:iso_flux}
 …
     + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]   \right\}
 \end{equation}
 where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells.
+where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells.
 This expression of the iso-neutral diffusion has been chosen in order to satisfy the following six properties:
 …
 % ================================================================
 % Skew flux formulation for Eddy Induced Velocity :
+% Skew flux formulation for Eddy Induced Velocity :
 % ================================================================
 \subsection{Eddy induced velocity and skew flux formulation}
 …
 the formulation of which depends on the slopes of iso-neutral surfaces.
 Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces,
 \ie\ \autoref{eq:ldfslp_geo} is used in $z$-coordinate,
 and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates.
 The eddy induced velocity is given by:
+\ie\ \autoref{eq:LDF_slp_geo} is used in $z$-coordinate,
+and the sum \autoref{eq:LDF_slp_geo} + \autoref{eq:LDF_slp_iso} in $z^*$ or $s$-coordinates.
+The eddy induced velocity is given by:
 \begin{equation}
   \label{eq:eiv_v}
 …
 (see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme.
 This is particularly useful for passive tracers where
 \emph{positivity} of the advection scheme is of paramount importance.
 % give here the expression using the triads. It is different from the one given in \autoref{eq:ldfeiv}
+\emph{positivity} of the advection scheme is of paramount importance.
+% give here the expression using the triads. It is different from the one given in \autoref{eq:LDF_eiv}
 % see just below a copy of this equation:
 %\begin{equation} \label{eq:ldfeiv}
 …
   \right)
 \end{equation}
 Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells.
+Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells.
 In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces must be added to
 $\mathbb{R}$ for the discret form to be exact.
+$\mathbb{R}$ for the discret form to be exact.
 Such a choice of discretisation is consistent with the iso-neutral operator as
 …
 $\ $\newpage      %force an empty line
 % ================================================================
 % Discrete Invariants of the iso-neutral diffrusion
+% Discrete Invariants of the iso-neutral diffrusion
 % ================================================================
 \subsection{Discrete invariants of the iso-neutral diffrusion}
 \label{subsec:Gf_operator}
 Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane.
+Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane.
 This part will be moved in an Appendix.
 …
   \int_D  D_l^T \; T \;dv   \leq 0
 \]
 The discrete form of its left hand side is obtained using \autoref{eq:iso_flux}
+The discrete form of its left hand side is obtained using \autoref{eq:TRIADS_iso_flux}
 \begin{align*}
 …
         \right\}
         \quad   \leq 0
 \end{align*}
+\end{align*}
 The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities.
 …
+                             %
                            &\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\}
 \end{align*}
+\end{align*}
 This means that the iso-neutral operator is self-adjoint.
 There is no need to develop a specific to obtain it.
 …
 \label{subsec:eiv_skew}
 Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.
+Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.
 This have to be moved in an Appendix.
 …
     &{\ \ \;_i^k  \mathbb{R}_{+1/2}^{+1/2}}   &\delta_{i+1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}]
     &\Bigr\}  \\
   \end{matrix}
+  \end{matrix}
 \end{align*}
 The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs,
 they cancel out.
+they cancel out.
 Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$.
 The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the same but both of opposite signs and

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_diff_opers.tex

-                      r11529
+                      r11543
 % Chapter Appendix B : Diffusive Operators
 % ================================================================
 \chapter{Appendix B : Diffusive Operators}
 \label{apdx:B}
+\chapter{Diffusive Operators}
+\label{apdx:DIFFOPERS}
 \chaptertoc
 …
 % ================================================================
 \section{Horizontal/Vertical $2^{nd}$ order tracer diffusive operators}
 \label{sec:B_1}
+\label{sec:DIFFOPERS_1}
 \subsubsection*{In z-coordinates}
 …
 In $z$-coordinates, the horizontal/vertical second order tracer diffusion operator is given by:
 \begin{align}
   \label{apdx:B1}
+  \label{eq:DIFFOPERS_1}
   &D^T = \frac{1}{e_1 \, e_2}      \left[
     \left. \frac{\partial}{\partial i} \left(   \frac{e_2}{e_1}A^{lT} \;\left. \frac{\partial T}{\partial i} \right|_z   \right)   \right|_z      \right.
 …
 \subsubsection*{In generalized vertical coordinates}
 In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and $\sigma_2$ by \autoref{apdx:A_s_slope} and
+In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and $\sigma_2$ by \autoref{eq:SCOORD_s_slope} and
 the vertical/horizontal ratio of diffusion coefficient by $\epsilon = A^{vT} / A^{lT}$.
 The diffusion operator is given by:
 \begin{equation}
   \label{apdx:B2}
+  \label{eq:DIFFOPERS_2}
   D^T = \left. \nabla \right|_s \cdot
   \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_s T  \right] \\
 …
   \begin{array}{*{20}l}
     D^T= \frac{1}{e_1\,e_2\,e_3 } & \left\{ \quad \quad \frac{\partial }{\partial i}  \left. \left[  e_2\,e_3 \, A^{lT}
                                \left( \  \frac{1}{e_1}\; \left. \frac{\partial T}{\partial i} \right|_s
+                               \left( \  \frac{1}{e_1}\; \left. \frac{\partial T}{\partial i} \right|_s
                                        -\frac{\sigma_1 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\
         &  \quad \  +   \            \left.   \frac{\partial }{\partial j}  \left. \left[  e_1\,e_3 \, A^{lT}
                                \left( \ \frac{1}{e_2 }\; \left. \frac{\partial T}{\partial j} \right|_s
+                               \left( \ \frac{1}{e_2 }\; \left. \frac{\partial T}{\partial j} \right|_s
                                        -\frac{\sigma_2 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\
         &  \quad \  +   \           \left.  e_1\,e_2\, \frac{\partial }{\partial s}  \left[ A^{lT} \; \left(
                      -\frac{\sigma_1 }{e_1 } \; \left. \frac{\partial T}{\partial i} \right|_s
                      -\frac{\sigma_2 }{e_2 } \; \left. \frac{\partial T}{\partial j} \right|_s
+        &  \quad \  +   \           \left.  e_1\,e_2\, \frac{\partial }{\partial s}  \left[ A^{lT} \; \left(
+                     -\frac{\sigma_1 }{e_1 } \; \left. \frac{\partial T}{\partial i} \right|_s
+                     -\frac{\sigma_2 }{e_2 } \; \left. \frac{\partial T}{\partial j} \right|_s
                           +\left( \varepsilon +\sigma_1^2+\sigma_2 ^2 \right) \; \frac{1}{e_3 } \; \frac{\partial T}{\partial s} \right) \; \right] \;  \right\} .
   \end{array}
 …
 \end{align*}
 \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption.
+\autoref{eq:DIFFOPERS_2} is obtained from \autoref{eq:DIFFOPERS_1} without any additional assumption.
 Indeed, for the special case $k=z$ and thus $e_3 =1$,
 we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A} and
 use \autoref{apdx:A_s_slope} and \autoref{apdx:A_s_chain_rule}.
 Since no cross horizontal derivative $\partial _i \partial _j $ appears in \autoref{apdx:B1},
+we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:SCOORD} and
+use \autoref{eq:SCOORD_s_slope} and \autoref{eq:SCOORD_s_chain_rule}.
+Since no cross horizontal derivative $\partial _i \partial _j $ appears in \autoref{eq:DIFFOPERS_1},
 the ($i$,$z$) and ($j$,$z$) planes are independent.
 The derivation can then be demonstrated for the ($i$,$z$)~$\to$~($j$,$s$) transformation without
 …
 % ================================================================
 \section{Iso/Diapycnal $2^{nd}$ order tracer diffusive operators}
 \label{sec:B_2}
+\label{sec:DIFFOPERS_2}
 \subsubsection*{In z-coordinates}
 …
 \begin{equation}
   \label{apdx:B3}
+  \label{eq:DIFFOPERS_3}
   \textbf {A}_{\textbf I} = \frac{A^{lT}}{\left( {1+a_1 ^2+a_2 ^2} \right)}
   \left[ {{
 …
 In practice, $\epsilon$ is small and isopycnal slopes are generally less than $10^{-2}$ in the ocean,
 so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{cox_OM87}. Keeping leading order terms\footnote{Apart from the (1,0)
+so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{cox_OM87}. Keeping leading order terms\footnote{Apart from the (1,0)
 and (0,1) elements which are set to zero. See \citet{griffies_bk04}, section 14.1.4.1 for a discussion of this point.}:
 \begin{subequations}
   \label{apdx:B4}
+  \label{eq:DIFFOPERS_4}
   \begin{equation}
     \label{apdx:B4a}
+    \label{eq:DIFFOPERS_4a}
     {\textbf{A}_{\textbf{I}}} \approx A^{lT}\;\Re\;\text{where} \;\Re =
     \left[ {{
 …
   and the iso/dianeutral diffusive operator in $z$-coordinates is then
   \begin{equation}
     \label{apdx:B4b}
+    \label{eq:DIFFOPERS_4b}
     D^T = \left. \nabla \right|_z \cdot
     \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_z T  \right]. \\
 …
 \end{subequations}
 Physically, the full tensor \autoref{apdx:B3} represents strong isoneutral diffusion on a plane parallel to
+Physically, the full tensor \autoref{eq:DIFFOPERS_3} represents strong isoneutral diffusion on a plane parallel to
 the isoneutral surface and weak dianeutral diffusion perpendicular to this plane.
 However,
 the approximate `weak-slope' tensor \autoref{apdx:B4a} represents strong diffusion along the isoneutral surface,
+the approximate `weak-slope' tensor \autoref{eq:DIFFOPERS_4a} represents strong diffusion along the isoneutral surface,
 with weak \emph{vertical} diffusion -- the principal axes of the tensor are no longer orthogonal.
 This simplification also decouples the ($i$,$z$) and ($j$,$z$) planes of the tensor.
 The weak-slope operator therefore takes the same form, \autoref{apdx:B4}, as \autoref{apdx:B2},
+The weak-slope operator therefore takes the same form, \autoref{eq:DIFFOPERS_4}, as \autoref{eq:DIFFOPERS_2},
 the diffusion operator for geopotential diffusion written in non-orthogonal $i,j,s$-coordinates.
 Written out explicitly,
 \begin{multline}
   \label{apdx:B_ldfiso}
+  \label{eq:DIFFOPERS_ldfiso}
   D^T=\frac{1}{e_1 e_2 }\left\{
     {\;\frac{\partial }{\partial i}\left[ {A_h \left( {\frac{e_2}{e_1}\frac{\partial T}{\partial i}-a_1 \frac{e_2}{e_3}\frac{\partial T}{\partial k}} \right)} \right]}
 …
 \end{multline}
 The isopycnal diffusion operator \autoref{apdx:B4},
 \autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its square.
 As \autoref{apdx:B4} is the divergence of a flux, the demonstration of the first property is trivial, providing that the flux normal to the boundary is zero
+The isopycnal diffusion operator \autoref{eq:DIFFOPERS_4},
+\autoref{eq:DIFFOPERS_ldfiso} conserves tracer quantity and dissipates its square.
+As \autoref{eq:DIFFOPERS_4} is the divergence of a flux, the demonstration of the first property is trivial, providing that the flux normal to the boundary is zero
 (as it is when $A_h$ is zero at the boundary). Let us demonstrate the second one:
 \[
 …
              j}-a_2 \frac{\partial T}{\partial k}} \right)^2}
              +\varepsilon \left(\frac{\partial T}{\partial k}\right) ^2\right]      \\
            & \geq 0 .
+           & \geq 0 .
   \end{array}
+             }
 …
 \subsubsection*{In generalized vertical coordinates}
 Because the weak-slope operator \autoref{apdx:B4},
 \autoref{apdx:B_ldfiso} is decoupled in the ($i$,$z$) and ($j$,$z$) planes,
+Because the weak-slope operator \autoref{eq:DIFFOPERS_4},
+\autoref{eq:DIFFOPERS_ldfiso} is decoupled in the ($i$,$z$) and ($j$,$z$) planes,
 it may be transformed into generalized $s$-coordinates in the same way as
 \autoref{sec:B_1} was transformed into \autoref{sec:B_2}.
+\autoref{sec:DIFFOPERS_1} was transformed into \autoref{sec:DIFFOPERS_2}.
 The resulting operator then takes the simple form
 \begin{equation}
   \label{apdx:B_ldfiso_s}
+  \label{eq:DIFFOPERS_ldfiso_s}
   D^T = \left. \nabla \right|_s \cdot
   \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_s T  \right] \\
 …
 \]
 To prove \autoref{apdx:B_ldfiso_s} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious.
 An easier way is first to note (by reversing the derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that
+To prove \autoref{eq:DIFFOPERS_ldfiso_s} by direct re-expression of \autoref{eq:DIFFOPERS_ldfiso} is straightforward, but laborious.
+An easier way is first to note (by reversing the derivation of \autoref{sec:DIFFOPERS_2} from \autoref{sec:DIFFOPERS_1} ) that
 the weak-slope operator may be \emph{exactly} reexpressed in non-orthogonal $i,j,\rho$-coordinates as
 \begin{equation}
   \label{apdx:B5}
+  \label{eq:DIFFOPERS_5}
   D^T = \left. \nabla \right|_\rho \cdot
   \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_\rho T  \right] \\
 …
 \end{equation}
 Then direct transformation from $i,j,\rho$-coordinates to $i,j,s$-coordinates gives
 \autoref{apdx:B_ldfiso_s} immediately.
+\autoref{eq:DIFFOPERS_ldfiso_s} immediately.
 Note that the weak-slope approximation is only made in transforming from
 …
 The further transformation into $i,j,s$-coordinates is exact, whatever the steepness of the $s$-surfaces,
 in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates in
 \autoref{sec:B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces.
+\autoref{sec:DIFFOPERS_1} onto $s$-coordinates is exact, however steep the $s$-surfaces.
 …
 % ================================================================
 \section{Lateral/Vertical momentum diffusive operators}
 \label{sec:B_3}
+\label{sec:DIFFOPERS_3}
 The second order momentum diffusion operator (Laplacian) in $z$-coordinates is found by
 applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian of a vector,
+applying \autoref{eq:MB_lap_vector}, the expression for the Laplacian of a vector,
 to the horizontal velocity vector:
 \begin{align*}
 …
   }} \right)
 \end{align*}
 Using \autoref{eq:PE_div}, the definition of the horizontal divergence,
+Using \autoref{eq:MB_div}, the definition of the horizontal divergence,
 the third component of the second vector is obviously zero and thus :
 \[
   \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \textbf{k} \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) .
+  \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \textbf{k} \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) .
 \]
 Note that this operator ensures a full separation between
 the vorticity and horizontal divergence fields (see \autoref{apdx:C}).
+the vorticity and horizontal divergence fields (see \autoref{apdx:INVARIANTS}).
 It is only equal to a Laplacian applied to each component in Cartesian coordinates, not on the sphere.
 …
 the $z$-coordinate therefore takes the following form:
 \begin{equation}
   \label{apdx:B_Lap_U}
+  \label{eq:DIFFOPERS_Lap_U}
+  {
     \textbf{D}}^{\textbf{U}} =
 …
 \end{align*}
 Note Bene: introducing a rotation in \autoref{apdx:B_Lap_U} does not lead to
+Note Bene: introducing a rotation in \autoref{eq:DIFFOPERS_Lap_U} does not lead to
 a useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate.
 Similarly, we did not found an expression of practical use for
 the geopotential horizontal/vertical Laplacian operator in the $s$-coordinate.
 Generally, \autoref{apdx:B_Lap_U} is used in both $z$- and $s$-coordinate systems,
+Generally, \autoref{eq:DIFFOPERS_Lap_U} is used in both $z$- and $s$-coordinate systems,
 that is a Laplacian diffusion is applied on momentum along the coordinate directions.

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex

-                      r11529
+                      r11543
 % ================================================================
 \chapter{Discrete Invariants of the Equations}
 \label{apdx:C}
+\label{apdx:INVARIANTS}
 \chaptertoc
 …
 % ================================================================
 \section{Introduction / Notations}
 \label{sec:C.0}
+\label{sec:INVARIANTS_0}
 Notation used in this appendix in the demonstations:
 …
 that is in a more compact form :
 \begin{flalign}
   \label{eq:Q2_flux}
+  \label{eq:INVARIANTS_Q2_flux}
   \partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right)
   =&                   \int_D { \frac{Q}{e_3}  \partial_t \left( e_3 \, Q \right) dv }
 …
 that is in a more compact form:
 \begin{flalign}
   \label{eq:Q2_vect}
+  \label{eq:INVARIANTS_Q2_vect}
   \partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right)
   =& \int_D {         Q   \,\partial_t Q  \;dv }
 …
 % ================================================================
 \section{Continuous conservation}
 \label{sec:C.1}
+\label{sec:INVARIANTS_1}
 The discretization of pimitive equation in $s$-coordinate (\ie\ time and space varying vertical coordinate)
 …
 The total energy (\ie\ kinetic plus potential energies) is conserved:
 \begin{flalign}
   \label{eq:Tot_Energy}
+  \label{eq:INVARIANTS_Tot_Energy}
   \partial_t \left( \int_D \left( \frac{1}{2} {\textbf{U}_h}^2 +  \rho \, g \, z \right) \;dv \right)  = & 0
 \end{flalign}
 …
 The transformation for the advection term depends on whether the vector invariant form or
 the flux form is used for the momentum equation.
 Using \autoref{eq:Q2_vect} and introducing \autoref{apdx:A_dyn_vect} in
 \autoref{eq:Tot_Energy} for the former form and
 using \autoref{eq:Q2_flux} and introducing \autoref{apdx:A_dyn_flux} in
 \autoref{eq:Tot_Energy} for the latter form leads to:
 % \label{eq:E_tot}
+Using \autoref{eq:INVARIANTS_Q2_vect} and introducing \autoref{eq:SCOORD_dyn_vect} in
+\autoref{eq:INVARIANTS_Tot_Energy} for the former form and
+using \autoref{eq:INVARIANTS_Q2_flux} and introducing \autoref{eq:SCOORD_dyn_flux} in
+\autoref{eq:INVARIANTS_Tot_Energy} for the latter form leads to:
+% \label{eq:INVARIANTS_E_tot}
 advection term (vector invariant form):
 \[
   % \label{eq:E_tot_vect_vor_1}
+  % \label{eq:INVARIANTS_E_tot_vect_vor_1}
   \int\limits_D  \zeta \; \left( \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0   \\
 \]
+%
 \[
   % \label{eq:E_tot_vect_adv_1}
+  % \label{eq:INVARIANTS_E_tot_vect_adv_1}
   \int\limits_D  \textbf{U}_h \cdot \nabla_h \left( \frac{{\textbf{U}_h}^2}{2} \right)     dv
   + \int\limits_D  \textbf{U}_h \cdot \nabla_z \textbf{U}_h  \;dv
 …
 advection term (flux form):
 \[
   % \label{eq:E_tot_flux_metric}
+  % \label{eq:INVARIANTS_E_tot_flux_metric}
   \int\limits_D  \frac{1} {e_1 e_2 } \left( v \,\partial_i e_2 - u \,\partial_j e_1  \right)\;
   \left(  \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0
 \]
 \[
   % \label{eq:E_tot_flux_adv}
+  % \label{eq:INVARIANTS_E_tot_flux_adv}
   \int\limits_D \textbf{U}_h \cdot     \left(                 {{
         \begin{array} {*{20}c}
 …
 coriolis term
 \[
   % \label{eq:E_tot_cor}
+  % \label{eq:INVARIANTS_E_tot_cor}
   \int\limits_D  f   \; \left( \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0
 \]
 pressure gradient:
 \[
   % \label{eq:E_tot_pg_1}
+  % \label{eq:INVARIANTS_E_tot_pg_1}
   - \int\limits_D  \left. \nabla p \right|_z \cdot \textbf{U}_h \;dv
   = - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
 …
 Vector invariant form:
 % \label{eq:E_tot_vect}
 \[
   % \label{eq:E_tot_vect_vor_2}
+% \label{eq:INVARIANTS_E_tot_vect}
+\[
+  % \label{eq:INVARIANTS_E_tot_vect_vor_2}
   \int\limits_D   \textbf{U}_h \cdot \text{VOR} \;dv   = 0
 \]
 \[
   % \label{eq:E_tot_vect_adv_2}
+  % \label{eq:INVARIANTS_E_tot_vect_adv_2}
   \int\limits_D  \textbf{U}_h \cdot \text{KEG}  \;dv
   + \int\limits_D  \textbf{U}_h \cdot \text{ZAD}  \;dv
 …
 \]
 \[
   % \label{eq:E_tot_pg_2}
+  % \label{eq:INVARIANTS_E_tot_pg_2}
   - \int\limits_D  \textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv
   = - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
 …
 Flux form:
 \begin{subequations}
   \label{eq:E_tot_flux}
+  \label{eq:INVARIANTS_E_tot_flux}
   \[
     % \label{eq:E_tot_flux_metric_2}
+    % \label{eq:INVARIANTS_E_tot_flux_metric_2}
     \int\limits_D  \textbf{U}_h \cdot \text {COR} \;  dv   = 0
   \]
   \[
     % \label{eq:E_tot_flux_adv_2}
+    % \label{eq:INVARIANTS_E_tot_flux_adv_2}
     \int\limits_D \textbf{U}_h \cdot \text{ADV}   \;dv
     +   \frac{1}{2} \int\limits_D {  {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t e_3  \;dv } =\;0
   \]
   \begin{equation}
     \label{eq:E_tot_pg_3}
+    \label{eq:INVARIANTS_E_tot_pg_3}
     - \int\limits_D  \textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv
     = - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
 …
 \end{subequations}
 \autoref{eq:E_tot_pg_3} is the balance between the conversion KE to PE and PE to KE.
 Indeed the left hand side of \autoref{eq:E_tot_pg_3} can be transformed as follows:
+\autoref{eq:INVARIANTS_E_tot_pg_3} is the balance between the conversion KE to PE and PE to KE.
+Indeed the left hand side of \autoref{eq:INVARIANTS_E_tot_pg_3} can be transformed as follows:
 \begin{flalign*}
   \partial_t  \left( \int\limits_D { \rho \, g \, z  \;dv} \right)
 …
 \end{flalign*}
 where the last equality is obtained by noting that the brackets is exactly the expression of $w$,
 the vertical velocity referenced to the fixe $z$-coordinate system (see \autoref{apdx:A_w_s}).
 The left hand side of \autoref{eq:E_tot_pg_3} can be transformed as follows:
+the vertical velocity referenced to the fixe $z$-coordinate system (see \autoref{eq:SCOORD_w_s}).
+The left hand side of \autoref{eq:INVARIANTS_E_tot_pg_3} can be transformed as follows:
 \begin{flalign*}
   - \int\limits_D  \left. \nabla p \right|_z & \cdot \textbf{U}_h \;dv
 …
 % ================================================================
 \section{Discrete total energy conservation: vector invariant form}
 \label{sec:C.2}
+\label{sec:INVARIANTS_2}
 % -------------------------------------------------------------------------------------------------------------
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Total energy conservation}
 \label{subsec:C_KE+PE_vect}
 The discrete form of the total energy conservation, \autoref{eq:Tot_Energy}, is given by:
+\label{subsec:INVARIANTS_KE+PE_vect}
+The discrete form of the total energy conservation, \autoref{eq:INVARIANTS_Tot_Energy}, is given by:
 \begin{flalign*}
   \partial_t  \left(  \sum\limits_{i,j,k} \biggl\{ \frac{u^2}{2} \,b_u + \frac{v^2}{2}\, b_v +  \rho \, g \, z_t \,b_t  \biggr\} \right) &=0
 …
 which in vector invariant forms, it leads to:
 \begin{equation}
   \label{eq:KE+PE_vect_discrete}
+  \label{eq:INVARIANTS_KE+PE_vect_discrete}
   \begin{split}
     \sum\limits_{i,j,k} \biggl\{   u\,                        \partial_t u         \;b_u
 …
 Substituting the discrete expression of the time derivative of the velocity either in vector invariant,
 leads to the discrete equivalent of the four equations \autoref{eq:E_tot_flux}.
+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:C_vor}
+\label{subsec:INVARIANTS_vor}
 Let $q$, located at $f$-points, be either the relative ($q=\zeta / e_{3f}$),
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsubsection{Vorticity term with ENE scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})}
 \label{subsec:C_vorENE}
+\label{subsec:INVARIANTS_vorENE}
 For the ENE scheme, the two components of the vorticity term are given by:
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsubsection{Vorticity term with EEN scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})}
 \label{subsec:C_vorEEN_vect}
+\label{subsec:INVARIANTS_vorEEN_vect}
 With the EEN scheme, the vorticity terms are represented as:
 \begin{equation}
   \tag{\ref{eq:dynvor_een}}
+  \label{eq:INVARIANTS_dynvor_een}
   \left\{ {
       \begin{aligned}
 …
 and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by:
 \begin{equation}
   \tag{\ref{eq:Q_triads}}
+  \label{eq:INVARIANTS_Q_triads}
   _i^j \mathbb{Q}^{i_p}_{j_p}
   = \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right)
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsubsection{Gradient of kinetic energy / Vertical advection}
 \label{subsec:C_zad}
+\label{subsec:INVARIANTS_zad}
 The change of Kinetic Energy (KE) due to the vertical advection is exactly balanced by the change of KE due to the horizontal gradient of KE~:
 …
+  %
   \intertext{The first term provides the discrete expression for the vertical advection of momentum (ZAD),
     while the second term corresponds exactly to \autoref{eq:KE+PE_vect_discrete}, therefore:}
+    while the second term corresponds exactly to \autoref{eq:INVARIANTS_KE+PE_vect_discrete}, therefore:}
   \equiv&                   \int\limits_D \textbf{U}_h \cdot \text{ZAD} \;dv
   + \frac{1}{2} \int_D { {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t  (e_3)  \;dv }    &&&\\
 …
 which is (over-)satified by defining the vertical scale factor as follows:
 \begin{flalign*}
   % \label{eq:e3u-e3v}
+  % \label{eq:INVARIANTS_e3u-e3v}
   e_{3u} = \frac{1}{e_{1u}\,e_{2u}}\;\overline{ e_{1t}^{ }\,e_{2t}^{ }\,e_{3t}^{ } }^{\,i+1/2}    \\
   e_{3v} = \frac{1}{e_{1v}\,e_{2v}}\;\overline{ e_{1t}^{ }\,e_{2t}^{ }\,e_{3t}^{ } }^{\,j+1/2}
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Pressure gradient term}
 \label{subsec:C.2.6}
+\label{subsec:INVARIANTS_2.6}
 \gmcomment{
 …
   \allowdisplaybreaks
   \intertext{Using successively \autoref{eq:DOM_di_adj}, \ie\ the skew symmetry property of
     the $\delta$ operator, \autoref{eq:wzv}, the continuity equation, \autoref{eq:dynhpg_sco},
+    the $\delta$ operator, \autoref{eq:DYN_wzv}, the continuity equation, \autoref{eq:DYN_hpg_sco},
     the hydrostatic equation in the $s$-coordinate, and $\delta_{k+1/2} \left[ z_t \right] \equiv e_{3w} $,
     which comes from the definition of $z_t$, it becomes: }
 …
+  %
 \end{flalign*}
 The first term is exactly the first term of the right-hand-side of \autoref{eq:KE+PE_vect_discrete}.
+The first term is exactly the first term of the right-hand-side of \autoref{eq:INVARIANTS_KE+PE_vect_discrete}.
 It remains to demonstrate that the last term,
 which is obviously a discrete analogue of $\int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to
 the last term of \autoref{eq:KE+PE_vect_discrete}.
+the last term of \autoref{eq:INVARIANTS_KE+PE_vect_discrete}.
 In other words, the following property must be satisfied:
 \begin{flalign*}
 …
 % ================================================================
 \section{Discrete total energy conservation: flux form}
 \label{sec:C.3}
+\label{sec:INVARIANTS_3}
 % -------------------------------------------------------------------------------------------------------------
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Total energy conservation}
 \label{subsec:C_KE+PE_flux}
 The discrete form of the total energy conservation, \autoref{eq:Tot_Energy}, is given by:
+\label{subsec:INVARIANTS_KE+PE_flux}
+The discrete form of the total energy conservation, \autoref{eq:INVARIANTS_Tot_Energy}, is given by:
 \begin{flalign*}
   \partial_t \left(  \sum\limits_{i,j,k} \biggl\{ \frac{u^2}{2} \,b_u + \frac{v^2}{2}\, b_v +  \rho \, g \, z_t \,b_t  \biggr\} \right) &=0  \\
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Coriolis and advection terms: flux form}
 \label{subsec:C.3.2}
+\label{subsec:INVARIANTS_3.2}
 % -------------------------------------------------------------------------------------------------------------
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsubsection{Coriolis plus ``metric'' term}
 \label{subsec:C.3.3}
+\label{subsec:INVARIANTS_3.3}
 In flux from the vorticity term reduces to a Coriolis term in which
 …
 Either the ENE or EEN scheme is then applied to obtain the vorticity term in flux form.
 It therefore conserves the total KE.
 The derivation is the same as for the vorticity term in the vector invariant form (\autoref{subsec:C_vor}).
+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:C.3.4}
+\label{subsec:INVARIANTS_3.4}
 The flux form operator of the momentum advection is evaluated using
 …
 \begin{equation}
   \label{eq:C_ADV_KE_flux}
+  \label{eq:INVARIANTS_ADV_KE_flux}
   -  \int_D \textbf{U}_h \cdot     \left(                 {{
         \begin{array} {*{20}c}
 …
 \]
 which is the discrete form of $ \frac{1}{2} \int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv $.
 \autoref{eq:C_ADV_KE_flux} is thus satisfied.
+\autoref{eq:INVARIANTS_ADV_KE_flux} is thus satisfied.
 When the UBS scheme is used to evaluate the flux form momentum advection,
 …
 % ================================================================
 \section{Discrete enstrophy conservation}
 \label{sec:C.4}
+\label{sec:INVARIANTS_4}
 % -------------------------------------------------------------------------------------------------------------
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsubsection{Vorticity term with ENS scheme  (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})}
 \label{subsec:C_vorENS}
+\label{subsec:INVARIANTS_vorENS}
 In the ENS scheme, the vorticity term is descretized as follows:
 \begin{equation}
   \tag{\ref{eq:dynvor_ens}}
+  \label{eq:INVARIANTS_dynvor_ens}
   \left\{
     \begin{aligned}
 …
 it can be shown that:
 \begin{equation}
   \label{eq:C_1.1}
+  \label{eq:INVARIANTS_1.1}
   \int_D {q\,\;{\textbf{k}}\cdot \frac{1} {e_3} \nabla \times \left( {e_3 \, q \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} \equiv 0
 \end{equation}
 where $dv=e_1\,e_2\,e_3 \; di\,dj\,dk$ is the volume element.
 Indeed, using \autoref{eq:dynvor_ens},
 the discrete form of the right hand side of \autoref{eq:C_1.1} can be transformed as follow:
+Indeed, using \autoref{eq:DYN_vor_ens},
+the discrete form of the right hand side of \autoref{eq:INVARIANTS_1.1} can be transformed as follow:
 \begin{flalign*}
   &\int_D q \,\; \textbf{k} \cdot \frac{1} {e_3 } \nabla \times
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsubsection{Vorticity Term with EEN scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})}
 \label{subsec:C_vorEEN}
+\label{subsec:INVARIANTS_vorEEN}
 With the EEN scheme, the vorticity terms are represented as:
 \begin{equation}
   \tag{\ref{eq:dynvor_een}}
+  \label{eq:INVARIANTS_dynvor_een}
   \left\{ {
       \begin{aligned}
 …
 and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by:
 \begin{equation}
   \tag{\ref{eq:Q_triads}}
+  \tag{\ref{eq:INVARIANTS_Q_triads}}
   _i^j \mathbb{Q}^{i_p}_{j_p}
   = \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right)
 …
 Let consider one of the vorticity triad, for example ${^{i}_j}\mathbb{Q}^{+1/2}_{+1/2} $,
 similar manipulation can be done for the 3 others.
 The discrete form of the right hand side of \autoref{eq:C_1.1} applied to
+The discrete form of the right hand side of \autoref{eq:INVARIANTS_1.1} applied to
 this triad only can be transformed as follow:
 …
 % ================================================================
 \section{Conservation properties on tracers}
 \label{sec:C.5}
+\label{sec:INVARIANTS_5}
 All the numerical schemes used in \NEMO\ are written such that the tracer content is conserved by
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Advection term}
 \label{subsec:C.5.1}
+\label{subsec:INVARIANTS_5.1}
 conservation of a tracer, $T$:
 …
 % ================================================================
 \section{Conservation properties on lateral momentum physics}
 \label{sec:dynldf_properties}
+\label{sec:INVARIANTS_dynldf_properties}
 The discrete formulation of the horizontal diffusion of momentum ensures
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Conservation of potential vorticity}
 \label{subsec:C.6.1}
+\label{subsec:INVARIANTS_6.1}
 The lateral momentum diffusion term conserves the potential vorticity:
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Dissipation of horizontal kinetic energy}
 \label{subsec:C.6.2}
+\label{subsec:INVARIANTS_6.2}
 The lateral momentum diffusion term dissipates the horizontal kinetic energy:
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Dissipation of enstrophy}
 \label{subsec:C.6.3}
+\label{subsec:INVARIANTS_6.3}
 The lateral momentum diffusion term dissipates the enstrophy when the eddy coefficients are horizontally uniform:
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Conservation of horizontal divergence}
 \label{subsec:C.6.4}
+\label{subsec:INVARIANTS_6.4}
 When the horizontal divergence of the horizontal diffusion of momentum (discrete sense) is taken,
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Dissipation of horizontal divergence variance}
 \label{subsec:C.6.5}
+\label{subsec:INVARIANTS_6.5}
 \begin{flalign*}
 …
 % ================================================================
 \section{Conservation properties on vertical momentum physics}
 \label{sec:C.7}
+\label{sec:INVARIANTS_7}
 As for the lateral momentum physics,
 …
 % ================================================================
 \section{Conservation properties on tracer physics}
 \label{sec:C.8}
+\label{sec:INVARIANTS_8}
 The numerical schemes used for tracer subgridscale physics are written such that
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Conservation of tracers}
 \label{subsec:C.8.1}
+\label{subsec:INVARIANTS_8.1}
 constraint of conservation of tracers:
 …
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Dissipation of tracer variance}
 \label{subsec:C.8.2}
+\label{subsec:INVARIANTS_8.2}
 constraint on the dissipation of tracer variance:

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_s_coord.tex

-                      r11529
+                      r11543
 % ================================================================
 \chapter{Curvilinear $s-$Coordinate Equations}
 \label{apdx:A}
+\label{apdx:SCOORD}
 \chaptertoc
 …
 % ================================================================
 \section{Chain rule for $s-$coordinates}
 \label{sec:A_chain}
+\label{sec:SCOORD_chain}
 In order to establish the set of Primitive Equation in curvilinear $s$-coordinates
 (\ie\ an orthogonal curvilinear coordinate in the horizontal and
 an Arbitrary Lagrangian Eulerian (ALE) coordinate in the vertical),
 we start from the set of equations established in \autoref{subsec:PE_zco_Eq} for
+we start from the set of equations established in \autoref{subsec:MB_zco_Eq} for
 the special case $k = z$ and thus $e_3 = 1$,
 and we introduce an arbitrary vertical coordinate $a = a(i,j,z,t)$.
 …
 the horizontal slope of $s-$surfaces by:
 \begin{equation}
   \label{apdx:A_s_slope}
+  \label{eq:SCOORD_s_slope}
   \sigma_1 =\frac{1}{e_1 } \; \left. {\frac{\partial z}{\partial i}} \right|_s
   \quad \text{and} \quad
 …
 The model fields (e.g. pressure $p$) can be viewed as functions of $(i,j,z,t)$ (e.g. $p(i,j,z,t)$) or as
 functions of $(i,j,s,t)$ (e.g. $p(i,j,s,t)$). The symbol $\bullet$ will be used to represent any one of
 these fields.  Any ``infinitesimal'' change in $\bullet$ can be written in two forms:
 \begin{equation}
   \label{apdx:A_s_infin_changes}
+functions of $(i,j,s,t)$ (e.g. $p(i,j,s,t)$). The symbol $\bullet$ will be used to represent any one of
+these fields.  Any ``infinitesimal'' change in $\bullet$ can be written in two forms:
+\begin{equation}
+  \label{eq:SCOORD_s_infin_changes}
   \begin{aligned}
     & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,s,t}
                 + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,s,t}
                 + \delta s \left. \frac{ \partial \bullet }{\partial s} \right|_{i,j,t}
+    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,s,t}
+                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,s,t}
+                + \delta s \left. \frac{ \partial \bullet }{\partial s} \right|_{i,j,t}
                 + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,s} , \\
     & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,z,t}
                 + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,z,t}
                 + \delta z \left. \frac{ \partial \bullet }{\partial z} \right|_{i,j,t}
+    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,z,t}
+                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,z,t}
+                + \delta z \left. \frac{ \partial \bullet }{\partial z} \right|_{i,j,t}
                 + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,z} .
   \end{aligned}
 …
 Using the first form and considering a change $\delta i$ with $j, z$ and $t$ held constant, shows that
 \begin{equation}
   \label{apdx:A_s_chain_rule}
+  \label{eq:SCOORD_s_chain_rule}
       \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,z,t}  =
       \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,s,t}
     + \left. {\frac{\partial s       }{\partial i}} \right|_{j,z,t} \;
       \left. {\frac{\partial \bullet }{\partial s}} \right|_{i,j,t} .
 \end{equation}
 The term $\left. \partial s / \partial i \right|_{j,z,t}$ can be related to the slope of constant $s$ surfaces,
 (\autoref{apdx:A_s_slope}), by applying the second of (\autoref{apdx:A_s_infin_changes}) with $\bullet$ set to
+    + \left. {\frac{\partial s       }{\partial i}} \right|_{j,z,t} \;
+      \left. {\frac{\partial \bullet }{\partial s}} \right|_{i,j,t} .
+\end{equation}
+The term $\left. \partial s / \partial i \right|_{j,z,t}$ can be related to the slope of constant $s$ surfaces,
+(\autoref{eq:SCOORD_s_slope}), by applying the second of (\autoref{eq:SCOORD_s_infin_changes}) with $\bullet$ set to
 $s$ and $j, t$ held constant
 \begin{equation}
 \label{apdx:a_delta_s}
 \delta s|_{j,t} =
          \delta i \left. \frac{ \partial s }{\partial i} \right|_{j,z,t}
+\label{eq:SCOORD_delta_s}
+\delta s|_{j,t} =
+         \delta i \left. \frac{ \partial s }{\partial i} \right|_{j,z,t}
        + \delta z \left. \frac{ \partial s }{\partial z} \right|_{i,j,t} .
 \end{equation}
 Choosing to look at a direction in the $(i,z)$ plane in which $\delta s = 0$ and using
 (\autoref{apdx:A_s_slope}) we obtain
 \begin{equation}
 \left. \frac{ \partial s }{\partial i} \right|_{j,z,t} =
+(\autoref{eq:SCOORD_s_slope}) we obtain
+\begin{equation}
+\left. \frac{ \partial s }{\partial i} \right|_{j,z,t} =
          -  \left. \frac{ \partial z }{\partial i} \right|_{j,s,t} \;
             \left. \frac{ \partial s }{\partial z} \right|_{i,j,t}
     = - \frac{e_1 }{e_3 }\sigma_1  .
 \label{apdx:a_ds_di_z}
 \end{equation}
 Another identity, similar in form to (\autoref{apdx:a_ds_di_z}), can be derived
 by choosing $\bullet$ to be $s$ and using the second form of (\autoref{apdx:A_s_infin_changes}) to consider
+\label{eq:SCOORD_ds_di_z}
+\end{equation}
+Another identity, similar in form to (\autoref{eq:SCOORD_ds_di_z}), can be derived
+by choosing $\bullet$ to be $s$ and using the second form of (\autoref{eq:SCOORD_s_infin_changes}) to consider
 changes in which $i , j$ and $s$ are constant. This shows that
 \begin{equation}
 \label{apdx:A_w_in_s}
 w_s = \left. \frac{ \partial z }{\partial t} \right|_{i,j,s} =
+\label{eq:SCOORD_w_in_s}
+w_s = \left. \frac{ \partial z }{\partial t} \right|_{i,j,s} =
 - \left. \frac{ \partial z }{\partial s} \right|_{i,j,t}
   \left. \frac{ \partial s }{\partial t} \right|_{i,j,z}
   = - e_3 \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} .
 \end{equation}
 In what follows, for brevity, indication of the constancy of the $i, j$ and $t$ indices is
 usually omitted. Using the arguments outlined above one can show that the chain rules needed to establish
+  \left. \frac{ \partial s }{\partial t} \right|_{i,j,z}
+  = - e_3 \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} .
+\end{equation}
+In what follows, for brevity, indication of the constancy of the $i, j$ and $t$ indices is
+usually omitted. Using the arguments outlined above one can show that the chain rules needed to establish
 the model equations in the curvilinear $s-$coordinate system are:
 \begin{equation}
   \label{apdx:A_s_chain_rule}
+  \label{eq:SCOORD_s_chain_rule}
   \begin{aligned}
     &\left. {\frac{\partial \bullet }{\partial t}} \right|_z  =
     \left. {\frac{\partial \bullet }{\partial t}} \right|_s
+    \left. {\frac{\partial \bullet }{\partial t}} \right|_s
     + \frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial t} , \\
     &\left. {\frac{\partial \bullet }{\partial i}} \right|_z  =
     \left. {\frac{\partial \bullet }{\partial i}} \right|_s
     +\frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial i}=
     \left. {\frac{\partial \bullet }{\partial i}} \right|_s
+    \left. {\frac{\partial \bullet }{\partial i}} \right|_s
     -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} , \\
     &\left. {\frac{\partial \bullet }{\partial j}} \right|_z  =
     \left. {\frac{\partial \bullet }{\partial j}} \right|_s
+    \left. {\frac{\partial \bullet }{\partial j}} \right|_s
     + \frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}=
     \left. {\frac{\partial \bullet }{\partial j}} \right|_s
+    \left. {\frac{\partial \bullet }{\partial j}} \right|_s
     - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} , \\
     &\;\frac{\partial \bullet }{\partial z}  \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s} .
 …
 % ================================================================
 \section{Continuity equation in $s-$coordinates}
 \label{sec:A_continuity}
 Using (\autoref{apdx:A_s_chain_rule}) and
+\label{sec:SCOORD_continuity}
+Using (\autoref{eq:SCOORD_s_chain_rule}) and
 the fact that the horizontal scale factors $e_1$ and $e_2$ do not depend on the vertical coordinate,
 the divergence of the velocity relative to the ($i$,$j$,$z$) coordinate system is transformed as follows in order to
 …
 \end{subequations}
 Here, $w$ is the vertical velocity relative to the $z-$coordinate system.
 Using the first form of (\autoref{apdx:A_s_infin_changes})
 and the definitions (\autoref{apdx:A_s_slope}) and (\autoref{apdx:A_w_in_s}) for $\sigma_1$, $\sigma_2$ and  $w_s$,
+Here, $w$ is the vertical velocity relative to the $z-$coordinate system.
+Using the first form of (\autoref{eq:SCOORD_s_infin_changes})
+and the definitions (\autoref{eq:SCOORD_s_slope}) and (\autoref{eq:SCOORD_w_in_s}) for $\sigma_1$, $\sigma_2$ and  $w_s$,
 one can show that the vertical velocity, $w_p$ of a point
 moving with the horizontal velocity of the fluid along an $s$ surface is given by
 \begin{equation}
 \label{apdx:A_w_p}
+moving with the horizontal velocity of the fluid along an $s$ surface is given by
+\begin{equation}
+\label{eq:SCOORD_w_p}
 \begin{split}
 w_p  = & \left. \frac{ \partial z }{\partial t} \right|_s
      + \frac{u}{e_1} \left. \frac{ \partial z }{\partial i} \right|_s
+     + \frac{u}{e_1} \left. \frac{ \partial z }{\partial i} \right|_s
      + \frac{v}{e_2} \left. \frac{ \partial z }{\partial j} \right|_s \\
      = & w_s + u \sigma_1 + v \sigma_2 .
 \end{split}
+\end{split}
 \end{equation}
  The vertical velocity across this surface is denoted by
 \begin{equation}
   \label{apdx:A_w_s}
   \omega  = w - w_p = w - ( w_s + \sigma_1 \,u + \sigma_2 \,v )  .
 \end{equation}
 Hence
 \begin{equation}
 \frac{1}{e_3 } \frac{\partial}{\partial s}   \left[  w -  u\;\sigma_1  - v\;\sigma_2  \right] =
 \frac{1}{e_3 } \frac{\partial}{\partial s} \left[  \omega + w_s \right] =
    \frac{1}{e_3 } \left[ \frac{\partial \omega}{\partial s}
  + \left. \frac{ \partial }{\partial t} \right|_s \frac{\partial z}{\partial s} \right] =
    \frac{1}{e_3 } \frac{\partial \omega}{\partial s} + \frac{1}{e_3 } \left. \frac{ \partial e_3}{\partial t} . \right|_s
 \end{equation}
 Using (\autoref{apdx:A_w_s}) in our expression for $\nabla \cdot {\mathrm {\mathbf U}}$ we obtain
+  \label{eq:SCOORD_w_s}
+  \omega  = w - w_p = w - ( w_s + \sigma_1 \,u + \sigma_2 \,v )  .
+\end{equation}
+Hence
+\begin{equation}
+\frac{1}{e_3 } \frac{\partial}{\partial s}   \left[  w -  u\;\sigma_1  - v\;\sigma_2  \right] =
+\frac{1}{e_3 } \frac{\partial}{\partial s} \left[  \omega + w_s \right] =
+   \frac{1}{e_3 } \left[ \frac{\partial \omega}{\partial s}
+ + \left. \frac{ \partial }{\partial t} \right|_s \frac{\partial z}{\partial s} \right] =
+   \frac{1}{e_3 } \frac{\partial \omega}{\partial s} + \frac{1}{e_3 } \left. \frac{ \partial e_3}{\partial t} . \right|_s
+\end{equation}
+Using (\autoref{eq:SCOORD_w_s}) in our expression for $\nabla \cdot {\mathrm {\mathbf U}}$ we obtain
 our final expression for the divergence of the velocity in the curvilinear $s-$coordinate system:
 \begin{equation}
 …
 \end{equation}
 As a result, the continuity equation \autoref{eq:PE_continuity} in the $s-$coordinates is:
 \begin{equation}
   \label{apdx:A_sco_Continuity}
+As a result, the continuity equation \autoref{eq:MB_PE_continuity} in the $s-$coordinates is:
+\begin{equation}
+  \label{eq:SCOORD_sco_Continuity}
   \frac{1}{e_3 } \frac{\partial e_3}{\partial t}
   + \frac{1}{e_1 \,e_2 \,e_3 }\left[
 …
 % ================================================================
 \section{Momentum equation in $s-$coordinate}
 \label{sec:A_momentum}
+\label{sec:SCOORD_momentum}
 Here we only consider the first component of the momentum equation,
 …
 $\bullet$ \textbf{Total derivative in vector invariant form}
 Let us consider \autoref{eq:PE_dyn_vect}, the first component of the momentum equation in the vector invariant form.
+Let us consider \autoref{eq:MB_dyn_vect}, the first component of the momentum equation in the vector invariant form.
 Its total $z-$coordinate time derivative,
 $\left. \frac{D u}{D t} \right|_z$ can be transformed as follows in order to obtain
 …
         +  w \;\frac{\partial u}{\partial z}      \\
+        %
       \intertext{introducing the chain rule (\autoref{apdx:A_s_chain_rule}) }
+      \intertext{introducing the chain rule (\autoref{eq:SCOORD_s_chain_rule}) }
+      %
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
 …
         \; \frac{\partial u}{\partial s} .  \\
+        %
       \intertext{Introducing $\omega$, the dia-s-surface velocity given by (\autoref{apdx:A_w_s}) }
+      \intertext{Introducing $\omega$, the dia-s-surface velocity given by (\autoref{eq:SCOORD_w_s}) }
+      %
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
 …
 \end{subequations}
+%
 Applying the time derivative chain rule (first equation of (\autoref{apdx:A_s_chain_rule})) to $u$ and
 using (\autoref{apdx:A_w_in_s}) provides the expression of the last term of the right hand side,
+Applying the time derivative chain rule (first equation of (\autoref{eq:SCOORD_s_chain_rule})) to $u$ and
+using (\autoref{eq:SCOORD_w_in_s}) provides the expression of the last term of the right hand side,
 \[
+  {
 …
 \ie\ the total $s-$coordinate time derivative :
 \begin{align}
   \label{apdx:A_sco_Dt_vect}
+  \label{eq:SCOORD_sco_Dt_vect}
   \left. \frac{D u}{D t} \right|_s
   = \left. {\frac{\partial u }{\partial t}} \right|_s
   - \left. \zeta \right|_s \;v
   + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s
   + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} .
+  + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} .
 \end{align}
 Therefore, the vector invariant form of the total time derivative has exactly the same mathematical form in
 …
 Let us start from the total time derivative in the curvilinear $s-$coordinate system we have just establish.
 Following the procedure used to establish (\autoref{eq:PE_flux_form}), it can be transformed into :
+Following the procedure used to establish (\autoref{eq:MB_flux_form}), it can be transformed into :
 % \begin{subequations}
 \begin{align*}
 …
 \end{align*}
+%
 Introducing the vertical scale factor inside the horizontal derivative of the first two terms
+Introducing the vertical scale factor inside the horizontal derivative of the first two terms
 (\ie\ the horizontal divergence), it becomes :
 \begin{align*}
 …
   \begin{array}{*{20}l}
     % \begin{align*} {\begin{array}{*{20}l}
     %     {\begin{array}{*{20}l} \left. \frac{D u}{D t} \right|_s
+    %     {\begin{array}{*{20}l} \left. \frac{D u}{D t} \right|_s
     &= \left. {\frac{\partial u }{\partial t}} \right|_s
     &+ \frac{1}{e_1\,e_2\,e_3}  \left(  \frac{\partial( e_2 e_3 \,u^2 )}{\partial i}
 …
+     %
     \intertext {Introducing a more compact form for the divergence of the momentum fluxes,
     and using (\autoref{apdx:A_sco_Continuity}), the $s-$coordinate continuity equation,
+    and using (\autoref{eq:SCOORD_sco_Continuity}), the $s-$coordinate continuity equation,
     it becomes : }
+  %
 …
+  }
 \end{align*}
 which leads to the $s-$coordinate flux formulation of the total $s-$coordinate time derivative,
+which leads to the $s-$coordinate flux formulation of the total $s-$coordinate time derivative,
 \ie\ the total $s-$coordinate time derivative in flux form:
 \begin{flalign}
   \label{apdx:A_sco_Dt_flux}
+  \label{eq:SCOORD_sco_Dt_flux}
   \left. \frac{D u}{D t} \right|_s   = \frac{1}{e_3}  \left. \frac{\partial ( e_3\,u)}{\partial t} \right|_s
   + \left.  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)    \right|_s
 …
 It has the same form as in the $z-$coordinate but for
 the vertical scale factor that has appeared inside the time derivative which
 comes from the modification of (\autoref{apdx:A_sco_Continuity}),
+comes from the modification of (\autoref{eq:SCOORD_sco_Continuity}),
 the continuity equation.
 …
 \]
 Applying similar manipulation to the second component and
 replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it becomes:
 \begin{equation}
   \label{apdx:A_grad_p_1}
+replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{eq:SCOORD_s_slope}, it becomes:
+\begin{equation}
+  \label{eq:SCOORD_grad_p_1}
   \begin{split}
     -\frac{1}{\rho_o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z
 …
 \end{equation}
 An additional term appears in (\autoref{apdx:A_grad_p_1}) which accounts for
+An additional term appears in (\autoref{eq:SCOORD_grad_p_1}) which accounts for
 the tilt of $s-$surfaces with respect to geopotential $z-$surfaces.
 …
 Therefore, $p$ and $p_h'$ are linked through:
 \begin{equation}
   \label{apdx:A_pressure}
+  \label{eq:SCOORD_pressure}
   p = \rho_o \; p_h' + \rho_o \, g \, ( \eta - z )
 \end{equation}
 …
 \]
 Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and
+Substituing \autoref{eq:SCOORD_pressure} in \autoref{eq:SCOORD_grad_p_1} and
 using the definition of the density anomaly it becomes an expression in two parts:
 \begin{equation}
   \label{apdx:A_grad_p_2}
+  \label{eq:SCOORD_grad_p_2}
   \begin{split}
     -\frac{1}{\rho_o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z
 …
 This formulation of the pressure gradient is characterised by the appearance of
 a term depending on the sea surface height only
 (last term on the right hand side of expression \autoref{apdx:A_grad_p_2}).
+(last term on the right hand side of expression \autoref{eq:SCOORD_grad_p_2}).
 This term will be loosely termed \textit{surface pressure gradient} whereas
 the first term will be termed the \textit{hydrostatic pressure gradient} by analogy to
 …
 The coriolis and forcing terms as well as the the vertical physics remain unchanged as
 they involve neither time nor space derivatives.
 The form of the lateral physics is discussed in \autoref{apdx:B}.
+The form of the lateral physics is discussed in \autoref{apdx:DIFFOPERS}.
 $\bullet$ \textbf{Full momentum equation}
 …
 the one in a curvilinear $z-$coordinate, except for the pressure gradient term:
 \begin{subequations}
   \label{apdx:A_dyn_vect}
+  \label{eq:SCOORD_dyn_vect}
   \begin{multline}
     \label{apdx:A_PE_dyn_vect_u}
+    \label{eq:SCOORD_PE_dyn_vect_u}
     \frac{\partial u}{\partial t}=
     +   \left( {\zeta +f} \right)\,v
 …
   \end{multline}
   \begin{multline}
     \label{apdx:A_dyn_vect_v}
+    \label{eq:SCOORD_dyn_vect_v}
     \frac{\partial v}{\partial t}=
     -   \left( {\zeta +f} \right)\,u
 …
 the formulation of both the time derivative and the pressure gradient term:
 \begin{subequations}
   \label{apdx:A_dyn_flux}
+  \label{eq:SCOORD_dyn_flux}
   \begin{multline}
     \label{apdx:A_PE_dyn_flux_u}
+    \label{eq:SCOORD_PE_dyn_flux_u}
     \frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t} =
     - \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)
 …
   \end{multline}
   \begin{multline}
     \label{apdx:A_dyn_flux_v}
+    \label{eq:SCOORD_dyn_flux_v}
     \frac{1}{e_3}\frac{\partial \left(  e_3\,v  \right) }{\partial t}=
     -  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,v}   \right)
 …
     -   \frac{1}{e_2 } \left(    \frac{\partial p_h'}{\partial j} + g\; d  \; \frac{\partial z}{\partial j}    \right)
     -   \frac{g}{e_2 } \frac{\partial \eta}{\partial j}
     +  D_v^{\vect{U}}  +   F_v^{\vect{U}} .
+    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} .
   \end{multline}
 \end{subequations}
 …
 hydrostatic pressure and density anomalies, $p_h'$ and $d=( \frac{\rho}{\rho_o}-1 )$:
 \begin{equation}
   \label{apdx:A_dyn_zph}
+  \label{eq:SCOORD_dyn_zph}
   \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 .
 \end{equation}
 …
 in particular the pressure gradient.
 By contrast, $\omega$ is not $w$, the third component of the velocity, but the dia-surface velocity component,
 \ie\ the volume flux across the moving $s$-surfaces per unit horizontal area.
+\ie\ the volume flux across the moving $s$-surfaces per unit horizontal area.
 …
 % ================================================================
 \section{Tracer equation}
 \label{sec:A_tracer}
+\label{sec:SCOORD_tracer}
 The tracer equation is obtained using the same calculation as for the continuity equation and then
 …
 \begin{multline}
   \label{apdx:A_tracer}
+  \label{eq:SCOORD_tracer}
   \frac{1}{e_3} \frac{\partial \left(  e_3 T  \right)}{\partial t}
   = -\frac{1}{e_1 \,e_2 \,e_3}
 …
 \end{multline}
 The expression for the advection term is a straight consequence of (\autoref{apdx:A_sco_Continuity}),
 the expression of the 3D divergence in the $s-$coordinates established above.
+The expression for the advection term is a straight consequence of (\autoref{eq:SCOORD_sco_Continuity}),
+the expression of the 3D divergence in the $s-$coordinates established above.
 \biblio

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex

-                      r11529
+                      r11543
 % ================================================================
 \chapter{Iso-Neutral Diffusion and Eddy Advection using Triads}
 \label{apdx:triad}
+\label{apdx:TRIADS}
 \chaptertoc
 …
 \begin{description}
 \item[\np{ln\_triad\_iso}]
   See \autoref{sec:taper}.
+  See \autoref{sec:TRIADS_taper}.
   If this is set false (the default),
   then `iso-neutral' mixing is accomplished within the surface mixed-layer along slopes linearly decreasing with
   depth from the value immediately below the mixed-layer to zero (flat) at the surface (\autoref{sec:lintaper}).
+  depth from the value immediately below the mixed-layer to zero (flat) at the surface (\autoref{sec:TRIADS_lintaper}).
   This is the same treatment as used in the default implementation
   \autoref{subsec:LDF_slp_iso}; \autoref{fig:eiv_slp}.
+  \autoref{subsec:LDF_slp_iso}; \autoref{fig:LDF_eiv_slp}.
   Where \np{ln\_triad\_iso} is set true,
   the vertical skew flux is further reduced to ensure no vertical buoyancy flux,
   giving an almost pure horizontal diffusive tracer flux within the mixed layer.
   This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:Gerdes-taper}
+  This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:TRIADS_Gerdes-taper}
 \item[\np{ln\_botmix\_triad}]
   See \autoref{sec:iso_bdry}.
+  See \autoref{sec:TRIADS_iso_bdry}.
   If this is set false (the default) then the lateral diffusive fluxes
   associated with triads partly masked by topography are neglected.
   If it is set true, however, then these lateral diffusive fluxes are applied,
+  associated with triads partly masked by topography are neglected.
+  If it is set true, however, then these lateral diffusive fluxes are applied,
   giving smoother bottom tracer fields at the cost of introducing diapycnal mixing.
 \item[\np{rn\_sw\_triad}]
 …
 \section{Triad formulation of iso-neutral diffusion}
 \label{sec:iso}
+\label{sec:TRIADS_iso}
 We have implemented into \NEMO\ a scheme inspired by \citet{griffies.gnanadesikan.ea_JPO98},
 …
 The iso-neutral second order tracer diffusive operator for small angles between
 iso-neutral surfaces and geopotentials is given by \autoref{eq:iso_tensor_1}:
+iso-neutral surfaces and geopotentials is given by \autoref{eq:TRIADS_iso_tensor_1}:
 \begin{subequations}
   \label{eq:iso_tensor_1}
+  \label{eq:TRIADS_iso_tensor_1}
   \begin{equation}
     D^{lT}=-\nabla \cdot\vect{f}^{lT}\equiv
 …
   \end{equation}
   \begin{equation}
     \label{eq:iso_tensor_2}
+    \label{eq:TRIADS_iso_tensor_2}
     \mbox{with}\quad \;\;\Re =
     \begin{pmatrix}
 …
 %  {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\
 % \end{array} }} \right)
 Here \autoref{eq:PE_iso_slopes}
+Here \autoref{eq:MB_iso_slopes}
 \begin{align*}
   r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i}
 …
 We will find it useful to consider the fluxes per unit area in $i,j,k$ space; we write
 \[
   % \label{eq:Fijk}
+  % \label{eq:TRIADS_Fijk}
   \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right).
 \]
 …
 The off-diagonal terms of the small angle diffusion tensor
 \autoref{eq:iso_tensor_1}, \autoref{eq:iso_tensor_2} produce skew-fluxes along
+\autoref{eq:TRIADS_iso_tensor_1}, \autoref{eq:TRIADS_iso_tensor_2} produce skew-fluxes along
 the $i$- and $j$-directions resulting from the vertical tracer gradient:
 \begin{align}
   \label{eq:i13c}
+  \label{eq:TRIADS_i13c}
   f_{13}=&+{A^{lT}} r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad f_{23}=+{A^{lT}} r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\
   \intertext{and in the k-direction resulting from the lateral tracer gradients}
   \label{eq:i31c}
+  \label{eq:TRIADS_i31c}
   f_{31}+f_{32}=& {A^{lT}} r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+{A^{lT}} r_2\frac{1}{e_1}\frac{\partial T}{\partial i}
 \end{align}
 …
 The vertical diffusive flux associated with the $_{33}$ component of the small angle diffusion tensor is
 \begin{equation}
   \label{eq:i33c}
+  \label{eq:TRIADS_i33c}
   f_{33}=-{A^{lT}}(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}.
 \end{equation}
 …
 The following description will describe the fluxes on the $i$-$k$ plane.
 There is no natural discretization for the $i$-component of the skew-flux, \autoref{eq:i13c},
+There is no natural discretization for the $i$-component of the skew-flux, \autoref{eq:TRIADS_i13c},
 as although it must be evaluated at $u$-points,
 it involves vertical gradients (both for the tracer and the slope $r_1$), defined at $w$-points.
 Similarly, the vertical skew flux, \autoref{eq:i31c},
+Similarly, the vertical skew flux, \autoref{eq:TRIADS_i31c},
 is evaluated at $w$-points but involves horizontal gradients defined at $u$-points.
 …
 The straightforward approach to discretize the lateral skew flux
 \autoref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA,
 \autoref{eq:tra_ldf_iso}, is to calculate a mean vertical gradient at the $u$-point from
+\autoref{eq:TRIADS_i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA,
+\autoref{eq:TRA_ldf_iso}, is to calculate a mean vertical gradient at the $u$-point from
 the average of the four surrounding vertical tracer gradients, and multiply this by a mean slope at the $u$-point,
 calculated from the averaged surrounding vertical density gradients.
 The total area-integrated skew-flux (flux per unit area in $ijk$ space) from tracer cell $i,k$ to $i+1,k$,
 noting that the $e_{{3}_{i+1/2}^k}$ in the area $e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with
 the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer gradient, is then \autoref{eq:tra_ldf_iso}
+the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer gradient, is then \autoref{eq:TRA_ldf_iso}
 \[
   \left(F_u^{13} \right)_{i+\frac{1}{2}}^k = {A}_{i+\frac{1}{2}}^k
 …
     \includegraphics[width=\textwidth]{Fig_GRIFF_triad_fluxes}
     \caption{
       \protect\label{fig:ISO_triad}
+      \protect\label{fig:TRIADS_ISO_triad}
       (a) Arrangement of triads $S_i$ and tracer gradients to
       give lateral tracer flux from box $i,k$ to $i+1,k$
 …
 the corresponding `triad' slope calculated from the lateral density gradient across the $u$-point divided by
 the vertical density gradient at the same $w$-point as the tracer gradient.
 See \autoref{fig:ISO_triad}a, where the thick lines denote the tracer gradients,
+See \autoref{fig:TRIADS_ISO_triad}a, where the thick lines denote the tracer gradients,
 and the thin lines the corresponding triads, with slopes $s_1, \dotsc s_4$.
 The total area-integrated skew-flux from tracer cell $i,k$ to $i+1,k$
 \begin{multline}
   \label{eq:i13}
+  \label{eq:TRIADS_i13}
   \left( F_u^{13}  \right)_{i+\frac{1}{2}}^k = {A}_{i+1}^k a_1 s_1
   \delta_{k+\frac{1}{2}} \left[ T^{i+1}
 …
 This discretization gives a much closer stencil, and disallows the two-point computational modes.
 The vertical skew flux \autoref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at
 the $w$-point $i,k+\frac{1}{2}$ is constructed similarly (\autoref{fig:ISO_triad}b) by
+The vertical skew flux \autoref{eq:TRIADS_i31c} from tracer cell $i,k$ to $i,k+1$ at
+the $w$-point $i,k+\frac{1}{2}$ is constructed similarly (\autoref{fig:TRIADS_ISO_triad}b) by
 multiplying lateral tracer gradients from each of the four surrounding $u$-points by the appropriate triad slope:
 \begin{multline}
   \label{eq:i31}
+  \label{eq:TRIADS_i31}
   \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} =  {A}_i^{k+1} a_{1}'
   s_{1}' \delta_{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1}
 …
 (appearing in both the vertical and lateral gradient),
 and the $u$- and $w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the triad as follows
 (see also \autoref{fig:ISO_triad}):
 \begin{equation}
   \label{eq:R}
+(see also \autoref{fig:TRIADS_ISO_triad}):
+\begin{equation}
+  \label{eq:TRIADS_R}
   _i^k \mathbb{R}_{i_p}^{k_p}
   =-\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}}
 …
     \includegraphics[width=\textwidth]{Fig_GRIFF_qcells}
     \caption{
       \protect\label{fig:qcells}
+      \protect\label{fig:TRIADS_qcells}
       Triad notation for quarter cells. $T$-cells are inside boxes,
       while the  $i+\fractext{1}{2},k$ $u$-cell is shaded in green and
 …
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:qcells}) with the quarter cell that is
+Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:TRIADS_qcells}) with the quarter cell that is
 the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ $u$-cell and the $i,k+k_p$ $w$-cell.
 Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:i13} and \autoref{eq:i31} in this notation,
+Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:TRIADS_i13} and \autoref{eq:TRIADS_i31} in this notation,
 we have \eg\ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$.
 Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to
 …
 and we notate these areas, similarly to the triad slopes,
 as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$,
 where \eg\ in \autoref{eq:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,
 and in \autoref{eq:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.
+where \eg\ in \autoref{eq:TRIADS_i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,
+and in \autoref{eq:TRIADS_i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.
 \subsection{Full triad fluxes}
 …
 tracer cell $i,k$ to $i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the form
 \begin{equation}
   \label{eq:i11}
+  \label{eq:TRIADS_i11}
   \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} =
   - \left( {A}_i^{k+1} a_{1} + {A}_i^{k+1} a_{2} + {A}_i^k
 …
   \frac{\delta_{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}},
 \end{equation}
 where the areas $a_i$ are as in \autoref{eq:i13}.
 In this case, separating the total lateral flux, the sum of \autoref{eq:i13} and \autoref{eq:i11},
+where the areas $a_i$ are as in \autoref{eq:TRIADS_i13}.
+In this case, separating the total lateral flux, the sum of \autoref{eq:TRIADS_i13} and \autoref{eq:TRIADS_i11},
 into triad components, a lateral tracer flux
 \begin{equation}
   \label{eq:latflux-triad}
+  \label{eq:TRIADS_latflux-triad}
   _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - {A}_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p}
   \left(
 …
 the lateral density flux associated with each triad separately disappears.
 \begin{equation}
   \label{eq:latflux-rho}
+  \label{eq:TRIADS_latflux-rho}
   {\mathbb{F}_u}_{i_p}^{k_p} (\rho)=-\alpha _i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (S)=0
 \end{equation}
 …
 tracer cell $i,k$ to $i+1,k$ must also vanish since it is a sum of four such triad fluxes.
 The squared slope $r_1^2$ in the expression \autoref{eq:i33c} for the $_{33}$ component is also expressed in
+The squared slope $r_1^2$ in the expression \autoref{eq:TRIADS_i33c} for the $_{33}$ component is also expressed in
 terms of area-weighted squared triad slopes,
 so the area-integrated vertical flux from tracer cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is
 \begin{equation}
   \label{eq:i33}
+  \label{eq:TRIADS_i33}
   \left( F_w^{33} \right) _i^{k+\frac{1}{2}} =
   - \left( {A}_i^{k+1} a_{1}' s_{1}'^2
 …
     + {A}_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right],
 \end{equation}
 where the areas $a'$ and slopes $s'$ are the same as in \autoref{eq:i31}.
 Then, separating the total vertical flux, the sum of \autoref{eq:i31} and \autoref{eq:i33},
+where the areas $a'$ and slopes $s'$ are the same as in \autoref{eq:TRIADS_i31}.
+Then, separating the total vertical flux, the sum of \autoref{eq:TRIADS_i31} and \autoref{eq:TRIADS_i33},
 into triad components, a vertical flux
 \begin{align}
   \label{eq:vertflux-triad}
+  \label{eq:TRIADS_vertflux-triad}
   _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T)
   &= {A}_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p}
 …
     \right) \\
   &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right)
     {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2}
+    {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:TRIADS_vertflux-triad2}
 \end{align}
 may be associated with each triad.
 …
 tracer cell $i,k$ to $i,k+1$ must also vanish since it is a sum of four such triad fluxes.
 We can explicitly identify (\autoref{fig:qcells}) the triads associated with the $s_i$, $a_i$,
 and $s'_i$, $a'_i$ used in the definition of the $u$-fluxes and $w$-fluxes in \autoref{eq:i31},
 \autoref{eq:i13}, \autoref{eq:i11} \autoref{eq:i33} and \autoref{fig:ISO_triad} to write out
+We can explicitly identify (\autoref{fig:TRIADS_qcells}) the triads associated with the $s_i$, $a_i$,
+and $s'_i$, $a'_i$ used in the definition of the $u$-fluxes and $w$-fluxes in \autoref{eq:TRIADS_i31},
+\autoref{eq:TRIADS_i13}, \autoref{eq:TRIADS_i11} \autoref{eq:TRIADS_i33} and \autoref{fig:TRIADS_ISO_triad} to write out
 the iso-neutral fluxes at $u$- and $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces:
 %(\autoref{fig:ISO_triad}):
+%(\autoref{fig:TRIADS_ISO_triad}):
 \begin{flalign}
   \label{eq:iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv
+  \label{eq:TRIADS_iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv
   \sum_{\substack{i_p,\,k_p}}
   \begin{pmatrix}
 …
 \subsection{Ensuring the scheme does not increase tracer variance}
 \label{subsec:variance}
+\label{subsec:TRIADS_variance}
 We now require that this operator should not increase the globally-integrated tracer variance.
 …
     &= -T_{i+i_p-1/2}^k{\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \quad + \quad  T_{i+i_p+1/2}^k
     {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\
     &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:dvar_iso_i}
+    &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:TRIADS_dvar_iso_i}
   \end{aligned}
 \end{multline}
 …
 the $T$-points $i,k+k_p-\fractext{1}{2}$ (above) and $i,k+k_p+\fractext{1}{2}$ (below) of
 \begin{equation}
   \label{eq:dvar_iso_k}
+  \label{eq:TRIADS_dvar_iso_k}
   _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
 \end{equation}
 The total variance tendency driven by the triad is the sum of these two.
 Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with
 \autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad}, it is
+\autoref{eq:TRIADS_latflux-triad} and \autoref{eq:TRIADS_vertflux-triad}, it is
 \begin{multline*}
   -{A}_i^k\left \{
 …
 be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by
 \begin{equation}
   \label{eq:V-A}
+  \label{eq:TRIADS_V-A}
   _i^k\mathbb{V}_{i_p}^{k_p}
   ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k}
 …
 the variance tendency reduces to the perfect square
 \begin{equation}
   \label{eq:perfect-square}
+  \label{eq:TRIADS_perfect-square}
   -{A}_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
   \left(
 …
   \right)^2\leq 0.
 \end{equation}
 Thus, the constraint \autoref{eq:V-A} ensures that the fluxes
 (\autoref{eq:latflux-triad}, \autoref{eq:vertflux-triad}) associated with
+Thus, the constraint \autoref{eq:TRIADS_V-A} ensures that the fluxes
+(\autoref{eq:TRIADS_latflux-triad}, \autoref{eq:TRIADS_vertflux-triad}) associated with
 a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase the net variance.
 Since the total fluxes are sums of such fluxes from the various triads, this constraint, applied to all triads,
 is sufficient to ensure that the globally integrated variance does not increase.
 The expression \autoref{eq:V-A} can be interpreted as a discretization of the global integral
 \begin{equation}
   \label{eq:cts-var}
+The expression \autoref{eq:TRIADS_V-A} can be interpreted as a discretization of the global integral
+\begin{equation}
+  \label{eq:TRIADS_cts-var}
   \frac{\partial}{\partial t}\int\!\fractext{1}{2} T^2\, dV =
   \int\!\mathbf{F}\cdot\nabla T\, dV,
 …
 \citet{griffies.gnanadesikan.ea_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells,
 defined in terms of the distances between $T$, $u$,$f$ and $w$-points.
 This is the natural discretization of \autoref{eq:cts-var}.
+This is the natural discretization of \autoref{eq:TRIADS_cts-var}.
 The \NEMO\ model, however, operates with scale factors instead of grid sizes,
 and scale factors for the quarter cells are not defined.
 Instead, therefore we simply choose
 \begin{equation}
   \label{eq:V-NEMO}
+  \label{eq:TRIADS_V-NEMO}
   _i^k\mathbb{V}_{i_p}^{k_p}=\fractext{1}{4} {b_u}_{i+i_p}^k,
 \end{equation}
 …
 the lateral flux from tracer cell $i,k$ to $i+1,k$ reduces to the classical form
 \begin{equation}
   \label{eq:lat-normal}
+  \label{eq:TRIADS_lat-normal}
   -\overline{A}_{\,i+1/2}^k\;
   \frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}}
 …
 In fact if the diffusive coefficient is defined at $u$-points,
 so that we employ ${A}_{i+i_p}^k$ instead of  ${A}_i^k$ in the definitions of the triad fluxes
 \autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad},
+\autoref{eq:TRIADS_latflux-triad} and \autoref{eq:TRIADS_vertflux-triad},
 we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above.
 …
 The iso-neutral fluxes at $u$- and $w$-points are the sums of the triad fluxes that
 cross the $u$- and $w$-faces \autoref{eq:iso_flux}:
+cross the $u$- and $w$-faces \autoref{eq:TRIADS_iso_flux}:
 \begin{subequations}
   % \label{eq:alltriadflux}
+  % \label{eq:TRIADS_alltriadflux}
   \begin{flalign*}
     % \label{eq:vect_isoflux}
+    % \label{eq:TRIADS_vect_isoflux}
     \vect{F}_{\mathrm{iso}}(T) &\equiv
     \sum_{\substack{i_p,\,k_p}}
 …
     \end{pmatrix},
   \end{flalign*}
   where \autoref{eq:latflux-triad}:
+  where \autoref{eq:TRIADS_latflux-triad}:
   \begin{align}
     \label{eq:triadfluxu}
+    \label{eq:TRIADS_triadfluxu}
     _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - {A}_i^k{
                                           \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}}
 …
                                           -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \
                                           \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
                                           \right),\label{eq:triadfluxw}
+                                          \right),\label{eq:TRIADS_triadfluxw}
   \end{align}
   with \autoref{eq:V-NEMO}
+  with \autoref{eq:TRIADS_V-NEMO}
   \[
     % \label{eq:V-NEMO2}
+    % \label{eq:TRIADS_V-NEMO2}
     _i^k{\mathbb{V}}_{i_p}^{k_p}=\fractext{1}{4} {b_u}_{i+i_p}^k.
   \]
 \end{subequations}
 The divergence of the expression \autoref{eq:iso_flux} for the fluxes gives the iso-neutral diffusion tendency at
+The divergence of the expression \autoref{eq:TRIADS_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at
 each tracer point:
 \[
   % \label{eq:iso_operator}
+  % \label{eq:TRIADS_iso_operator}
   D_l^T = \frac{1}{b_T}
   \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k
 …
 \item[$\bullet$ horizontal diffusion]
   The discretization of the diffusion operator recovers the traditional five-point Laplacian
   \autoref{eq:lat-normal} in the limit of flat iso-neutral direction:
+  \autoref{eq:TRIADS_lat-normal} in the limit of flat iso-neutral direction:
   \[
     % \label{eq:iso_property0}
+    % \label{eq:TRIADS_iso_property0}
     D_l^T = \frac{1}{b_T} \
     \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \;
 …
 \item[$\bullet$ implicit treatment in the vertical]
   Only tracer values associated with a single water column appear in the expression \autoref{eq:i33} for
+  Only tracer values associated with a single water column appear in the expression \autoref{eq:TRIADS_i33} for
   the $_{33}$ fluxes, vertical fluxes driven by vertical gradients.
   This is of paramount importance since it means that a time-implicit algorithm can be used to
 …
 \item[$\bullet$ pure iso-neutral operator]
   The iso-neutral flux of locally referenced potential density is zero.
   See \autoref{eq:latflux-rho} and \autoref{eq:vertflux-triad2}.
+  See \autoref{eq:TRIADS_latflux-rho} and \autoref{eq:TRIADS_vertflux-triad2}.
 \item[$\bullet$ conservation of tracer]
   The iso-neutral diffusion conserves tracer content, \ie
   \[
     % \label{eq:iso_property1}
+    % \label{eq:TRIADS_iso_property1}
     \sum_{i,j,k} \left\{ D_l^T \      b_T \right\} = 0
   \]
 …
   The iso-neutral diffusion does not increase the tracer variance, \ie
   \[
     % \label{eq:iso_property2}
+    % \label{eq:TRIADS_iso_property2}
     \sum_{i,j,k} \left\{ T \ D_l^T      \ b_T \right\} \leq 0
   \]
   The property is demonstrated in \autoref{subsec:variance} above.
+  The property is demonstrated in \autoref{subsec:TRIADS_variance} above.
   It is a key property for a diffusion term.
   It means that it is also a dissipation term,
 …
   The iso-neutral diffusion operator is self-adjoint, \ie
   \begin{equation}
     \label{eq:iso_property3}
+    \label{eq:TRIADS_iso_property3}
     \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}
   \end{equation}
 …
   We just have to apply the same routine.
   This property can be demonstrated similarly to the proof of the `no increase of tracer variance' property.
   The contribution by a single triad towards the left hand side of \autoref{eq:iso_property3},
   can be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:dvar_iso_i} and \autoref{eq:dvar_iso_k}.
   This results in a term similar to \autoref{eq:perfect-square},
+  The contribution by a single triad towards the left hand side of \autoref{eq:TRIADS_iso_property3},
+  can be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:TRIADS_dvar_iso_i} and \autoref{eq:TRIADS_dvar_iso_k}.
+  This results in a term similar to \autoref{eq:TRIADS_perfect-square},
   \[
     % \label{eq:TScovar}
+    % \label{eq:TRIADS_TScovar}
     - {A}_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
     \left(
 …
   \]
 This is symmetrical in $T $ and $S$, so exactly the same term arises from
 the discretization of this triad's contribution towards the RHS of \autoref{eq:iso_property3}.
+the discretization of this triad's contribution towards the RHS of \autoref{eq:TRIADS_iso_property3}.
 \end{description}
 \subsection{Treatment of the triads at the boundaries}
 \label{sec:iso_bdry}
+\label{sec:TRIADS_iso_bdry}
 The triad slope can only be defined where both the grid boxes centred at the end of the arms exist.
 Triads that would poke up through the upper ocean surface into the atmosphere,
 or down into the ocean floor, must be masked out.
 See \autoref{fig:bdry_triads}.
+See \autoref{fig:TRIADS_bdry_triads}.
 Surface layer triads \triad{i}{1}{R}{1/2}{-1/2} (magenta) and \triad{i+1}{1}{R}{-1/2}{-1/2} (blue) that
 require density to be specified above the ocean surface are masked (\autoref{fig:bdry_triads}a):
+require density to be specified above the ocean surface are masked (\autoref{fig:TRIADS_bdry_triads}a):
 this ensures that lateral tracer gradients produce no flux through the ocean surface.
 However, to prevent surface noise, it is customary to retain the $_{11}$ contributions towards
 the lateral triad fluxes \triad[u]{i}{1}{F}{1/2}{-1/2} and \triad[u]{i+1}{1}{F}{-1/2}{-1/2};
 this drives diapycnal tracer fluxes.
 Similar comments apply to triads that would intersect the ocean floor (\autoref{fig:bdry_triads}b).
+Similar comments apply to triads that would intersect the ocean floor (\autoref{fig:TRIADS_bdry_triads}b).
 Note that both near bottom triad slopes \triad{i}{k}{R}{1/2}{1/2} and \triad{i+1}{k}{R}{-1/2}{1/2} are masked when
 either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, \ie\ the $i,k+1$ $u$-point is masked.
 …
     \includegraphics[width=\textwidth]{Fig_GRIFF_bdry_triads}
     \caption{
       \protect\label{fig:bdry_triads}
+      \protect\label{fig:TRIADS_bdry_triads}
       (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer points (black dots),
       and $i+1/2,1$ $u$-point (blue square).
 …
 \subsection{ Limiting of the slopes within the interior}
 \label{sec:limit}
+\label{sec:TRIADS_limit}
 As discussed in \autoref{subsec:LDF_slp_iso},
 …
 It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to geopotentials
 (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to geopotentials)
 \autoref{eq:PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate surfaces, so we require
+\autoref{eq:MB_slopes_eiv} rather than the slope $r_i$ relative to coordinate surfaces, so we require
 \[
   |\tilde{r}_i|\leq \tilde{r}_\mathrm{max}=0.01.
 …
 Each individual triad slope
 \begin{equation}
   \label{eq:Rtilde}
+  \label{eq:TRIADS_Rtilde}
   _i^k\tilde{\mathbb{R}}_{i_p}^{k_p} = {}_i^k\mathbb{R}_{i_p}^{k_p}  + \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}
 \end{equation}
 …
 \subsection{Tapering within the surface mixed layer}
 \label{sec:taper}
+\label{sec:TRIADS_taper}
 Additional tapering of the iso-neutral fluxes is necessary within the surface mixed layer.
 …
 \subsubsection{Linear slope tapering within the surface mixed layer}
 \label{sec:lintaper}
+\label{sec:TRIADS_lintaper}
 This is the option activated by the default choice \np{ln\_triad\_iso}\forcode{ = .false.}.
 Slopes $\tilde{r}_i$ relative to geopotentials are tapered linearly from their value immediately below
 the mixed layer to zero at the surface, as described in option (c) of \autoref{fig:eiv_slp}, to values
 \begin{equation}
   \label{eq:rmtilde}
+the mixed layer to zero at the surface, as described in option (c) of \autoref{fig:LDF_eiv_slp}, to values
+\begin{equation}
+  \label{eq:TRIADS_rmtilde}
   \rMLt = -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h,
 \end{equation}
 and then the $r_i$ relative to vertical coordinate surfaces are appropriately adjusted to
 \[
   % \label{eq:rm}
+  % \label{eq:TRIADS_rm}
   \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h.
 \]
 Thus the diffusion operator within the mixed layer is given by:
 \[
   % \label{eq:iso_tensor_ML}
+  % \label{eq:TRIADS_iso_tensor_ML}
   D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
   \mbox{with}\quad \;\;\Re =\left( {{
 …
 in isopycnal layers immediately below, in the thermocline.
 It is consistent with the way the $\tilde{r}_i$ are tapered within the mixed layer
 (see \autoref{sec:taperskew} below) so as to ensure a uniform GM eddy-induced velocity throughout the mixed layer.
+(see \autoref{sec:TRIADS_taperskew} below) so as to ensure a uniform GM eddy-induced velocity throughout the mixed layer.
 However, it gives a downwards density flux and so acts so as to reduce potential energy in the same way as
 does the slope limiting discussed above in \autoref{sec:limit}.
 As in \autoref{sec:limit} above, the tapering \autoref{eq:rmtilde} is applied separately to
+does the slope limiting discussed above in \autoref{sec:TRIADS_limit}.
+As in \autoref{sec:TRIADS_limit} above, the tapering \autoref{eq:TRIADS_rmtilde} is applied separately to
 each triad $_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the $_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted.
 For clarity, we assume $z$-coordinates in the following;
 the conversion from $\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as
 described above by \autoref{eq:Rtilde}.
+described above by \autoref{eq:TRIADS_Rtilde}.
 \begin{enumerate}
 \item
   Mixed-layer depth is defined so as to avoid including regions of weak vertical stratification in
   the slope definition.
   At each $i,j$ (simplified to $i$ in \autoref{fig:MLB_triad}),
+  At each $i,j$ (simplified to $i$ in \autoref{fig:TRIADS_MLB_triad}),
   we define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer,
   $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) such that
   the potential density ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$,
   where $i,k_{10}$ is the tracer gridbox within which the depth reaches 10~m.
   See the left side of \autoref{fig:MLB_triad}.
+  See the left side of \autoref{fig:TRIADS_MLB_triad}.
   We use the $k_{10}$-gridbox instead of the surface gridbox to avoid problems \eg\ with thin daytime mixed-layers.
   Currently we use the same $\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is used to
 …
   This is to ensure that the vertical density gradients associated with
   these basal triad slopes ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are representative of the thermocline.
   The four basal triads defined in the bottom part of \autoref{fig:MLB_triad} are then
+  The four basal triads defined in the bottom part of \autoref{fig:TRIADS_MLB_triad} are then
   \begin{align*}
     {\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p} &=
                                                        {\:}^{k_{\mathrm{ML}}-k_p-1/2}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p},
                                                        % \label{eq:Rbase}
+                                                       % \label{eq:TRIADS_Rbase}
     \\
     \intertext{with \eg\ the green triad}
 …
 the $w$-point $i,k_{\mathrm{ML}}-1/2$ lying \emph{below} the $i,k_{\mathrm{ML}}$ tracer point, so it is this depth
 \[
   % \label{eq:zbase}
+  % \label{eq:TRIADS_zbase}
   {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2}
 \]
 one gridbox deeper than the diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper the slopes in
 \autoref{eq:rmtilde}.
+\autoref{eq:TRIADS_rmtilde}.
 \item
   Finally, we calculate the adjusted triads ${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within
 …
     {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p} &=
                                                        \frac{{z_w}_{k+k_p}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}.
                                                        % \label{eq:RML}
+                                                       % \label{eq:TRIADS_RML}
   \end{align*}
 \end{enumerate}
 …
 %  \fcapside {
   \caption{
     \protect\label{fig:MLB_triad}
+    \protect\label{fig:TRIADS_MLB_triad}
     Definition of mixed-layer depth and calculation of linearly tapered triads.
     The figure shows a water column at a given $i,j$ (simplified to $i$), with the ocean surface at the top.
 …
 \subsubsection{Additional truncation of skew iso-neutral flux components}
 \label{subsec:Gerdes-taper}
+\label{subsec:TRIADS_Gerdes-taper}
 The alternative option is activated by setting \np{ln\_triad\_iso} = true.
 …
 but replaces the $\rML$ in the skew term by
 \begin{equation}
   \label{eq:rm*}
+  \label{eq:TRIADS_rm*}
   \rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i,
 \end{equation}
 giving a ML diffusive operator
 \[
   % \label{eq:iso_tensor_ML2}
+  % \label{eq:TRIADS_iso_tensor_ML2}
   D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
   \mbox{with}\quad \;\;\Re =\left( {{
 …
 % ================================================================
 \section{Eddy induced advection formulated as a skew flux}
 \label{sec:skew-flux}
+\label{sec:TRIADS_skew-flux}
 \subsection{Continuous skew flux formulation}
 \label{sec:continuous-skew-flux}
+\label{sec:TRIADS_continuous-skew-flux}
 When Gent and McWilliams's [1990] diffusion is used, an additional advection term is added.
 …
 the formulation of which depends on the slopes of iso-neutral surfaces.
 Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces,
 \ie\ \autoref{eq:ldfslp_geo} is used in $z$-coordinate,
 and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates.
+\ie\ \autoref{eq:LDF_slp_geo} is used in $z$-coordinate,
+and the sum \autoref{eq:LDF_slp_geo} + \autoref{eq:LDF_slp_iso} in $z^*$ or $s$-coordinates.
 The eddy induced velocity is given by:
 \begin{subequations}
   % \label{eq:eiv}
+  % \label{eq:TRIADS_eiv}
   \begin{equation}
     \label{eq:eiv_v}
+    \label{eq:TRIADS_eiv_v}
     \begin{split}
       u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\
 …
   where the streamfunctions $\psi_i$ are given by
   \begin{equation}
     \label{eq:eiv_psi}
+    \label{eq:TRIADS_eiv_psi}
     \begin{split}
       \psi_1 & = A_{e} \; \tilde{r}_1,   \\
 …
 we end up with the skew form of the eddy induced advective fluxes per unit area in $ijk$ space:
 \begin{equation}
   \label{eq:eiv_skew_ijk}
+  \label{eq:TRIADS_eiv_skew_ijk}
   \textbf{F}_\mathrm{eiv}^T =
   \begin{pmatrix}
 …
 The total fluxes per unit physical area are then
 \begin{equation}
   \label{eq:eiv_skew_physical}
+  \label{eq:TRIADS_eiv_skew_physical}
   \begin{split}
     f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T   \\
 …
 \end{split}
 \end{equation}
 Note that \autoref{eq:eiv_skew_physical} takes the same form whatever the vertical coordinate,
 though of course the slopes $\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:eiv_psi} are relative to
+Note that \autoref{eq:TRIADS_eiv_skew_physical} takes the same form whatever the vertical coordinate,
+though of course the slopes $\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:TRIADS_eiv_psi} are relative to
 geopotentials.
 The tendency associated with eddy induced velocity is then simply the convergence of the fluxes
 (\autoref{eq:eiv_skew_ijk}, \autoref{eq:eiv_skew_physical}), so
 \[
   % \label{eq:skew_eiv_conv}
+(\autoref{eq:TRIADS_eiv_skew_ijk}, \autoref{eq:TRIADS_eiv_skew_physical}), so
+\[
+  % \label{eq:TRIADS_skew_eiv_conv}
   \frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[
     \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right)
 …
 \subsection{Discrete skew flux formulation}
 The skew fluxes in (\autoref{eq:eiv_skew_physical}, \autoref{eq:eiv_skew_ijk}),
 like the off-diagonal terms (\autoref{eq:i13c}, \autoref{eq:i31c}) of the small angle diffusion tensor,
 are best expressed in terms of the triad slopes, as in \autoref{fig:ISO_triad} and
 (\autoref{eq:i13}, \autoref{eq:i31});
+The skew fluxes in (\autoref{eq:TRIADS_eiv_skew_physical}, \autoref{eq:TRIADS_eiv_skew_ijk}),
+like the off-diagonal terms (\autoref{eq:TRIADS_i13c}, \autoref{eq:TRIADS_i31c}) of the small angle diffusion tensor,
+are best expressed in terms of the triad slopes, as in \autoref{fig:TRIADS_ISO_triad} and
+(\autoref{eq:TRIADS_i13}, \autoref{eq:TRIADS_i31});
 but now in terms of the triad slopes $\tilde{\mathbb{R}}$ relative to geopotentials instead of
 the $\mathbb{R}$ relative to coordinate surfaces.
 The discrete form of \autoref{eq:eiv_skew_ijk} using the slopes \autoref{eq:R} and
+The discrete form of \autoref{eq:TRIADS_eiv_skew_ijk} using the slopes \autoref{eq:TRIADS_R} and
 defining $A_e$ at $T$-points is then given by:
 \begin{subequations}
   % \label{eq:allskewflux}
+  % \label{eq:TRIADS_allskewflux}
   \begin{flalign*}
     % \label{eq:vect_skew_flux}
+    % \label{eq:TRIADS_vect_skew_flux}
     \vect{F}_{\mathrm{eiv}}(T) &\equiv    \sum_{\substack{i_p,\,k_p}}
     \begin{pmatrix}
 …
     \end{pmatrix},
   \end{flalign*}
   where the skew flux in the $i$-direction associated with a given triad is (\autoref{eq:latflux-triad},
   \autoref{eq:triadfluxu}):
+  where the skew flux in the $i$-direction associated with a given triad is (\autoref{eq:TRIADS_latflux-triad},
+  \autoref{eq:TRIADS_triadfluxu}):
   \begin{align}
     \label{eq:skewfluxu}
+    \label{eq:TRIADS_skewfluxu}
     _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \fractext{1}{4} {A_e}_i^k{
                                           \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}}
 …
                                           \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }, \\
     \intertext{
     and \autoref{eq:triadfluxw} in the $k$-direction, changing the sign
     to be consistent with \autoref{eq:eiv_skew_ijk}:
+    and \autoref{eq:TRIADS_triadfluxw} in the $k$-direction, changing the sign
+    to be consistent with \autoref{eq:TRIADS_eiv_skew_ijk}:
+    }
     _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T)
                                         &= -\fractext{1}{4} {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}}
                                           {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:skewfluxw}
+                                          {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:TRIADS_skewfluxw}
   \end{align}
 \end{subequations}
 …
 This can be seen %either from Appendix \autoref{apdx:eiv_skew} or
 by considering the fluxes associated with a given triad slope $_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$.
 For, following \autoref{subsec:variance} and \autoref{eq:dvar_iso_i},
+For, following \autoref{subsec:TRIADS_variance} and \autoref{eq:TRIADS_dvar_iso_i},
 the associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ drives a net rate of change of variance,
 summed over the two $T$-points $i+i_p-\fractext{1}{2},k$ and $i+i_p+\fractext{1}{2},k$, of
 \begin{equation}
   \label{eq:dvar_eiv_i}
+  \label{eq:TRIADS_dvar_eiv_i}
   _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k],
 \end{equation}
 …
 the $T$-points $i,k+k_p-\fractext{1}{2}$ (above) and $i,k+k_p+\fractext{1}{2}$ (below) of
 \begin{equation}
   \label{eq:dvar_eiv_k}
+  \label{eq:TRIADS_dvar_eiv_k}
   _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
 \end{equation}
 Inspection of the definitions (\autoref{eq:skewfluxu}, \autoref{eq:skewfluxw}) shows that
 these two variance changes (\autoref{eq:dvar_eiv_i}, \autoref{eq:dvar_eiv_k}) sum to zero.
+Inspection of the definitions (\autoref{eq:TRIADS_skewfluxu}, \autoref{eq:TRIADS_skewfluxw}) shows that
+these two variance changes (\autoref{eq:TRIADS_dvar_eiv_i}, \autoref{eq:TRIADS_dvar_eiv_k}) sum to zero.
 Hence the two fluxes associated with each triad make no net contribution to the variance budget.
 …
 For the change in gravitational PE driven by the $k$-flux is
 \begin{align}
   \label{eq:vert_densityPE}
+  \label{eq:TRIADS_vert_densityPE}
   g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho)
   &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k
     {\mathbb{S}_w}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k
     {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\
   \intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from \autoref{eq:skewfluxw}, gives}
+  \intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from \autoref{eq:TRIADS_skewfluxw}, gives}
   % and separating out
   % $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$,
 …
     \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}},
 \end{align}
 using the definition of the triad slope $\rtriad{R}$, \autoref{eq:R} to
+using the definition of the triad slope $\rtriad{R}$, \autoref{eq:TRIADS_R} to
 express $-\alpha _i^k\delta_{i+ i_p}[T^k]+\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of
 $-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$.
 …
 Where the coordinates slope, the $i$-flux gives a PE change
 \begin{multline}
   \label{eq:lat_densityPE}
+  \label{eq:TRIADS_lat_densityPE}
   g \delta_{i+i_p}[z_T^k]
   \left[
 …
   \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}},
 \end{multline}
 (using \autoref{eq:skewfluxu}) and so the total PE change \autoref{eq:vert_densityPE} +
 \autoref{eq:lat_densityPE} associated with the triad fluxes is
+(using \autoref{eq:TRIADS_skewfluxu}) and so the total PE change \autoref{eq:TRIADS_vert_densityPE} +
+\autoref{eq:TRIADS_lat_densityPE} associated with the triad fluxes is
 \begin{multline*}
   % \label{eq:tot_densityPE}
+  % \label{eq:TRIADS_tot_densityPE}
   g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) +
   g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\
 …
 \subsection{Treatment of the triads at the boundaries}
 \label{sec:skew_bdry}
 Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes are masked at the boundaries
 in exactly the same way as are the triad slopes \rtriad{R} used for the iso-neutral diffusive fluxes,
 as described in \autoref{sec:iso_bdry} and \autoref{fig:bdry_triads}.
 Thus surface layer triads $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are masked,
 and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when
 either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, \ie\ the $i,k+1$ $u$-point is masked.
+\label{sec:TRIADS_skew_bdry}
+Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes are masked at the boundaries
+in exactly the same way as are the triad slopes \rtriad{R} used for the iso-neutral diffusive fluxes,
+as described in \autoref{sec:TRIADS_iso_bdry} and \autoref{fig:TRIADS_bdry_triads}.
+Thus surface layer triads $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are masked,
+and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when
+either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, \ie\ the $i,k+1$ $u$-point is masked.
 The namelist parameter \np{ln\_botmix\_triad} has no effect on the eddy-induced skew-fluxes.
 \subsection{Limiting of the slopes within the interior}
 \label{sec:limitskew}
 Presently, the iso-neutral slopes $\tilde{r}_i$ relative to geopotentials are limited to be less than $1/100$,
 exactly as in calculating the iso-neutral diffusion, \S \autoref{sec:limit}.
+\label{sec:TRIADS_limitskew}
+Presently, the iso-neutral slopes $\tilde{r}_i$ relative to geopotentials are limited to be less than $1/100$,
+exactly as in calculating the iso-neutral diffusion, \S \autoref{sec:TRIADS_limit}.
 Each individual triad \rtriadt{R} is so limited.
 \subsection{Tapering within the surface mixed layer}
 \label{sec:taperskew}
 The slopes $\tilde{r}_i$ relative to geopotentials (and thus the individual triads \rtriadt{R})
 are always tapered linearly from their value immediately below the mixed layer to zero at the surface
 \autoref{eq:rmtilde}, as described in \autoref{sec:lintaper}.
 This is option (c) of \autoref{fig:eiv_slp}.
 This linear tapering for the slopes used to calculate the eddy-induced fluxes is unaffected by
+\label{sec:TRIADS_taperskew}
+The slopes $\tilde{r}_i$ relative to geopotentials (and thus the individual triads \rtriadt{R})
+are always tapered linearly from their value immediately below the mixed layer to zero at the surface
+\autoref{eq:TRIADS_rmtilde}, as described in \autoref{sec:TRIADS_lintaper}.
+This is option (c) of \autoref{fig:LDF_eiv_slp}.
+This linear tapering for the slopes used to calculate the eddy-induced fluxes is unaffected by
 the value of \np{ln\_triad\_iso}.
 …
 the horizontal (the most commonly used options in \NEMO: see \autoref{sec:LDF_coef}),
 it is equivalent to a horizontal eiv (eddy-induced velocity) that is uniform within the mixed layer
 \autoref{eq:eiv_v}.
+\autoref{eq:TRIADS_eiv_v}.
 This ensures that the eiv velocities do not restratify the mixed layer \citep{treguier.held.ea_JPO97,danabasoglu.ferrari.ea_JC08}.
 Equivantly, in terms of the skew-flux formulation we use here,
 …
 \subsection{Streamfunction diagnostics}
 \label{sec:sfdiag}
+\label{sec:TRIADS_sfdiag}
 Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.},
 …
 Each time step, streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
 $uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ (integer $i$, integer +1/2 $j$, integer +1/2 $k$)
 points (see Table \autoref{tab:cell}) respectively.
+points (see Table \autoref{tab:DOM_cell}) respectively.
 We follow \citep{griffies_bk04} and calculate the streamfunction at a given $uw$-point from
 the surrounding four triads according to:
 \[
   % \label{eq:sfdiagi}
+  % \label{eq:TRIADS_sfdiagi}
   {\psi_1}_{i+1/2}^{k+1/2}={\fractext{1}{4}}\sum_{\substack{i_p,\,k_p}}
   {A_e}_{i+1/2-i_p}^{k+1/2-k_p}\:\triadd{i+1/2-i_p}{k+1/2-k_p}{R}{i_p}{k_p}.
 \]
 The streamfunction $\psi_1$ is calculated similarly at $vw$ points.
 The eddy-induced velocities are then calculated from the straightforward discretisation of \autoref{eq:eiv_v}:
 \[
   % \label{eq:eiv_v_discrete}
+The eddy-induced velocities are then calculated from the straightforward discretisation of \autoref{eq:TRIADS_eiv_v}:
+\[
+  % \label{eq:TRIADS_eiv_v_discrete}
   \begin{split}
     {u^*}_{i+1/2}^{k} & = - \frac{1}{{e_{3u}}_{i}^{k}}\left({\psi_1}_{i+1/2}^{k+1/2}-{\psi_1}_{i+1/2}^{k+1/2}\right),   \\

NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex

-                      r11537
+                      r11543
 additional tendency terms to the prognostic equations:
 \begin{align*}
   % \label{eq:wa_traj_iau}
+  % \label{eq:ASM_wa_traj_iau}
   {\mathbf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\mathbf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\mathbf x}^{a}
 \end{align*}
 …
 The first function (namelist option \np{niaufn}=0) employs constant weights,
 \begin{align}
   \label{eq:F1_i}
+  \label{eq:ASM_F1_i}
   F^{(1)}_{i}
   =\left\{
 …
 with the weighting reduced linearly to a small value at the window end-points:
 \begin{align}
   \label{eq:F2_i}
+  \label{eq:ASM_F2_i}
   F^{(2)}_{i}
   =\left\{
 …
 \end{align}
 where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even.
 The weights described by \autoref{eq:F2_i} provide a smoother transition of the analysis trajectory from
 one assimilation cycle to the next than that described by \autoref{eq:F1_i}.
+The weights described by \autoref{eq:ASM_F2_i} provide a smoother transition of the analysis trajectory from
+one assimilation cycle to the next than that described by \autoref{eq:ASM_F1_i}.
 %==========================================================================
 …
 \begin{equation}
   \label{eq:asm_dmp}
+  \label{eq:ASM_dmp}
   \left\{
     \begin{aligned}
 …
 \[
   % \label{eq:asm_div}
+  % \label{eq:ASM_div}
   \chi^{n-1}_I = \frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
   \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u^{n-1}_I} \right]
 …
 \]
 By the application of \autoref{eq:asm_dmp} the divergence is filtered in each iteration,
+By the application of \autoref{eq:ASM_dmp} the divergence is filtered in each iteration,
 and the vorticity is left unchanged.
 In the presence of coastal boundaries with zero velocity increments perpendicular to the coast

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex

-                      r11537
+                      r11543
 \begin{table}
   \scriptsize
   \begin{tabularx}{\textwidth}{|X|c|c|c|}
+  \begin{tabular}{|l|c|c|}
     \hline
     tag ids affected by automatic definition of some of their attributes &
     name attribute                                                       &
     attribute value                      \\
+    attribute value                                                      \\
     \hline
     \hline
     field\_definition                                                    &
     freq\_op                                                             &
     \np{rn\_rdt}                         \\
+    \np{rn\_rdt}                                                         \\
     \hline
     SBC                                                                  &
     freq\_op                                                             &
     \np{rn\_rdt} $\times$ \np{nn\_fsbc}  \\
+    \np{rn\_rdt} $\times$ \np{nn\_fsbc}                                  \\
     \hline
     ptrc\_T                                                              &
     freq\_op                                                             &
     \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\
+    \np{rn\_rdt} $\times$ \np{nn\_dttrc}                                 \\
     \hline
     diad\_T                                                              &
     freq\_op                                                             &
     \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\
+    \np{rn\_rdt} $\times$ \np{nn\_dttrc}                                 \\
     \hline
     EqT, EqU, EqW                                                        &
     jbegin, ni,                                                          &
     according to the grid                \\
+    &
+    according to the grid                                                \\
+                                                                         &
     name\_suffix                                                         &
     \\
+                                                                         \\
     \hline
     TAO, RAMA and PIRATA moorings                                        &
     zoom\_ibegin, zoom\_jbegin,                                          &
     according to the grid                \\
+    &
+    according to the grid                                                \\
+                                                                         &
     name\_suffix                                                         &
     \\
     \hline
   \end{tabularx}
+                                                                         \\
+    \hline
+  \end{tabular}
 \end{table}
 …
 \subsection{XML reference tables}
 \label{subsec:IOM_xmlref}
+\label{subsec:DIA_IOM_xmlref}
 \begin{enumerate}
 …
 the CF metadata standard.
 Therefore while a user may wish to add their own metadata to the output files (as demonstrated in example 4 of
 section \autoref{subsec:IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard.
+section \autoref{subsec:DIA_IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard.
 Some metadata that may significantly increase the file size (horizontal cell areas and vertices) are controlled by
 …
 the mono-processor case (\ie\ global domain of {\small\ttfamily 182x149x31}).
 An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in
 table \autoref{tab:NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with
+table \autoref{tab:DIA_NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with
 a 4x2 mpi partitioning.
 Note the variation in the compression ratio achieved which reflects chiefly the dry to wet volume ratio of
 …
   \end{tabular}
   \caption{
     \protect\label{tab:NC4}
+    \protect\label{tab:DIA_NC4}
     Filesize comparison between NetCDF3 and NetCDF4 with chunking and compression
+  }
 …
 \section[FLO: On-Line Floats trajectories (\texttt{\textbf{key\_floats}})]
 {FLO: On-Line Floats trajectories (\protect\key{floats})}
 \label{sec:FLO}
+\label{sec:DIA_FLO}
 %--------------------------------------------namflo-------------------------------------------------------
 …
     \mathcal{V} &=  \mathcal{A}  \;\bar{\eta}
   \end{split}
   \label{eq:MV_nBq}
+  \label{eq:DIA_MV_nBq}
 \end{equation}
 …
   \frac{1}{e_3} \partial_t ( e_3\,\rho) + \nabla( \rho \, \textbf{U} )
   = \left. \frac{\textit{emp}}{e_3}\right|_\textit{surface}
   \label{eq:Co_nBq}
+  \label{eq:DIA_Co_nBq}
 \end{equation}
 …
 \begin{equation}
   \partial_t \mathcal{M} = \mathcal{A} \;\overline{\textit{emp}}
   \label{eq:Mass_nBq}
+  \label{eq:DIA_Mass_nBq}
 \end{equation}
 where $\overline{\textit{emp}} = \int_S \textit{emp}\,ds$ is the net mass flux through the ocean surface.
 Bringing \autoref{eq:Mass_nBq} and the time derivative of \autoref{eq:MV_nBq} together leads to
+Bringing \autoref{eq:DIA_Mass_nBq} and the time derivative of \autoref{eq:DIA_MV_nBq} together leads to
 the evolution equation of the mean sea level
 …
   \partial_t \bar{\eta} = \frac{\overline{\textit{emp}}}{ \bar{\rho}}
   - \frac{\mathcal{V}}{\mathcal{A}} \;\frac{\partial_t \bar{\rho} }{\bar{\rho}}
   \label{eq:ssh_nBq}
+  \label{eq:DIA_ssh_nBq}
 \end{equation}
 The first term in equation \autoref{eq:ssh_nBq} alters sea level by adding or subtracting mass from the ocean.
+The first term in equation \autoref{eq:DIA_ssh_nBq} alters sea level by adding or subtracting mass from the ocean.
 The second term arises from temporal changes in the global mean density; \ie\ from steric effects.
 In a Boussinesq fluid, $\rho$ is replaced by $\rho_o$ in all the equation except when $\rho$ appears multiplied by
 the gravity (\ie\ in the hydrostatic balance of the primitive Equations).
 In particular, the mass conservation equation, \autoref{eq:Co_nBq}, degenerates into the incompressibility equation:
+In particular, the mass conservation equation, \autoref{eq:DIA_Co_nBq}, degenerates into the incompressibility equation:
 \[
   \frac{1}{e_3} \partial_t ( e_3 ) + \nabla( \textbf{U} ) = \left. \frac{\textit{emp}}{\rho_o \,e_3}\right|_ \textit{surface}
   % \label{eq:Co_Bq}
+  % \label{eq:DIA_Co_Bq}
 \]
 …
 \[
   \partial_t \mathcal{V} = \mathcal{A} \;\frac{\overline{\textit{emp}}}{\rho_o}
   % \label{eq:V_Bq}
+  % \label{eq:DIA_V_Bq}
 \]
 …
 \begin{equation}
   \mathcal{M}_o = \mathcal{M} + \rho_o \,\eta_s \,\mathcal{A}
   \label{eq:M_Bq}
+  \label{eq:DIA_M_Bq}
 \end{equation}
 …
 Introducing the total density anomaly, $\mathcal{D}= \int_D d_a \,dv$,
 where $d_a = (\rho -\rho_o ) / \rho_o$ is the density anomaly used in \NEMO\ (cf. \autoref{subsec:TRA_eos})
 in \autoref{eq:M_Bq} leads to a very simple form for the steric height:
+in \autoref{eq:DIA_M_Bq} leads to a very simple form for the steric height:
 \begin{equation}
   \eta_s = - \frac{1}{\mathcal{A}} \mathcal{D}
   \label{eq:steric_Bq}
+  \label{eq:DIA_steric_Bq}
 \end{equation}
 …
 (wetting and drying of grid point is not allowed).
 Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of free surface which is considered.
+Third, the discretisation of \autoref{eq:DIA_steric_Bq} depends on the type of free surface which is considered.
 In the non linear free surface case, \ie\ \np{ln\_linssh}\forcode{=.true.}, it is given by
 \[
   \eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t} e_{2t} e_{3t} }{ \sum_{i,\,j,\,k}       e_{1t} e_{2t} e_{3t} }
   % \label{eq:discrete_steric_Bq_nfs}
+  % \label{eq:DIA_discrete_steric_Bq_nfs}
 \]
 …
   \eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t}e_{2t}e_{3t} + \sum_{i,\,j} d_a\; e_{1t}e_{2t} \eta }
                   { \sum_{i,\,j,\,k}       e_{1t}e_{2t}e_{3t} + \sum_{i,\,j}       e_{1t}e_{2t} \eta }
   % \label{eq:discrete_steric_Bq_fs}
+  % \label{eq:DIA_discrete_steric_Bq_fs}
 \]
 …
 \[
   \eta_s = - \frac{1}{\mathcal{A}} \int_D d_a(T,S_o,p_o) \,dv
   % \label{eq:thermosteric_Bq}
+  % \label{eq:DIA_thermosteric_Bq}
 \]
 …
     \includegraphics[width=\textwidth]{Fig_mask_subasins}
     \caption{
       \protect\label{fig:mask_subasins}
+      \protect\label{fig:DIA_mask_subasins}
       Decomposition of the World Ocean (here ORCA2) into sub-basin used in to
       compute the heat and salt transports as well as the meridional stream-function:
 …
 Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean.
 The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays,
 the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}).
+the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:DIA_mask_subasins}).
 %------------------------------------------namptr-----------------------------------------
 …
 \[
   C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}}
   % \label{eq:CFL}
+  % \label{eq:DIA_CFL}
 \]

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex

-                      r11537
+                      r11543
     {\em
       Compatibility changes Major simplification has moved many of the options to external domain configuration tools.
       (see \autoref{apdx:DOMAINcfg})
+      (see \autoref{apdx:DOMCFG})
     }                                                                                            \\
     {\em 3.x} & {\em Rachid Benshila, Gurvan Madec \& S\'{e}bastien Masson} &
 …
 \newpage
 Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretisation \autoref{chap:STP},
+Having defined the continuous equations in \autoref{chap:MB} and chosen a time discretisation \autoref{chap:TD},
 we need to choose a grid for spatial discretisation and related numerical algorithms.
 In the present chapter, we provide a general description of the staggered grid used in \NEMO,
 …
     \includegraphics[width=\textwidth]{Fig_cell}
     \caption{
       \protect\label{fig:cell}
+      \protect\label{fig:DOM_cell}
       Arrangement of variables.
       $t$ indicates scalar points where temperature, salinity, density, pressure and
 …
 The arrangement of variables is the same in all directions.
 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in
 the centre of each face of the cells (\autoref{fig:cell}).
+the centre of each face of the cells (\autoref{fig:DOM_cell}).
 This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification
 \citep{mesinger.arakawa_bk76}.
 …
 The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by the transformation that
 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}.
+The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:DOM_cell}.
 In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of
 the grid-point where the scale factors are defined.
 Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}.
+Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}.
 As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and
 $\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity.
 …
 centred finite difference approximation, not from their analytical expression.
 This preserves the symmetry of the discrete set of equations and therefore satisfies many of
 the continuous properties (see \autoref{apdx:C}).
+the continuous properties (see \autoref{apdx:INVARIANTS}).
 A similar, related remark can be made about the domain size:
 when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors
 …
     \end{tabular}
     \caption{
       \protect\label{tab:cell}
+      \protect\label{tab:DOM_cell}
       Location of grid-points as a function of integer or integer and a half value of the column, line or level.
       This indexing is only used for the writing of the semi -discrete equations.
 …
 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids
 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}.
 An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}.
+An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]
 …
     \includegraphics[width=\textwidth]{Fig_zgr_e3}
     \caption{
       \protect\label{fig:zgr_e3}
+      \protect\label{fig:DOM_zgr_e3}
       Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,
       and (b) analytically derived grid-point position and scale factors.
 …
 the midpoint between them are:
 \begin{alignat*}{2}
   % \label{eq:di_mi}
+  % \label{eq:DOM_di_mi}
   \delta_i [q]      &= &       &q (i + 1/2) - q (i - 1/2) \\
   \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2
 …
 Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$.
 Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at a $t$-point has
+Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, the gradient of a variable $q$ defined at a $t$-point has
 its three components defined at $u$-, $v$- and $w$-points while its Laplacian is defined at the $t$-point.
 These operators have the following discrete forms in the curvilinear $s$-coordinates system:
 …
 \end{multline*}
 Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at
+Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at
 vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and
 its divergence defined at $t$-points:
 …
 In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and
 $(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively.
 These two properties will be used extensively in the \autoref{apdx:C} to
+These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to
 demonstrate integral conservative properties of the discrete formulation chosen.
 …
     \includegraphics[width=\textwidth]{Fig_index_hor}
     \caption{
       \protect\label{fig:index_hor}
+      \protect\label{fig:DOM_index_hor}
       Horizontal integer indexing used in the \fortran code.
       The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices
 …
 \label{subsec:DOM_Num_Index_hor}
 The indexing in the horizontal plane has been chosen as shown in \autoref{fig:index_hor}.
+The indexing in the horizontal plane has been chosen as shown in \autoref{fig:DOM_index_hor}.
 For an increasing $i$ index ($j$ index),
 the $t$-point and the eastward $u$-point (northward $v$-point) have the same index
 (see the dashed area in \autoref{fig:index_hor}).
+(see the dashed area in \autoref{fig:DOM_index_hor}).
 A $t$-point and its nearest north-east $f$-point have the same $i$-and $j$-indices.
 …
 given in \autoref{subsec:DOM_cell}.
 The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below
 (\autoref{fig:index_vert}).
+(\autoref{fig:DOM_index_vert}).
 The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while
 the last $t$-level is always outside the ocean domain (\autoref{fig:index_vert}).
+the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}).
 Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index
 (\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction),
 in contrast to the indexing on the horizontal plane where the $t$-point has the same index as
 the nearest velocity points in the positive direction of the respective horizontal axis index
 (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}).
+(compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}).
 Since the scale factors are chosen to be strictly positive,
 a \textit{minus sign} is included in the \fortran implementations of
 …
     \includegraphics[width=\textwidth]{Fig_index_vert}
     \caption{
       \protect\label{fig:index_vert}
+      \protect\label{fig:DOM_index_vert}
       Vertical integer indexing used in the \fortran code.
       Note that the $k$-axis is oriented downward.
 …
 the model domain itself can be altered by runtime selections.
 The code previously used to perform vertical discretisation has been incorporated into an external tool
 (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMAINcfg}.
+(\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}.
 The next subsections summarise the parameter and fields related to the configuration of the whole model domain.
 …
 The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to
 the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$,
 evaluated at the values as specified in \autoref{tab:cell} for the respective grid-point position.
+evaluated at the values as specified in \autoref{tab:DOM_cell} for the respective grid-point position.
 The calculation of the values of the horizontal scale factor arrays in general additionally involves
 partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$,
 …
     \includegraphics[width=\textwidth]{Fig_z_zps_s_sps}
     \caption{
       \protect\label{fig:z_zps_s_sps}
+      \protect\label{fig:DOM_z_zps_s_sps}
       The ocean bottom as seen by the model:
       (a) $z$-coordinate with full step,
 …
 By default a non-linear free surface is used (\np{ln\_linssh} set to \forcode{=.false.} in \nam{dom}):
 the coordinate follow the time-variation of the free surface so that the transformation is time dependent:
 $z(i,j,k,t)$ (\eg\ \autoref{fig:z_zps_s_sps}f).
+$z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f).
 When a linear free surface is assumed (\np{ln\_linssh} set to \forcode{=.true.} in \nam{dom}),
 the vertical coordinates are fixed in time, but the seawater can move up and down across the $z_0$ surface
 …
 \medskip
 The decision on these choices must be made when the \np{cn\_domcfg} file is constructed.
 Three main choices are offered (\autoref{fig:z_zps_s_sps}a-c):
+Three main choices are offered (\autoref{fig:DOM_z_zps_s_sps}a-c):
 \begin{itemize}
 …
 Additionally, hybrid combinations of the three main coordinates are available:
 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps}d and \autoref{fig:z_zps_s_sps}e).
+$s-z$ or $s-zps$ coordinate (\autoref{fig:DOM_z_zps_s_sps}d and \autoref{fig:DOM_z_zps_s_sps}e).
 A further choice related to vertical coordinate concerns
 …
 \section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]
 {Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})}
 \label{sec:DTA_tsd}
+\label{sec:DOM_DTA_tsd}
 %-----------------------------------------namtsd-------------------------------------------
 \nlst{namtsd}
 …
   Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}.
   The default version sets horizontally uniform T and profiles as used in the GYRE configuration
   (see \autoref{sec:CFG_gyre}).
+  (see \autoref{sec:CFGS_gyre}).
 \end{description}

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex

-                      r11537
+                      r11543
 %           Horizontal divergence and relative vorticity
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]
+{Horizontal divergence and relative vorticity (\protect\mdl{divcur})}
+\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})}
 \label{subsec:DYN_divcur}
 The vorticity is defined at an $f$-point (\ie\ corner point) as follows:
 \begin{equation}
   \label{eq:divcur_cur}
+  \label{eq:DYN_divcur_cur}
   \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right]
       -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right)
 …
 It is given by:
 \[
   % \label{eq:divcur_div}
+  % \label{eq:DYN_divcur_div}
   \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
   \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right]
 …
 %           Sea Surface Height evolution
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]
+{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}
+\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}
 \label{subsec:DYN_sshwzv}
 The sea surface height is given by:
 \begin{equation}
   \label{eq:dynspg_ssh}
+  \label{eq:DYN_spg_ssh}
   \begin{aligned}
     \frac{\partial \eta }{\partial t}
 …
 \textit{emp} can be written as the evaporation minus precipitation, minus the river runoff.
 The sea-surface height is evaluated using exactly the same time stepping scheme as
 the tracer equation \autoref{eq:tra_nxt}:
+the tracer equation \autoref{eq:TRA_nxt}:
 a leapfrog scheme in combination with an Asselin time filter,
 \ie\ the velocity appearing in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity).
+\ie\ the velocity appearing in \autoref{eq:DYN_spg_ssh} is centred in time (\textit{now} velocity).
 This is of paramount importance.
 Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to
 …
 taking into account the change of the thickness of the levels:
 \begin{equation}
   \label{eq:wzv}
+  \label{eq:DYN_wzv}
   \left\{
     \begin{aligned}
 …
 re-orientated downward.
 \gmcomment{not sure of this...  to be modified with the change in emp setting}
 In the case of a linear free surface, the time derivative in \autoref{eq:wzv} disappears.
+In the case of a linear free surface, the time derivative in \autoref{eq:DYN_wzv} disappears.
 The upper boundary condition applies at a fixed level $z=0$.
 The top vertical velocity is thus equal to the divergence of the barotropic transport
 (\ie\ the first term in the right-hand-side of \autoref{eq:dynspg_ssh}).
+(\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}).
 Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates,
 …
 in the second case, $w$ is the velocity normal to the $s$-surfaces.
 Note also that the $k$-axis is re-orientated downwards in the \fortran code compared to
 the indexing used in the semi-discrete equations such as \autoref{eq:wzv}
+the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv}
 (see \autoref{subsec:DOM_Num_Index_vertical}).
 …
 %        Vorticity term
 % -------------------------------------------------------------------------------------------------------------
+\subsection[Vorticity term (\textit{dynvor.F90})]
+{Vorticity term (\protect\mdl{dynvor})}
+\subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})}
 \label{subsec:DYN_vor}
 %------------------------------------------nam_dynvor----------------------------------------------------
 …
 horizontal kinetic energy for the planetary vorticity term (MIX scheme);
 or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy
 (EEN scheme) (see \autoref{subsec:C_vorEEN}).
+(EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}).
 In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of
 vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{=.true.}).
 …
 %                 enstrophy conserving scheme
 %-------------------------------------------------------------
+\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens=.true.})]
+{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{=.true.})}
+\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens = .true.})]{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})}
 \label{subsec:DYN_vor_ens}
 …
 It is given by:
 \begin{equation}
   \label{eq:dynvor_ens}
+  \label{eq:DYN_vor_ens}
   \left\{
     \begin{aligned}
 …
 %                 energy conserving scheme
 %-------------------------------------------------------------
+\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene=.true.})]
+{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{=.true.})}
+\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene = .true.})]{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})}
 \label{subsec:DYN_vor_ene}
 …
 It is given by:
 \begin{equation}
   \label{eq:dynvor_ene}
+  \label{eq:DYN_vor_ene}
   \left\{
     \begin{aligned}
 …
 %                 mix energy/enstrophy conserving scheme
 %-------------------------------------------------------------
+\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix=.true.})]
+{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{=.true.})}
+\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix = .true.})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.})}
 \label{subsec:DYN_vor_mix}
 For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used.
 It consists of the ENS scheme (\autoref{eq:dynvor_ens}) for the relative vorticity term,
 and of the ENE scheme (\autoref{eq:dynvor_ene}) applied to the planetary vorticity term.
 \[
   % \label{eq:dynvor_mix}
+It consists of the ENS scheme (\autoref{eq:DYN_vor_ens}) for the relative vorticity term,
+and of the ENE scheme (\autoref{eq:DYN_vor_ene}) applied to the planetary vorticity term.
+\[
+  % \label{eq:DYN_vor_mix}
   \left\{ {
       \begin{aligned}
 …
 %                 energy and enstrophy conserving scheme
 %-------------------------------------------------------------
+\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een=.true.})]
+{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{=.true.})}
+\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een = .true.})]{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})}
 \label{subsec:DYN_vor_een}
 …
 The idea is to get rid of the double averaging by considering triad combinations of vorticity.
 It is noteworthy that this solution is conceptually quite similar to the one proposed by
 \citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}).
+\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:INVARIANTS}).
 The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified
 …
 First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:
 \[
   % \label{eq:pot_vor}
+  % \label{eq:DYN_pot_vor}
   q  = \frac{\zeta +f} {e_{3f} }
 \]
 where the relative vorticity is defined by (\autoref{eq:divcur_cur}),
+where the relative vorticity is defined by (\autoref{eq:DYN_divcur_cur}),
 the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is:
 \begin{equation}
   \label{eq:een_e3f}
+  \label{eq:DYN_een_e3f}
   e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2}
 \end{equation}
 …
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 A key point in \autoref{eq:een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.
+A key point in \autoref{eq:DYN_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.
 It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks
 (\np{nn\_een\_e3f}\forcode{=1}), or just by $4$ (\np{nn\_een\_e3f}\forcode{=.true.}).
 …
 (\autoref{fig:DYN_een_triad}):
 \begin{equation}
   \label{eq:Q_triads}
+  \label{eq:DYN_Q_triads}
   _i^j \mathbb{Q}^{i_p}_{j_p}
   = \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right)
 …
 Finally, the vorticity terms are represented as:
 \begin{equation}
   \label{eq:dynvor_een}
+  \label{eq:DYN_vor_een}
   \left\{ {
       \begin{aligned}
 …
 This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes.
 It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow
 (\ie\ $\chi$=$0$) (see \autoref{subsec:C_vorEEN}).
+(\ie\ $\chi$=$0$) (see \autoref{subsec:INVARIANTS_vorEEN}).
 Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of
 the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}.
 …
 %           Kinetic Energy Gradient term
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]
+{Kinetic energy gradient term (\protect\mdl{dynkeg})}
+\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})}
 \label{subsec:DYN_keg}
 As demonstrated in \autoref{apdx:C},
+As demonstrated in \autoref{apdx:INVARIANTS},
 there is a single discrete formulation of the kinetic energy gradient term that,
 together with the formulation chosen for the vertical advection (see below),
 conserves the total kinetic energy:
 \[
   % \label{eq:dynkeg}
+  % \label{eq:DYN_keg}
   \left\{
     \begin{aligned}
 …
 %           Vertical advection term
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Vertical advection term (\textit{dynzad.F90})]
+{Vertical advection term (\protect\mdl{dynzad})}
+\subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})}
 \label{subsec:DYN_zad}
 …
 conserves the total kinetic energy.
 Indeed, the change of KE due to the vertical advection is exactly balanced by
 the change of KE due to the gradient of KE (see \autoref{apdx:C}).
 \[
   % \label{eq:dynzad}
+the change of KE due to the gradient of KE (see \autoref{apdx:INVARIANTS}).
+\[
+  % \label{eq:DYN_zad}
   \left\{
     \begin{aligned}
 …
 %           Coriolis plus curvature metric terms
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]
+{Coriolis plus curvature metric terms (\protect\mdl{dynvor})}
+\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})}
 \label{subsec:DYN_cor_flux}
 …
 It is given by:
 \begin{multline*}
   % \label{eq:dyncor_metric}
+  % \label{eq:DYN_cor_metric}
   f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right)  \\
   \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right]
 …
 \end{multline*}
 Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) schemes can be used to
+Any of the (\autoref{eq:DYN_vor_ens}), (\autoref{eq:DYN_vor_ene}) and (\autoref{eq:DYN_vor_een}) schemes can be used to
 compute the product of the Coriolis parameter and the vorticity.
 However, the energy-conserving scheme (\autoref{eq:dynvor_een}) has exclusively been used to date.
+However, the energy-conserving scheme (\autoref{eq:DYN_vor_een}) has exclusively been used to date.
 This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity).
 …
 %           Flux form Advection term
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Flux form advection term (\textit{dynadv.F90})]
+{Flux form advection term (\protect\mdl{dynadv})}
+\subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})}
 \label{subsec:DYN_adv_flux}
 The discrete expression of the advection term is given by:
 \[
   % \label{eq:dynadv}
+  % \label{eq:DYN_adv}
   \left\{
     \begin{aligned}
 …
 %                 2nd order centred scheme
 %-------------------------------------------------------------
+\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2=.true.})]
+{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{=.true.})}
+\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2 = .true.})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})}
 \label{subsec:DYN_adv_cen2}
 In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points:
 \begin{equation}
   \label{eq:dynadv_cen2}
+  \label{eq:DYN_adv_cen2}
   \left\{
     \begin{aligned}
 …
 %                 UBS scheme
 %-------------------------------------------------------------
+\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs=.true.})]
+{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{=.true.})}
+\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs = .true.})]{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})}
 \label{subsec:DYN_adv_ubs}
 …
 For example, the evaluation of $u_T^{ubs} $ is done as follows:
 \begin{equation}
   \label{eq:dynadv_ubs}
+  \label{eq:DYN_adv_ubs}
   u_T^{ubs} =\overline u ^i-\;\frac{1}{6}
   \begin{cases}
 …
 The UBS scheme is not used in all directions.
 In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie\ $u_{uw}^{ubs}$ and
 $u_{vw}^{ubs}$ in \autoref{eq:dynadv_cen2} are used.
+$u_{vw}^{ubs}$ in \autoref{eq:DYN_adv_cen2} are used.
 UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm  pursue the
 sentence:Since vertical mixing of momentum is a source term of the TKE equation...  }
 For stability reasons, the first term in (\autoref{eq:dynadv_ubs}),
+For stability reasons, the first term in (\autoref{eq:DYN_adv_ubs}),
 which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time),
 while the second term, which is the diffusion part of the scheme,
 …
 Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by
 one coefficient.
 Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
+Replacing $1/6$ by $1/8$ in (\autoref{eq:DYN_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
 This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded.
 Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme.
 …
 %           Hydrostatic pressure gradient term
 % ================================================================
+\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]
+{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
+\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
 \label{sec:DYN_hpg}
 %------------------------------------------nam_dynhpg---------------------------------------------------
 …
 %           z-coordinate with full step
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco=.true.})]
+{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{=.true.})}
+\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco = .true.})]{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})}
 \label{subsec:DYN_hpg_zco}
 …
 for $k=km$ (surface layer, $jk=1$ in the code)
 \begin{equation}
   \label{eq:dynhpg_zco_surf}
+  \label{eq:DYN_hpg_zco_surf}
   \left\{
     \begin{aligned}
 …
 for $1<k<km$ (interior layer)
 \begin{equation}
   \label{eq:dynhpg_zco}
+  \label{eq:DYN_hpg_zco}
   \left\{
     \begin{aligned}
 …
 \end{equation}
 Note that the $1/2$ factor in (\autoref{eq:dynhpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as
+Note that the $1/2$ factor in (\autoref{eq:DYN_hpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as
 the vertical derivative of the scale factor at the surface level ($z=0$).
 Note also that in case of variable volume level (\texttt{vvl?} defined),
 the surface pressure gradient is included in \autoref{eq:dynhpg_zco_surf} and
 \autoref{eq:dynhpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$.
+the surface pressure gradient is included in \autoref{eq:DYN_hpg_zco_surf} and
+\autoref{eq:DYN_hpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$.
 %--------------------------------------------------------------------------------------------------------------
 %           z-coordinate with partial step
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps=.true.})]
+{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{=.true.})}
+\subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps = .true.})]{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})}
 \label{subsec:DYN_hpg_zps}
 …
 $\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{=.true.})
 \begin{equation}
   \label{eq:dynhpg_sco}
+  \label{eq:DYN_hpg_sco}
   \left\{
     \begin{aligned}
 …
 Where the first term is the pressure gradient along coordinates,
 computed as in \autoref{eq:dynhpg_zco_surf} - \autoref{eq:dynhpg_zco},
+computed as in \autoref{eq:DYN_hpg_zco_surf} - \autoref{eq:DYN_hpg_zco},
 and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point
 ($e_{3w}$).
 …
 (\np{ln\_dynhpg\_djc}\forcode{=.true.}) (currently disabled; under development)
 Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is activated
+Note that expression \autoref{eq:DYN_hpg_sco} is commonly used when the variable volume formulation is activated
 (\texttt{vvl?}) because in that case, even with a flat bottom,
 the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}.
 …
 \subsection{Ice shelf cavity}
 \label{subsec:DYN_hpg_isf}
 Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and
 the pressure gradient due to the ocean load (\np{ln\_dynhpg\_isf}\forcode{=.true.}).\\
 …
 A detailed description of this method is described in \citet{losch_JGR08}.\\
 The pressure gradient due to ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in
+The pressure gradient due to ocean load is computed using the expression \autoref{eq:DYN_hpg_sco} described in
 \autoref{subsec:DYN_hpg_sco}.
 …
 %           Time-scheme
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Time-scheme (\forcode{ln_dynhpg_imp={.true.,.false.}})]
+{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{=.true.,.false.})}
+\subsection[Time-scheme (\forcode{ln_dynhpg_imp = .{true,false}.})]{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .\{true,false\}}.)}
 \label{subsec:DYN_hpg_imp}
 …
 \begin{equation}
   \label{eq:dynhpg_lf}
+  \label{eq:DYN_hpg_lf}
   \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
   -\frac{1}{\rho_o \,e_{1u} }\delta_{i+1/2} \left[ {p_h^t } \right]
 …
 $\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}\forcode{=.true.}):
 \begin{equation}
   \label{eq:dynhpg_imp}
+  \label{eq:DYN_hpg_imp}
   \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
   -\frac{1}{4\,\rho_o \,e_{1u} } \delta_{i+1/2} \left[ p_h^{t+\rdt} +2\,p_h^t +p_h^{t-\rdt}  \right]
 \end{equation}
 The semi-implicit time scheme \autoref{eq:dynhpg_imp} is made possible without
+The semi-implicit time scheme \autoref{eq:DYN_hpg_imp} is made possible without
 significant additional computation since the density can be updated to time level $t+\rdt$ before
 computing the horizontal hydrostatic pressure gradient.
 It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using
 \autoref{eq:dynhpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:dynhpg_lf}.
 Note that \autoref{eq:dynhpg_imp} is equivalent to applying a time filter to the pressure gradient to
+\autoref{eq:DYN_hpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:DYN_hpg_lf}.
+Note that \autoref{eq:DYN_hpg_imp} is equivalent to applying a time filter to the pressure gradient to
 eliminate high frequency IGWs.
 Obviously, when using \autoref{eq:dynhpg_imp},
+Obviously, when using \autoref{eq:DYN_hpg_imp},
 the doubling of the time-step is achievable only if no other factors control the time-step,
 such as the stability limits associated with advection or diffusion.
 …
 The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows:
 \[
   % \label{eq:rho_flt}
+  % \label{eq:DYN_rho_flt}
   \rho^t = \rho( \widetilde{T},\widetilde {S},z_t)
   \quad    \text{with}  \quad
 …
 % Surface Pressure Gradient
 % ================================================================
+\section[Surface pressure gradient (\textit{dynspg.F90})]
+{Surface pressure gradient (\protect\mdl{dynspg})}
+\section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})}
 \label{sec:DYN_spg}
 %-----------------------------------------nam_dynspg----------------------------------------------------
 …
 Options are defined through the \nam{dyn\_spg} namelist variables.
 The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:PE_hor_pg}).
+The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:MB_hor_pg}).
 The main distinction is between the fixed volume case (linear free surface) and
 the variable volume case (nonlinear free surface, \texttt{vvl?} is defined).
 In the linear free surface case (\autoref{subsec:PE_free_surface})
+In the linear free surface case (\autoref{subsec:MB_free_surface})
 the vertical scale factors $e_{3}$ are fixed in time,
 while they are time-dependent in the nonlinear case (\autoref{subsec:PE_free_surface}).
+while they are time-dependent in the nonlinear case (\autoref{subsec:MB_free_surface}).
 With both linear and nonlinear free surface, external gravity waves are allowed in the equations,
 which imposes a very small time step when an explicit time stepping is used.
 Two methods are proposed to allow a longer time step for the three-dimensional equations:
 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt?}),
+the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:MB_flt?}),
 and the split-explicit free surface described below.
 The extra term introduced in the filtered method is calculated implicitly,
 …
 The form of the surface pressure gradient term depends on how the user wants to
 handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:PE_hor_pg}).
+handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:MB_hor_pg}).
 Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx):
 an explicit formulation which requires a small time step;
 …
 % Explicit free surface formulation
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Explicit free surface (\texttt{ln\_dynspg\_exp}\forcode{=.true.})]
+{Explicit free surface (\protect\np{ln\_dynspg\_exp}\forcode{=.true.})}
+\subsection[Explicit free surface (\texttt{ln\_dynspg\_exp}\forcode{ = .true.})]{Explicit free surface (\protect\np{ln\_dynspg\_exp}\forcode{ = .true.})}
 \label{subsec:DYN_spg_exp}
 …
 is thus simply given by :
 \begin{equation}
   \label{eq:dynspg_exp}
+  \label{eq:DYN_spg_exp}
   \left\{
     \begin{aligned}
 …
 % Split-explict free surface formulation
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Split-explicit free surface (\texttt{ln\_dynspg\_ts}\forcode{=.true.})]
+{Split-explicit free surface (\protect\np{ln\_dynspg\_ts}\forcode{=.true.})}
+\subsection[Split-explicit free surface (\texttt{ln\_dynspg\_ts}\forcode{ = .true.})]{Split-explicit free surface (\protect\np{ln\_dynspg\_ts}\forcode{ = .true.})}
 \label{subsec:DYN_spg_ts}
 %------------------------------------------namsplit-----------------------------------------------------------
 …
 The general idea is to solve the free surface equation and the associated barotropic velocity equations with
 a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables
 (\autoref{fig:DYN_dynspg_ts}).
+(\autoref{fig:DYN_spg_ts}).
 The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through
 the \np{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$.
 …
 The barotropic mode solves the following equations:
 % \begin{subequations}
 %  \label{eq:BT}
 \begin{equation}
   \label{eq:BT_dyn}
+%  \label{eq:DYN_BT}
+\begin{equation}
+  \label{eq:DYN_BT_dyn}
   \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}=
   -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h}
 …
 \end{equation}
 \[
   % \label{eq:BT_ssh}
+  % \label{eq:DYN_BT_ssh}
   \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E
 \]
 …
 where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing coupling term between modes,
 surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency.
 The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress
 (see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration.
+The third term on the right hand side of \autoref{eq:DYN_BT_dyn} represents the bottom stress
+(see section \autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration.
 Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm
 detailed in \citet{shchepetkin.mcwilliams_OM05}.
 …
     \includegraphics[width=\textwidth]{Fig_DYN_dynspg_ts}
     \caption{
       \protect\label{fig:DYN_dynspg_ts}
+      \protect\label{fig:DYN_spg_ts}
       Schematic of the split-explicit time stepping scheme for the external and internal modes.
       Time increases to the right. In this particular exemple,
 …
 In the default case (\np{ln\_bt\_fw}\forcode{=.true.}),
 the external mode is integrated between \textit{now} and \textit{after} baroclinic time-steps
 (\autoref{fig:DYN_dynspg_ts}a).
+(\autoref{fig:DYN_spg_ts}a).
 To avoid aliasing of fast barotropic motions into three dimensional equations,
 time filtering is eventually applied on barotropic quantities (\np{ln\_bt\_av}\forcode{=.true.}).
 …
 % Filtered free surface formulation
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Filtered free surface (\texttt{dynspg\_flt?})]
+{Filtered free surface (\protect\texttt{dynspg\_flt?})}
+\subsection[Filtered free surface (\texttt{dynspg\_flt?})]{Filtered free surface (\protect\texttt{dynspg\_flt?})}
 \label{subsec:DYN_spg_fltp}
 The filtered formulation follows the \citet{roullet.madec_JGR00} implementation.
 The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly.
+The extra term introduced in the equations (see \autoref{subsec:MB_free_surface}) is solved implicitly.
 The elliptic solvers available in the code are documented in \autoref{chap:MISC}.
 …
 \gmcomment{               %%% copy from chap-model basics
   \[
     % \label{eq:spg_flt}
+    % \label{eq:DYN_spg_flt}
     \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t}= {\mathrm {\mathbf M}}
     - g \nabla \left( \tilde{\rho} \ \eta \right)
 …
   $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density,
   and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient,
   non-linear and viscous terms in \autoref{eq:PE_dyn}.
+  non-linear and viscous terms in \autoref{eq:MB_dyn}.
 }   %end gmcomment
 …
 % Lateral diffusion term
 % ================================================================
+\section[Lateral diffusion term and operators (\textit{dynldf.F90})]
+{Lateral diffusion term and operators (\protect\mdl{dynldf})}
+\section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})}
 \label{sec:DYN_ldf}
 %------------------------------------------nam_dynldf----------------------------------------------------
 …
 \ie\ the velocity appearing in its expression is the \textit{before} velocity in time,
 except for the pure vertical component that appears when a tensor of rotation is used.
 This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}).
+This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}).
 At the lateral boundaries either free slip,
 …
 % ================================================================
+\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap=.true.})]
+{Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{=.true.})}
+\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap = .true.})]{Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})}
 \label{subsec:DYN_ldf_lap}
 For lateral iso-level diffusion, the discrete operator is:
 \begin{equation}
   \label{eq:dynldf_lap}
+  \label{eq:DYN_ldf_lap}
   \left\{
     \begin{aligned}
 …
 \end{equation}
 As explained in \autoref{subsec:PE_ldf},
+As explained in \autoref{subsec:MB_ldf},
 this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and
 ensures a complete separation between the vorticity and divergence parts of the momentum diffusion.
 …
 %           Rotated laplacian operator
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Rotated laplacian (\forcode{ln_dynldf_iso=.true.})]
+{Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{=.true.})}
+\subsection[Rotated laplacian (\forcode{ln_dynldf_iso = .true.})]{Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})}
 \label{subsec:DYN_ldf_iso}
 …
 The resulting discrete representation is:
 \begin{equation}
   \label{eq:dyn_ldf_iso}
+  \label{eq:DYN_ldf_iso}
   \begin{split}
     D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\
 …
 %           Iso-level bilaplacian operator
 %--------------------------------------------------------------------------------------------------------------
+\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap=.true.})]
+{Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{=.true.})}
+\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap = .true.})]{Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})}
 \label{subsec:DYN_ldf_bilap}
 The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:dynldf_lap} twice.
+The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice.
 It requires an additional assumption on boundary conditions:
 the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen,
 …
 %           Vertical diffusion term
 % ================================================================
+\section[Vertical diffusion term (\textit{dynzdf.F90})]
+{Vertical diffusion term (\protect\mdl{dynzdf})}
+\section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})}
 \label{sec:DYN_zdf}
 %----------------------------------------------namzdf------------------------------------------------------
 …
 (\np{ln\_zdfexp}\forcode{=.true.}) using a time splitting technique (\np{nn\_zdfexp} $>$ 1) or
 $(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp}\forcode{=.false.})
 (see \autoref{chap:STP}).
+(see \autoref{chap:TD}).
 Note that namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.
 The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is.
 The vertical diffusion operators given by \autoref{eq:PE_zdf} take the following semi-discrete space form:
 \[
   % \label{eq:dynzdf}
+The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form:
+\[
+  % \label{eq:DYN_zdf}
   \left\{
     \begin{aligned}
 …
 the vertical turbulent momentum fluxes,
 \begin{equation}
   \label{eq:dynzdf_sbc}
+  \label{eq:DYN_zdf_sbc}
   \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}
   = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v }
 …
 The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation
 (see \autoref{sec:ZDF_bfr})
+(see \autoref{sec:ZDF_drg})
 % ================================================================
 …
 \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).
 …
 %   Iterative limiters
 %-----------------------------------------------------------------------------------------
+\subsection[Directional limiter (\textit{wet\_dry.F90})]
+{Directional limiter (\mdl{wet\_dry})}
+\subsection[Directional limiter (\textit{wet\_dry.F90})]{Directional limiter (\mdl{wet\_dry})}
 \label{subsec:DYN_wd_directional_limiter}
 The principal idea of the directional limiter is that
 water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than \np{rn\_wdmin1}).
 …
 %-----------------------------------------------------------------------------------------
+\subsection[Iterative limiter (\textit{wet\_dry.F90})]
+{Iterative limiter (\mdl{wet\_dry})}
+\subsection[Iterative limiter (\textit{wet\_dry.F90})]{Iterative limiter (\mdl{wet\_dry})}
 \label{subsec:DYN_wd_iterative_limiter}
+\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]
+{Iterative flux limiter (\mdl{wet\_dry})}
+\label{subsubsec:DYN_wd_il_spg_limiter}
+\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]{Iterative flux limiter (\mdl{wet\_dry})}
+\label{subsec:DYN_wd_il_spg_limiter}
 The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry''
 …
 The continuity equation for the total water depth in a column
+\begin{equation} \label{dyn_wd_continuity}
+ \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .
+\begin{equation}
+  \label{eq:DYN_wd_continuity}
+  \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .
 \end{equation}
 can be written in discrete form  as
+\begin{align} \label{dyn_wd_continuity_2}
+\frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) )
+&= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\
+&= \mathrm{zzflx}_{i,j} .
+\begin{align}
+  \label{eq:DYN_wd_continuity_2}
+  \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) )
+  &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\
+  &= \mathrm{zzflx}_{i,j} .
 \end{align}
 …
 (zzflxp) and fluxes that are into the cell (zzflxn).  Clearly
+\begin{equation} \label{dyn_wd_zzflx_p_n_1}
+\mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .
+\begin{equation}
+  \label{eq:DYN_wd_zzflx_p_n_1}
+  \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .
 \end{equation}
 …
 $\mathrm{zcoef}_{i,j}^{(m)}$ such that:
+\begin{equation} \label{dyn_wd_continuity_coef}
+\begin{split}
+\mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\
+\mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}
+\end{split}
+\begin{equation}
+  \label{eq:DYN_wd_continuity_coef}
+  \begin{split}
+    \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\
+    \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}
+  \end{split}
 \end{equation}
 …
 The iteration is initialised by setting
+\begin{equation} \label{dyn_wd_zzflx_initial}
+\mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} .
+\begin{equation}
+  \label{eq:DYN_wd_zzflx_initial}
+  \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} .
 \end{equation}
 The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the
 cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell
 times the timestep divided by the cell area. Using (\ref{dyn_wd_continuity_2}) this
+times the timestep divided by the cell area. Using (\autoref{eq:DYN_wd_continuity_2}) this
 condition is
+\begin{equation} \label{dyn_wd_continuity_if}
+h_{i,j}(t_e)  - \mathrm{rn\_wdmin1} <  \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) .
+\end{equation}
+Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum
+\begin{equation}
+  \label{eq:DYN_wd_continuity_if}
+  h_{i,j}(t_e)  - \mathrm{rn\_wdmin1} <  \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) .
+\end{equation}
+Rearranging (\autoref{eq:DYN_wd_continuity_if}) we can obtain an expression for the maximum
 outward flux that can be allowed and still maintain the minimum wet depth:
+\begin{equation} \label{dyn_wd_max_flux}
+\begin{split}
+\mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
+\phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big]
+\end{split}
+\begin{equation}
+  \label{eq:DYN_wd_max_flux}
+  \begin{split}
+    \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
+    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big]
+  \end{split}
 \end{equation}
 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is
 this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an
+this necessary/desirable?]}. Substituting from (\autoref{eq:DYN_wd_continuity_coef}) gives an
 expression for the coefficient needed to multiply the outward flux at this cell in order
 to avoid drying.
+\begin{equation} \label{dyn_wd_continuity_nxtcoef}
+\begin{split}
+\mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
+\phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} }
+\end{split}
+\begin{equation}
+  \label{eq:DYN_wd_continuity_nxtcoef}
+  \begin{split}
+    \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
+    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} }
+  \end{split}
 \end{equation}
 …
 %      Surface pressure gradients
 %----------------------------------------------------------------------------------------
+\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]
+{Modification of surface pressure gradients (\mdl{dynhpg})}
+\label{subsubsec:DYN_wd_il_spg}
+\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]{Modification of surface pressure gradients (\mdl{dynhpg})}
+\label{subsec:DYN_wd_il_spg}
 At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the
 …
 neighbouring $(i+1,j)$ and $(i,j)$ tracer points.  zcpx is calculated using two logicals
 variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid
 column.  The three possible combinations are illustrated in figure \ref{Fig_WAD_dynhpg}.
+column.  The three possible combinations are illustrated in figure \autoref{fig:DYN_WAD_dynhpg}.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht] \begin{center}
 \includegraphics[width=\textwidth]{Fig_WAD_dynhpg}
+\caption{ \label{Fig_WAD_dynhpg}
+Illustrations of the three possible combinations of the logical variables controlling the
+limiting of the horizontal pressure gradient in wetting and drying regimes}
+\end{center}\end{figure}
+\caption{
+  \label{fig:DYN_WAD_dynhpg}
+  Illustrations of the three possible combinations of the logical variables controlling the
+  limiting of the horizontal pressure gradient in wetting and drying regimes}
+\end{center}
+\end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 of the topography at the two points:
+\begin{equation} \label{dyn_ll_tmp1}
+\begin{split}
+\mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\
+\begin{equation}
+  \label{eq:DYN_ll_tmp1}
+  \begin{split}
+    \mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\
                      & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\  .and.} \\
 & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\
 & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\
 & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }
 \end{split}
+                     & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\
+                     & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\
+                     & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }
+  \end{split}
 \end{equation}
 …
 at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$
+\begin{equation} \label{dyn_ll_tmp2}
+\begin{split}
+\mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\
+& \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\
+& \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .
+\end{split}
+\begin{equation}
+  \label{eq:DYN_ll_tmp2}
+  \begin{split}
+    \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\
+    & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\
+    & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .
+  \end{split}
 \end{equation}
 …
 conditions.
+\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]
+{Additional considerations (\mdl{usrdef\_zgr})}
+\label{subsubsec:WAD_additional}
+\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]{Additional considerations (\mdl{usrdef\_zgr})}
+\label{subsec:DYN_WAD_additional}
 In the very shallow water where wetting and drying occurs the parametrisation of
 …
 %      The WAD test cases
 %----------------------------------------------------------------------------------------
+\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]
+{The WAD test cases (\mdl{usrdef\_zgr})}
+\label{WAD_test_cases}
+\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]{The WAD test cases (\mdl{usrdef\_zgr})}
+\label{subsec:DYN_WAD_test_cases}
 See the WAD tests MY\_DOC documention for details of the WAD test cases.
 …
 % Time evolution term
 % ================================================================
+\section[Time evolution term (\textit{dynnxt.F90})]
+{Time evolution term (\protect\mdl{dynnxt})}
+\section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})}
 \label{sec:DYN_nxt}
 …
 Options are defined through the \nam{dom} namelist variables.
 The general framework for dynamics time stepping is a leap-frog scheme,
 \ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:STP}).
+\ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:TD}).
 The scheme is applied to the velocity, except when
 using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux})
 in the variable volume case (\texttt{vvl?} defined),
 where it has to be applied to the thickness weighted velocity (see \autoref{sec:A_momentum})
+where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum})
 $\bullet$ vector invariant form or linear free surface
 (\np{ln\_dynhpg\_vec}\forcode{=.true.} ; \texttt{vvl?} not defined):
 \[
   % \label{eq:dynnxt_vec}
+  % \label{eq:DYN_nxt_vec}
   \left\{
     \begin{aligned}
 …
 (\np{ln\_dynhpg\_vec}\forcode{=.false.} ; \texttt{vvl?} defined):
 \[
   % \label{eq:dynnxt_flux}
+  % \label{eq:DYN_nxt_flux}
   \left\{
     \begin{aligned}

NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex

-                      r11537
+                      r11543
 \[
   % \label{eq:lbc_aaaa}
+  % \label{eq:LBC_aaaa}
   \frac{A^{lT} }{e_1 }\frac{\partial T}{\partial i}\equiv \frac{A_u^{lT}
   }{e_{1u} } \; \delta_{i+1 / 2} \left[ T \right]\;\;mask_u
 …
   the no-slip boundary condition, simply by multiplying it by the mask$_{f}$ :
   \[
     % \label{eq:lbc_bbbb}
+    % \label{eq:LBC_bbbb}
     \zeta \equiv \frac{1}{e_{1f} {\kern 1pt}e_{2f} }\left( {\delta_{i+1/2}
         \left[ {e_{2v} \,v} \right]-\delta_{j+1/2} \left[ {e_{1u} \,u} \right]}
 …
 The north fold boundary condition has been introduced in order to handle the north boundary of
 a three-polar ORCA grid.
 Such a grid has two poles in the northern hemisphere (\autoref{fig:MISC_ORCA_msh},
 and thus requires a specific treatment illustrated in \autoref{fig:North_Fold_T}.
+Such a grid has two poles in the northern hemisphere (\autoref{fig:CFGS_ORCA_msh},
+and thus requires a specific treatment illustrated in \autoref{fig:LBC_North_Fold_T}.
 Further information can be found in \mdl{lbcnfd} module which applies the north fold boundary condition.
 …
     \includegraphics[width=\textwidth]{Fig_North_Fold_T}
     \caption{
       \protect\label{fig:North_Fold_T}
+      \protect\label{fig:LBC_North_Fold_T}
       North fold boundary with a $T$-point pivot and cyclic east-west boundary condition ($jperio=4$),
       as used in ORCA 2, 1/4, and 1/12.
 …
 %-----------------------------------------------------------------------------------------------
+For massively parallel processing (mpp), a domain decomposition method is used. The basic idea of the method is to split the large computation domain of a numerical experiment into several smaller domains and solve the set of equations by addressing independent local problems. Each processor has its own local memory and computes the model equation over a subdomain of the whole model domain. The subdomain boundary conditions are specified through communications between processors which are organized by explicit statements (message passing method). The present implementation is largely inspired by Guyon's work [Guyon 1995].
+For massively parallel processing (mpp), a domain decomposition method is used.
+The basic idea of the method is to split the large computation domain of a numerical experiment into several smaller domains and
+solve the set of equations by addressing independent local problems.
+Each processor has its own local memory and computes the model equation over a subdomain of the whole model domain.
+The subdomain boundary conditions are specified through communications between processors which are organized by
+explicit statements (message passing method).
+The present implementation is largely inspired by Guyon's work [Guyon 1995].
 The parallelization strategy is defined by the physical characteristics of the ocean model.
 …
 each processor sends to its neighbouring processors the update values of the points corresponding to
 the interior overlapping area to its neighbouring sub-domain (\ie\ the innermost of the two overlapping rows).
+Communications are first done according to the east-west direction and next according to the north-south direction. There is no specific communications for the corners. The communication is done through the Message Passing Interface (MPI) and requires \key{mpp\_mpi}. Use also \key{mpi2} if MPI3 is not available on your computer.
+Communications are first done according to the east-west direction and next according to the north-south direction.
+There is no specific communications for the corners.
+The communication is done through the Message Passing Interface (MPI) and requires \key{mpp\_mpi}.
+Use also \key{mpi2} if MPI3 is not available on your computer.
 The data exchanges between processors are required at the very place where
 lateral domain boundary conditions are set in the mono-domain computation:
 …
     \includegraphics[width=\textwidth]{Fig_mpp}
     \caption{
       \protect\label{fig:mpp}
+      \protect\label{fig:LBC_mpp}
       Positioning of a sub-domain when massively parallel processing is used.
+    }
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+In \NEMO, the splitting is regular and arithmetic. The total number of subdomains corresponds to the number of MPI processes allocated to \NEMO\ when the model is launched (\ie\ mpirun -np x ./nemo will automatically give x subdomains). The i-axis is divided by \np{jpni} and the j-axis by \np{jpnj}. These parameters are defined in \nam{mpp} namelist. If \np{jpni} and \np{jpnj} are < 1, they will be automatically redefined in the code to give the best domain decomposition (see bellow).
+Each processor is independent and without message passing or synchronous process, programs run alone and access just its own local memory. For this reason, the main model dimensions are now the local dimensions of the subdomain (pencil) that are named \jp{jpi}, \jp{jpj}, \jp{jpk}.
+These dimensions include the internal domain and the overlapping rows. The number of rows to exchange (known as the halo) is usually set to one (nn\_hls=1, in \mdl{par\_oce}, and must be kept to one until further notice). The whole domain dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}. The relationship between the whole domain and a sub-domain is:
+\[
+  jpi = ( jpiglo-2\times nn\_hls + (jpni-1) ) / jpni + 2\times nn\_hls
+\]
+\[
+In \NEMO, the splitting is regular and arithmetic.
+The total number of subdomains corresponds to the number of MPI processes allocated to \NEMO\ when the model is launched
+(\ie\ mpirun -np x ./nemo will automatically give x subdomains).
+The i-axis is divided by \np{jpni} and the j-axis by \np{jpnj}.
+These parameters are defined in \nam{mpp} namelist.
+If \np{jpni} and \np{jpnj} are < 1, they will be automatically redefined in the code to give the best domain decomposition
+(see bellow).
+Each processor is independent and without message passing or synchronous process, programs run alone and access just its own local memory.
+For this reason,
+the main model dimensions are now the local dimensions of the subdomain (pencil) that are named \jp{jpi}, \jp{jpj}, \jp{jpk}.
+These dimensions include the internal domain and the overlapping rows.
+The number of rows to exchange (known as the halo) is usually set to one (nn\_hls=1, in \mdl{par\_oce},
+and must be kept to one until further notice).
+The whole domain dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}.
+The relationship between the whole domain and a sub-domain is:
+\begin{gather*}
+  jpi = ( jpiglo-2\times nn\_hls + (jpni-1) ) / jpni + 2\times nn\_hls \\
   jpj = ( jpjglo-2\times nn\_hls + (jpnj-1) ) / jpnj + 2\times nn\_hls
 \]
 One also defines variables nldi and nlei which correspond to the internal domain bounds, and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain (\autoref{fig:mpp}). Note that since the version 4, there is no more extra-halo area as defined in \autoref{fig:mpp} so \jp{jpi} is now always equal to nlci and \jp{jpj} equal to nlcj.
+\end{gather*}
+One also defines variables nldi and nlei which correspond to the internal domain bounds, and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain (\autoref{fig:LBC_mpp}). Note that since the version 4, there is no more extra-halo area as defined in \autoref{fig:LBC_mpp} so \jp{jpi} is now always equal to nlci and \jp{jpj} equal to nlcj.
 An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$,
 a global array (whole domain) by the relationship:
 \[
   % \label{eq:lbc_nimpp}
+  % \label{eq:LBC_nimpp}
   T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k),
 \]
 …
 If the domain decomposition is automatically defined (when \np{jpni} and \np{jpnj} are < 1), the decomposition chosen by the model will minimise the sub-domain size (defined as $max_{all domains}(jpi \times jpj)$) and maximize the number of eliminated land subdomains. This means that no other domain decomposition (a set of \np{jpni} and \np{jpnj} values) will use less processes than $(jpni  \times  jpnj - N_{land})$ and get a smaller subdomain size.
 In order to specify $N_{mpi}$ properly (minimize $N_{useless}$), you must run the model once with \np{ln\_list} activated. In this case, the model will start the initialisation phase, print the list of optimum decompositions ($N_{mpi}$, \np{jpni} and \np{jpnj}) in \texttt{ocean.output} and directly abort. The maximum value of $N_{mpi}$ tested in this list is given by $max(N_{MPI\_tasks}, \np{jpni} \times \np{jpnj})$. For example, run the model on 40 nodes with ln\_list activated and $\np{jpni} = 10000$ and $\np{jpnj} = 1$, will print the list of optimum domains decomposition from 1 to about 10000.
 Processors are numbered from 0 to $N_{mpi} - 1$. Subdomains containning some ocean points are numbered first from 0 to $jpni * jpnj - N_{land} -1$. The remaining $N_{useless}$ land subdomains are numbered next, which means that, for a given (\np{jpni}, \np{jpnj}), the numbers attributed to he ocean subdomains do not vary with $N_{useless}$.
+In order to specify $N_{mpi}$ properly (minimize $N_{useless}$), you must run the model once with \np{ln\_list} activated. In this case, the model will start the initialisation phase, print the list of optimum decompositions ($N_{mpi}$, \np{jpni} and \np{jpnj}) in \texttt{ocean.output} and directly abort. The maximum value of $N_{mpi}$ tested in this list is given by $max(N_{MPI\_tasks}, \np{jpni} \times \np{jpnj})$. For example, run the model on 40 nodes with ln\_list activated and $\np{jpni} = 10000$ and $\np{jpnj} = 1$, will print the list of optimum domains decomposition from 1 to about 10000.
+Processors are numbered from 0 to $N_{mpi} - 1$. Subdomains containning some ocean points are numbered first from 0 to $jpni * jpnj - N_{land} -1$. The remaining $N_{useless}$ land subdomains are numbered next, which means that, for a given (\np{jpni}, \np{jpnj}), the numbers attributed to he ocean subdomains do not vary with $N_{useless}$.
 When land processors are eliminated, the value corresponding to these locations in the model output files is undefined. \np{ln\_mskland} must be activated in order avoid Not a Number values in output files. Note that it is better to not eliminate land processors when creating a meshmask file (\ie\ when setting a non-zero value to \np{nn\_msh}).
 …
   \begin{center}
     \includegraphics[width=\textwidth]{Fig_mppini2}
     \caption {
       \protect\label{fig:mppini2}
+    \caption[Atlantic domain]{
+      \protect\label{fig:LBC_mppini2}
       Example of Atlantic domain defined for the CLIPPER projet.
       Initial grid is composed of 773 x 1236 horizontal points.
 …
 %----------------------------------------------
 \subsection{Namelists}
 \label{subsec:BDY_namelist}
+\label{subsec:LBC_bdy_namelist}
 The BDY module is activated by setting \np{ln\_bdy}\forcode{=.true.} .
 …
 In the example above, there are two boundary sets, the first of which is defined via a file and
 the second is defined in the namelist.
 For more details of the definition of the boundary geometry see section \autoref{subsec:BDY_geometry}.
+For more details of the definition of the boundary geometry see section \autoref{subsec:LBC_bdy_geometry}.
 For each boundary set a boundary condition has to be chosen for the barotropic solution
 …
 %----------------------------------------------
 \subsection{Flow relaxation scheme}
 \label{subsec:BDY_FRS_scheme}
+\label{subsec:LBC_bdy_FRS_scheme}
 The Flow Relaxation Scheme (FRS) \citep{davies_QJRMS76,engedahl_T95},
 …
 Given a model prognostic variable $\Phi$
 \[
   % \label{eq:bdy_frs1}
+  % \label{eq:LBC_bdy_frs1}
   \Phi(d) = \alpha(d)\Phi_{e}(d) + (1-\alpha(d))\Phi_{m}(d)\;\;\;\;\; d=1,N
 \]
 …
 the prognostic equation for $\Phi$ of the form:
 \[
   % \label{eq:bdy_frs2}
+  % \label{eq:LBC_bdy_frs2}
   -\frac{1}{\tau}\left(\Phi - \Phi_{e}\right)
 \]
 where the relaxation time scale $\tau$ is given by a function of $\alpha$ and the model time step $\Delta t$:
 \[
   % \label{eq:bdy_frs3}
+  % \label{eq:LBC_bdy_frs3}
   \tau = \frac{1-\alpha}{\alpha}  \,\rdt
 \]
 …
 The function $\alpha$ is specified as a $tanh$ function:
 \[
   % \label{eq:bdy_frs4}
+  % \label{eq:LBC_bdy_frs4}
   \alpha(d) = 1 - \tanh\left(\frac{d-1}{2}\right),       \quad d=1,N
 \]
 …
 %----------------------------------------------
 \subsection{Flather radiation scheme}
 \label{subsec:BDY_flather_scheme}
+\label{subsec:LBC_bdy_flather_scheme}
 The \citet{flather_JPO94} scheme is a radiation condition on the normal,
 depth-mean transport across the open boundary.
 It takes the form
+\begin{equation}  \label{eq:bdy_fla1}
+U = U_{e} + \frac{c}{h}\left(\eta - \eta_{e}\right),
+\begin{equation}
+  \label{eq:LBC_bdy_fla1}
+  U = U_{e} + \frac{c}{h}\left(\eta - \eta_{e}\right),
 \end{equation}
 where $U$ is the depth-mean velocity normal to the boundary and $\eta$ is the sea surface height,
 …
 the external depth-mean normal velocity,
 plus a correction term that allows gravity waves generated internally to exit the model boundary.
 Note that the sea-surface height gradient in \autoref{eq:bdy_fla1} is a spatial gradient across the model boundary,
+Note that the sea-surface height gradient in \autoref{eq:LBC_bdy_fla1} is a spatial gradient across the model boundary,
 so that $\eta_{e}$ is defined on the $T$ points with $nbr=1$ and $\eta$ is defined on the $T$ points with $nbr=2$.
 $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:BDY_orlanski_scheme}
+\label{subsec:LBC_bdy_orlanski_scheme}
 The Orlanski scheme is based on the algorithm described by \citep{marchesiello.mcwilliams.ea_OM01}, hereafter MMS.
 …
 The adaptive Orlanski condition solves a wave plus relaxation equation at the boundary:
 \begin{equation}
+\frac{\partial\phi}{\partial t} + c_x \frac{\partial\phi}{\partial x} + c_y \frac{\partial\phi}{\partial y} =
                                                 -\frac{1}{\tau}(\phi - \phi^{ext})
+\label{eq:wave_continuous}
+  \label{eq:LBC_wave_continuous}
+  \frac{\partial\phi}{\partial t} + c_x \frac{\partial\phi}{\partial x} + c_y \frac{\partial\phi}{\partial y} =
+  -\frac{1}{\tau}(\phi - \phi^{ext})
 \end{equation}
 …
 velocities are diagnosed from the model fields as:
+\begin{equation} \label{eq:cx}
+c_x = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial x}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2}
+\begin{equation}
+  \label{eq:LBC_cx}
+  c_x = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial x}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2}
 \end{equation}
 \begin{equation}
 \label{eq:cy}
 c_y = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial y}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2}
+  \label{eq:LBC_cy}
+  c_y = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial y}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2}
 \end{equation}
 (As noted by MMS, this is a circular diagnosis of the phase speeds which only makes sense on a discrete grid).
 Equation (\autoref{eq:wave_continuous}) is defined adaptively depending on the sign of the phase velocity normal to the boundary $c_x$.
+Equation (\autoref{eq:LBC_wave_continuous}) is defined adaptively depending on the sign of the phase velocity normal to the boundary $c_x$.
 For $c_x$ outward, we have
 …
 \begin{equation}
+\tau = \tau_{in}\,\,\,;\,\,\, c_x = c_y = 0
+\label{eq:tau_in}
+  \label{eq:LBC_tau_in}
+  \tau = \tau_{in}\,\,\,;\,\,\, c_x = c_y = 0
 \end{equation}
 Generally the relaxation time scale at inward propagation points (\np{rn\_time\_dmp}) is set much shorter than the time scale at outward propagation
 points (\np{rn\_time\_dmp\_out}) so that the solution is constrained more strongly by the external data at inward propagation points.
 See \autoref{subsec:BDY_relaxation} for detailed on the spatial shape of the scaling.\\
+See \autoref{subsec:LBC_bdy_relaxation} for detailed on the spatial shape of the scaling.\\
 The ``normal propagation of oblique radiation'' or NPO approximation (called \forcode{'orlanski_npo'}) involves assuming
 that $c_y$ is zero in equation (\autoref{eq:wave_continuous}), but including
 this term in the denominator of equation (\autoref{eq:cx}). Both versions of the scheme are options in BDY. Equations
 (\autoref{eq:wave_continuous}) - (\autoref{eq:tau_in}) correspond to equations (13) - (15) and (2) - (3) in MMS.\\
+that $c_y$ is zero in equation (\autoref{eq:LBC_wave_continuous}), but including
+this term in the denominator of equation (\autoref{eq:LBC_cx}). Both versions of the scheme are options in BDY. Equations
+(\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:BDY_relaxation}
+\label{subsec:LBC_bdy_relaxation}
 In addition to a specific boundary condition specified as \np{cn\_tra} and \np{cn\_dyn3d}, relaxation on baroclinic velocities and tracers variables are available.
 …
 %----------------------------------------------
 \subsection{Boundary geometry}
 \label{subsec:BDY_geometry}
+\label{subsec:LBC_bdy_geometry}
 Each open boundary set is defined as a list of points.
 …
 %----------------------------------------------
 \subsection{Input boundary data files}
 \label{subsec:BDY_data}
+\label{subsec:LBC_bdy_data}
 The data files contain the data arrays in the order in which the points are defined in the $nbi$ and $nbj$ arrays.
 …
 %----------------------------------------------
 \subsection{Volume correction}
 \label{subsec:BDY_vol_corr}
+\label{subsec:LBC_bdy_vol_corr}
 There is an option to force the total volume in the regional model to be constant.
 …
 %----------------------------------------------
 \subsection{Tidal harmonic forcing}
 \label{subsec:BDY_tides}
+\label{subsec:LBC_bdy_tides}
 %-----------------------------------------nambdy_tide--------------------------------------------

NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex

-                      r11537
+                      r11543
 \newpage
 The lateral physics terms in the momentum and tracer equations have been described in \autoref{eq:PE_zdf} and
+The lateral physics terms in the momentum and tracer equations have been described in \autoref{eq:MB_zdf} and
 their discrete formulation in \autoref{sec:TRA_ldf} and \autoref{sec:DYN_ldf}).
 In this section we further discuss each lateral physics option.
 …
 Note that this chapter describes the standard implementation of iso-neutral tracer mixing.
 Griffies's implementation, which is used if \np{ln\_traldf\_triad}\forcode{=.true.},
 is described in \autoref{apdx:triad}
+is described in \autoref{apdx:TRIADS}
 %-----------------------------------namtra_ldf - namdyn_ldf--------------------------------------------
 …
 the cell of the quantity to be diffused.
 For a tracer, this leads to the following four slopes:
 $r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \autoref{eq:tra_ldf_iso}),
+$r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \autoref{eq:TRA_ldf_iso}),
 while for momentum the slopes are  $r_{1t}$, $r_{1uw}$, $r_{2f}$, $r_{2uw}$ for $u$ and
 $r_{1f}$, $r_{1vw}$, $r_{2t}$, $r_{2vw}$ for $v$.
 …
 In $s$-coordinates, geopotential mixing (\ie\ horizontal mixing) $r_1$ and $r_2$ are the slopes between
 the geopotential and computational surfaces.
 Their discrete formulation is found by locally solving \autoref{eq:tra_ldf_iso} when
+Their discrete formulation is found by locally solving \autoref{eq:TRA_ldf_iso} when
 the diffusive fluxes in the three directions are set to zero and $T$ is assumed to be horizontally uniform,
 \ie\ a linear function of $z_T$, the depth of a $T$-point.
 …
 \begin{equation}
   \label{eq:ldfslp_geo}
+  \label{eq:LDF_slp_geo}
   \begin{aligned}
     r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)}
 …
 Their discrete formulation is found using the fact that the diffusive fluxes of
 locally referenced potential density (\ie\ $in situ$ density) vanish.
 So, substituting $T$ by $\rho$ in \autoref{eq:tra_ldf_iso} and setting the diffusive fluxes in
+So, substituting $T$ by $\rho$ in \autoref{eq:TRA_ldf_iso} and setting the diffusive fluxes in
 the three directions to zero leads to the following definition for the neutral slopes:
 \begin{equation}
   \label{eq:ldfslp_iso}
+  \label{eq:LDF_slp_iso}
   \begin{split}
     r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac{\delta_{i+1/2}[\rho]}
 …
 %gm% rewrite this as the explanation is not very clear !!!
 %In practice, \autoref{eq:ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth.
 %By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
 %In the $z$-coordinate, the derivative of the  \autoref{eq:ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation.
 As the mixing is performed along neutral surfaces, the gradient of $\rho$ in \autoref{eq:ldfslp_iso} has to
+%In practice, \autoref{eq:LDF_slp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:LDF_slp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth.
+%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{eq:LDF_slp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
+%In the $z$-coordinate, the derivative of the  \autoref{eq:LDF_slp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation.
+As the mixing is performed along neutral surfaces, the gradient of $\rho$ in \autoref{eq:LDF_slp_iso} has to
 be evaluated at the same local pressure
 (which, in decibars, is approximated by the depth in meters in the model).
 Therefore \autoref{eq:ldfslp_iso} cannot be used as such,
+Therefore \autoref{eq:LDF_slp_iso} cannot be used as such,
 but further transformation is needed depending on the vertical coordinate used:
 …
 \item[$z$-coordinate with full step: ]
   in \autoref{eq:ldfslp_iso} the densities appearing in the $i$ and $j$ derivatives  are taken at the same depth,
+  in \autoref{eq:LDF_slp_iso} the densities appearing in the $i$ and $j$ derivatives  are taken at the same depth,
   thus the $in situ$ density can be used.
   This is not the case for the vertical derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$,
 …
   in the current release of \NEMO, iso-neutral mixing is only employed for $s$-coordinates if
   the Griffies scheme is used (\np{ln\_traldf\_triad}\forcode{=.true.};
   see \autoref{apdx:triad}).
+  see \autoref{apdx:TRIADS}).
   In other words, iso-neutral mixing will only be accurately represented with a linear equation of state
   (\np{ln\_seos}\forcode{=.true.}).
   In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:ldfslp_iso}
+  In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:LDF_slp_iso}
   will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes.
 …
 \[
   % \label{eq:ldfslp_iso2}
+  % \label{eq:LDF_slp_iso2}
   \begin{split}
     r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac
 …
 To overcome this problem, several techniques have been proposed in which the numerical schemes of
 the ocean model are modified \citep{weaver.eby_JPO97, griffies.gnanadesikan.ea_JPO98}.
 Griffies's scheme is now available in \NEMO\ if \np{ln\_traldf\_triad}\forcode{=.true.}; see \autoref{apdx:triad}.
+Griffies's scheme is now available in \NEMO\ if \np{ln\_traldf\_triad}\forcode{ = .true.}; see \autoref{apdx:TRIADS}.
 Here, another strategy is presented \citep{lazar_phd97}:
 a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of
 …
     \includegraphics[width=\textwidth]{Fig_eiv_slp}
     \caption{
       \protect\label{fig:eiv_slp}
+      \protect\label{fig:LDF_eiv_slp}
       Vertical profile of the slope used for lateral mixing in the mixed layer:
       \textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior,
 …
 The iso-neutral diffusion operator on momentum is the same as the one used on tracers but
 applied to each component of the velocity separately
 (see \autoref{eq:dyn_ldf_iso} in section~\autoref{subsec:DYN_ldf_iso}).
+(see \autoref{eq:DYN_ldf_iso} in section~\autoref{subsec:DYN_ldf_iso}).
 The slopes between the surface along which the diffusion operator acts and the surface of computation
 ($z$- or $s$-surfaces) are defined at $T$-, $f$-, and \textit{uw}- points for the $u$-component, and $T$-, $f$- and
 \textit{vw}- points for the $v$-component.
 They are computed from the slopes used for tracer diffusion,
 \ie\ \autoref{eq:ldfslp_geo} and \autoref{eq:ldfslp_iso}:
+\ie\ \autoref{eq:LDF_slp_geo} and \autoref{eq:LDF_slp_iso}:
 \[
   % \label{eq:ldfslp_dyn}
+  % \label{eq:LDF_slp_dyn}
   \begin{aligned}
     &r_{1t}\ \ = \overline{r_{1u}}^{\,i}       &&&    r_{1f}\ \ &= \overline{r_{1u}}^{\,i+1/2} \\
 …
 \begin{equation}
   \label{eq:constantah}
+  \label{eq:LDF_constantah}
   A_o^l = \left\{
     \begin{aligned}
 …
 In the vertically varying case, a hyperbolic variation of the lateral mixing coefficient is introduced in which
 the surface value is given by \autoref{eq:constantah}, the bottom value is 1/4 of the surface value,
+the surface value is given by \autoref{eq:LDF_constantah}, the bottom value is 1/4 of the surface value,
 and the transition takes place around z=500~m with a width of 200~m.
 This profile is hard coded in module \mdl{ldfc1d\_c2d}, but can be easily modified by users.
 …
 the type of operator used:
 \begin{equation}
   \label{eq:title}
+  \label{eq:LDF_title}
   A_l = \left\{
     \begin{aligned}
 …
 model configurations presenting large changes in grid spacing such as global ocean models.
 Indeed, in such a case, a constant mixing coefficient can lead to a blow up of the model due to
 large coefficient compare to the smallest grid size (see \autoref{sec:STP_forward_imp}),
+large coefficient compare to the smallest grid size (see \autoref{sec:TD_forward_imp}),
 especially when using a bilaplacian operator.
 …
 \begin{equation}
   \label{eq:flowah}
+  \label{eq:LDF_flowah}
   A_l = \left\{
     \begin{aligned}
 …
 \begin{equation}
   \label{eq:smag1}
+  \label{eq:LDF_smag1}
   \begin{split}
     T_{smag}^{-1} & = \sqrt{\left( \partial_x u - \partial_y v\right)^2 + \left( \partial_y u + \partial_x v\right)^2  } \\
 …
 \begin{equation}
   \label{eq:smag2}
+  \label{eq:LDF_smag2}
   A_{smag} = \left\{
     \begin{aligned}
 …
 \end{equation}
 For stability reasons, upper and lower limits are applied on the resulting coefficient (see \autoref{sec:STP_forward_imp}) so that:
 \begin{equation}
   \label{eq:smag3}
+For stability reasons, upper and lower limits are applied on the resulting coefficient (see \autoref{sec:TD_forward_imp}) so that:
+\begin{equation}
+  \label{eq:LDF_smag3}
     \begin{aligned}
       & C_{min} \frac{1}{2}   \lvert U \rvert  e    < A_{smag} < C_{max} \frac{e^2}{   8\rdt}                 & \text{for laplacian operator } \\
 …
 (1) the momentum diffusion operator acting along model level surfaces is written in terms of curl and
 divergent components of the horizontal current (see \autoref{subsec:PE_ldf}).
+divergent components of the horizontal current (see \autoref{subsec:MB_ldf}).
 Although the eddy coefficient could be set to different values in these two terms,
 this option is not currently available.
 …
 (2) with an horizontally varying viscosity, the quadratic integral constraints on enstrophy and on the square of
 the horizontal divergence for operators acting along model-surfaces are no longer satisfied
 (\autoref{sec:dynldf_properties}).
+(\autoref{sec:INVARIANTS_dynldf_properties}).
 % ================================================================
 …
 the formulation of which depends on the slopes of iso-neutral surfaces.
 Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces,
 \ie\ \autoref{eq:ldfslp_geo} is used in $z$-coordinates,
 and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $s$-coordinates.
+\ie\ \autoref{eq:LDF_slp_geo} is used in $z$-coordinates,
+and the sum \autoref{eq:LDF_slp_geo} + \autoref{eq:LDF_slp_iso} in $s$-coordinates.
 If isopycnal mixing is used in the standard way, \ie\ \np{ln\_traldf\_triad}\forcode{=.false.}, the eddy induced velocity is given by:
 \begin{equation}
   \label{eq:ldfeiv}
+  \label{eq:LDF_eiv}
   \begin{split}
     u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\
 …
 \colorbox{yellow}{CASE \np{nn\_aei\_ijk\_t} = 21 to be added}
 In case of setting \np{ln\_traldf\_triad}\forcode{=.true.}, a skew form of the eddy induced advective fluxes is used, which is described in \autoref{apdx:triad}.
+In case of setting \np{ln\_traldf\_triad}\forcode{ = .true.}, a skew form of the eddy induced advective fluxes is used, which is described in \autoref{apdx:TRIADS}.
 % ================================================================

NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex

-                      r11435
+                      r11543
 Examples of the weights calculated for an observation with rectangular and radial footprints are shown in
 \autoref{fig:obsavgrec} and~\autoref{fig:obsavgrad}.
+\autoref{fig:OBS_avgrec} and~\autoref{fig:OBS_avgrad}.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
     \includegraphics[width=\textwidth]{Fig_OBS_avg_rec}
     \caption{
       \protect\label{fig:obsavgrec}
+      \protect\label{fig:OBS_avgrec}
       Weights associated with each model grid box (blue lines and numbers)
       for an observation at -170.5\deg{E}, 56.0\deg{N} with a rectangular footprint of 1\deg x 1\deg.
 …
     \includegraphics[width=\textwidth]{Fig_OBS_avg_rad}
     \caption{
       \protect\label{fig:obsavgrad}
+      \protect\label{fig:OBS_avgrad}
       Weights associated with each model grid box (blue lines and numbers)
       for an observation at -170.5\deg{E}, 56.0\deg{N} with a radial footprint with diameter 1\deg.
 …
           ({\phi_{}}_{\mathrm D}  \;  - \; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\
   \end{array}
   % \label{eq:cross}
+  % \label{eq:OBS_cross}
 \end{align*}
 point in the opposite direction to the unit normal $\widehat{\mathbf k}$
 …
     \includegraphics[width=\textwidth]{Fig_ASM_obsdist_local}
     \caption{
       \protect\label{fig:obslocal}
+      \protect\label{fig:OBS_local}
       Example of the distribution of observations with the geographical distribution of observational data.
+    }
 …
 This is the simplest option in which the observations are distributed according to
 the domain of the grid-point parallelization.
 \autoref{fig:obslocal} shows an example of the distribution of the {\em in situ} data on processors with
+\autoref{fig:OBS_local} shows an example of the distribution of the {\em in situ} data on processors with
 a different colour for each observation on a given processor for a 4 $\times$ 2 decomposition with ORCA2.
 The grid-point domain decomposition is clearly visible on the plot.
 …
     \includegraphics[width=\textwidth]{Fig_ASM_obsdist_global}
     \caption{
       \protect\label{fig:obsglobal}
+      \protect\label{fig:OBS_global}
       Example of the distribution of observations with the round-robin distribution of observational data.
+    }
 …
 use message passing in order to retrieve the stencil for interpolation.
 The simplest distribution of the observations is to distribute them using a round-robin scheme.
 \autoref{fig:obsglobal} shows the distribution of the {\em in situ} data on processors for
+\autoref{fig:OBS_global} shows the distribution of the {\em in situ} data on processors for
 the round-robin distribution of observations with a different colour for each observation on a given processor for
 a 4 $\times$ 2 decomposition with ORCA2 for the same input data as in \autoref{fig:obslocal}.
+a 4 $\times$ 2 decomposition with ORCA2 for the same input data as in \autoref{fig:OBS_local}.
 The observations are now clearly randomly distributed on the globe.
 In order to be able to perform horizontal interpolation in this case,
 …
 \end{minted}
 \autoref{fig:obsdataplotmain} shows the main window which is launched when dataplot starts.
+\autoref{fig:OBS_dataplotmain} shows the main window which is launched when dataplot starts.
 This is split into three parts.
 At the top there is a menu bar which contains a variety of drop down menus.
 …
     \includegraphics[width=\textwidth]{Fig_OBS_dataplot_main}
     \caption{
       \protect\label{fig:obsdataplotmain}
+      \protect\label{fig:OBS_dataplotmain}
       Main window of dataplot.
+    }
 …
 If a profile point is clicked with the mouse button a plot of the observation and background values as
 a function of depth (\autoref{fig:obsdataplotprofile}).
+a function of depth (\autoref{fig:OBS_dataplotprofile}).
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
     \includegraphics[width=\textwidth]{Fig_OBS_dataplot_prof}
     \caption{
       \protect\label{fig:obsdataplotprofile}
+      \protect\label{fig:OBS_dataplotprofile}
       Profile plot from dataplot produced by right clicking on a point in the main window.
+    }

NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex

-                      r11537
+                      r11543
 Next, the scheme for interpolation on the fly is described.
 Finally, the different options that further modify the fluxes applied to the ocean are discussed.
 One of these is modification by icebergs (see \autoref{sec:ICB_icebergs}),
+One of these is modification by icebergs (see \autoref{sec:SBC_ICB_icebergs}),
 which act as drifting sources of fresh water.
 Another example of modification is that due to the ice shelf melting/freezing (see \autoref{sec:SBC_isf}),
 …
 The surface ocean stress is the stress exerted by the wind and the sea-ice on the ocean.
 It is applied in \mdl{dynzdf} module as a surface boundary condition of the computation of
 the momentum vertical mixing trend (see \autoref{eq:dynzdf_sbc} in \autoref{sec:DYN_zdf}).
+the momentum vertical mixing trend (see \autoref{eq:DYN_zdf_sbc} in \autoref{sec:DYN_zdf}).
 As such, it has to be provided as a 2D vector interpolated onto the horizontal velocity ocean mesh,
 \ie\ resolved onto the model (\textbf{i},\textbf{j}) direction at $u$- and $v$-points.
 …
 It is applied in \mdl{trasbc} module as a surface boundary condition trend of
 the first level temperature time evolution equation
 (see \autoref{eq:tra_sbc} and \autoref{eq:tra_sbc_lin} in \autoref{subsec:TRA_sbc}).

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

Legend:

NEMO/trunk/doc/latex/NEMO/main/appendices.tex

NEMO/trunk/doc/latex/NEMO/main/chapters.tex

NEMO/trunk/doc/latex/NEMO/main/introduction.tex

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_diff_opers.tex

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_s_coord.tex

NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex