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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9407 for branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex – NEMO

Ignore:
Timestamp:
2018-03-15T17:40:35+01:00 (6 years ago)
Author:
nicolasmartin
Message:

Complete refactoring of cross-referencing

  • Use of \autoref instead of simple \ref for contextual text depending on target type
  • creation of few prefixes for marker to identify the type reference: apdx|chap|eq|fig|sec|subsec|tab
File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex

    r9393 r9407  
    44% Iso-neutral diffusion : 
    55% ================================================================ 
    6 \chapter{Iso-neutral diffusion and eddy advection using triads} 
    7 \label{sec:triad} 
     6\chapter[Iso-Neutral Diffusion and Eddy Advection using Triads] 
     7         {\texorpdfstring{Iso-Neutral Diffusion and\\ Eddy Advection using Triads}{Iso-Neutral Diffusion and Eddy Advection using Triads}} 
     8\label{apdx:triad} 
    89\minitoc 
    910\pagebreak 
    … …  
    1819of iso-neutral diffusion and the eddy-induced advective skew (GM) fluxes.  
    1920If the namelist logical \np{ln\_traldf\_iso} is set true,  
    20 the filtered version of Cox's original scheme (the Standard scheme) is employed (\S\ref{LDF_slp}).  
     21the filtered version of Cox's original scheme (the Standard scheme) is employed (\autoref{sec:LDF_slp}).  
    2122In the present implementation of the Griffies scheme,  
    2223the advective skew fluxes are implemented even if \np{ln\_traldf\_eiv} is false. 
    2324 
    2425Values of iso-neutral diffusivity and GM coefficient are set as 
    25 described in \S\ref{LDF_coef}. Note that when GM fluxes are used,  
     26described in \autoref{sec:LDF_coef}. Note that when GM fluxes are used,  
    2627the eddy-advective (GM) velocities are output for diagnostic purposes using xIOS,  
    2728even though the eddy advection is accomplished by means of the skew fluxes. 
    … …  
    3031The options specific to the Griffies scheme include: 
    3132\begin{description}[font=\normalfont] 
    32 \item[\np{ln\_triad\_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then 
     33\item[\np{ln\_triad\_iso}] See \autoref{sec:taper}. If this is set false (the default), then 
    3334  `iso-neutral' mixing is accomplished within the surface mixed-layer 
    3435  along slopes linearly decreasing with depth from the value immediately below 
    35   the mixed-layer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}).  
    36   This is the same treatment as used in the default implementation \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}.   
     36  the mixed-layer to zero (flat) at the surface (\autoref{sec:lintaper}).  
     37  This is the same treatment as used in the default implementation \autoref{subsec:LDF_slp_iso}; \autoref{fig:eiv_slp}.   
    3738  Where \np{ln\_triad\_iso} is set true, the vertical skew flux is further reduced  
    3839  to ensure no vertical buoyancy flux, giving an almost pure 
    3940  horizontal diffusive tracer flux within the mixed layer. This is similar to 
    40   the tapering suggested by \citet{Gerdes1991}. See \S\ref{sec:triad:Gerdes-taper} 
    41 \item[\np{ln\_botmix\_triad}] See \S\ref{sec:triad:iso_bdry}.  
     41  the tapering suggested by \citet{Gerdes1991}. See \autoref{subsec:Gerdes-taper} 
     42\item[\np{ln\_botmix\_triad}] See \autoref{sec:iso_bdry}.  
    4243  If this is set false (the default) then the lateral diffusive fluxes 
    4344  associated with triads partly masked by topography are neglected.  
    … …  
    5354 
    5455\section{Triad formulation of iso-neutral diffusion} 
    55 \label{sec:triad:iso} 
     56\label{sec:iso} 
    5657We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98},  
    5758but formulated within the \NEMO framework, using scale factors rather than grid-sizes. 
    … …  
    6061The iso-neutral second order tracer diffusive operator for small 
    6162angles between iso-neutral surfaces and geopotentials is given by 
    62 \eqref{Eq_PE_iso_tensor}: 
    63 \begin{subequations} \label{eq:triad:PE_iso_tensor} 
     63\autoref{eq:PE_iso_tensor}: 
     64\begin{subequations} \label{eq:PE_iso_tensor} 
    6465  \begin{equation} 
    6566    D^{lT}=-\Div\vect{f}^{lT}\equiv 
    … …  
    7273  \end{equation} 
    7374  \begin{equation} 
    74     \label{eq:triad:PE_iso_tensor:c} 
     75    \label{eq:PE_iso_tensor:c} 
    7576    \mbox{with}\quad \;\;\Re = 
    7677    \begin{pmatrix} 
    … …  
    9293%  {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ 
    9394% \end{array} }} \right) 
    94  Here \eqref{Eq_PE_iso_slopes}  
     95 Here \autoref{eq:PE_iso_slopes}  
    9596\begin{align*} 
    9697  r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i} 
    … …  
    108109space; we write 
    109110\begin{equation} 
    110   \label{eq:triad:Fijk} 
     111  \label{eq:Fijk} 
    111112  \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right). 
    112113\end{equation} 
    … …  
    117118 
    118119The off-diagonal terms of the small angle diffusion tensor 
    119 \eqref{Eq_PE_iso_tensor}, \eqref{eq:triad:PE_iso_tensor:c} produce skew-fluxes along the 
     120\autoref{eq:PE_iso_tensor}, \autoref{eq:PE_iso_tensor:c} produce skew-fluxes along the 
    120121$i$- and $j$-directions resulting from the vertical tracer gradient: 
    121122\begin{align} 
    122   \label{eq:triad:i13c} 
     123  \label{eq:i13c} 
    123124  f_{13}=&+\Alt r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad f_{23}=+\Alt r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\ 
    124125\intertext{and in the k-direction resulting from the lateral tracer gradients} 
    125   \label{eq:triad:i31c} 
     126  \label{eq:i31c} 
    126127 f_{31}+f_{32}=& \Alt r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\Alt r_2\frac{1}{e_1}\frac{\partial T}{\partial i} 
    127128\end{align} 
    … …  
    130131component of the small angle diffusion tensor is 
    131132\begin{equation} 
    132   \label{eq:triad:i33c} 
     133  \label{eq:i33c} 
    133134  f_{33}=-\Alt(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}. 
    134135\end{equation} 
    … …  
    141142 
    142143There is no natural discretization for the $i$-component of the 
    143 skew-flux, \eqref{eq:triad:i13c}, as 
     144skew-flux, \autoref{eq:i13c}, as 
    144145although it must be evaluated at $u$-points, it involves vertical 
    145146gradients (both for the tracer and the slope $r_1$), defined at 
    146 $w$-points. Similarly, the vertical skew flux, \eqref{eq:triad:i31c}, is evaluated at 
     147$w$-points. Similarly, the vertical skew flux, \autoref{eq:i31c}, is evaluated at 
    147148$w$-points but involves horizontal gradients defined at $u$-points. 
    148149 
    149150\subsection{Standard discretization} 
    150151The straightforward approach to discretize the lateral skew flux 
    151 \eqref{eq:triad:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 
    152 into OPA, \eqref{Eq_tra_ldf_iso}, is to calculate a mean vertical 
     152\autoref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 
     153into OPA, \autoref{eq:tra_ldf_iso}, is to calculate a mean vertical 
    153154gradient at the $u$-point from the average of the four surrounding 
    154155vertical tracer gradients, and multiply this by a mean slope at the 
    … …  
    159160$e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with 
    160161the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer 
    161 gradient, is then \eqref{Eq_tra_ldf_iso} 
     162gradient, is then \autoref{eq:tra_ldf_iso} 
    162163\begin{equation*} 
    163164  \left(F_u^{13} \right)_{i+\hhalf}^k = \Alts_{i+\hhalf}^k 
    … …  
    181182operator to a tracer does not guarantee the decrease of its 
    182183global-average variance. To correct this, we introduced a smoothing of 
    183 the slopes of the iso-neutral surfaces (see \S\ref{LDF}). This 
     184the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). This 
    184185technique works for $T$ and $S$ in so far as they are active tracers 
    185186($i.e.$ they enter the computation of density), but it does not work 
    … …  
    194195\begin{figure}[tb] \begin{center} 
    195196    \includegraphics[width=1.05\textwidth]{Fig_GRIFF_triad_fluxes} 
    196     \caption{ \protect\label{fig:triad:ISO_triad} 
     197    \caption{ \protect\label{fig:ISO_triad} 
    197198      (a) Arrangement of triads $S_i$ and tracer gradients to 
    198199           give lateral tracer flux from box $i,k$ to $i+1,k$ 
    … …  
    205206slope calculated from the lateral density gradient across the $u$-point 
    206207divided by the vertical density gradient at the same $w$-point as the 
    207 tracer gradient. See Fig.~\ref{fig:triad:ISO_triad}a, where the thick lines 
     208tracer gradient. See \autoref{fig:ISO_triad}a, where the thick lines 
    208209denote the tracer gradients, and the thin lines the corresponding 
    209210triads, with slopes $s_1, \dotsc s_4$. The total area-integrated 
    210211skew-flux from tracer cell $i,k$ to $i+1,k$ 
    211212\begin{multline} 
    212   \label{eq:triad:i13} 
     213  \label{eq:i13} 
    213214  \left( F_u^{13}  \right)_{i+\frac{1}{2}}^k = \Alts_{i+1}^k a_1 s_1 
    214215  \delta _{k+\frac{1}{2}} \left[ T^{i+1} 
    … …  
    225226stencil, and disallows the two-point computational modes. 
    226227 
    227  The vertical skew flux \eqref{eq:triad:i31c} from tracer cell $i,k$ to $i,k+1$ at the 
    228 $w$-point $i,k+\hhalf$ is constructed similarly (Fig.~\ref{fig:triad:ISO_triad}b) 
     228 The vertical skew flux \autoref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at the 
     229$w$-point $i,k+\hhalf$ is constructed similarly (\autoref{fig:ISO_triad}b) 
    229230by multiplying lateral tracer gradients from each of the four 
    230231surrounding $u$-points by the appropriate triad slope: 
    231232\begin{multline} 
    232   \label{eq:triad:i31} 
     233  \label{eq:i31} 
    233234  \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} =  \Alts_i^{k+1} a_{1}' 
    234235  s_{1}' \delta _{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1} 
    … …  
    241242(appearing in both the vertical and lateral gradient), and the $u$- and 
    242243$w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the 
    243 triad as follows (see also Fig.~\ref{fig:triad:ISO_triad}): 
    244 \begin{equation} 
    245   \label{eq:triad:R} 
     244triad as follows (see also \autoref{fig:ISO_triad}): 
     245\begin{equation} 
     246  \label{eq:R} 
    246247  _i^k \mathbb{R}_{i_p}^{k_p} 
    247248  =-\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} 
    … …  
    258259\begin{figure}[tb] \begin{center} 
    259260    \includegraphics[width=0.80\textwidth]{Fig_GRIFF_qcells} 
    260     \caption{   \protect\label{fig:triad:qcells} 
     261    \caption{   \protect\label{fig:qcells} 
    261262    Triad notation for quarter cells. $T$-cells are inside 
    262263      boxes, while the  $i+\half,k$ $u$-cell is shaded in green and the 
    … …  
    265266% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    266267 
    267 Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter 
     268Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:qcells}) with the quarter 
    268269cell 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.  
    269 Expressing the slopes $s_i$ and $s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation,  
     270Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:i13} and \autoref{eq:i31} in this notation,  
    270271we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$.  
    271272Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$)  
    … …  
    276277of a unique triad, and we notate these areas, similarly to the triad slopes,  
    277278as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$,  
    278 where $e.g.$ in \eqref{eq:triad:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,  
    279 and in \eqref{eq:triad:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 
     279where $e.g.$ in \autoref{eq:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,  
     280and in \autoref{eq:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 
    280281 
    281282\subsection{Full triad fluxes} 
    … …  
    287288form 
    288289\begin{equation} 
    289   \label{eq:triad:i11} 
     290  \label{eq:i11} 
    290291  \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} = 
    291292  - \left( \Alts_i^{k+1} a_{1} + \Alts_i^{k+1} a_{2} + \Alts_i^k 
    … …  
    293294  \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}}, 
    294295\end{equation} 
    295 where the areas $a_i$ are as in \eqref{eq:triad:i13}. In this case, 
    296 separating the total lateral flux, the sum of \eqref{eq:triad:i13} and 
    297 \eqref{eq:triad:i11}, into triad components, a lateral tracer 
     296where the areas $a_i$ are as in \autoref{eq:i13}. In this case, 
     297separating the total lateral flux, the sum of \autoref{eq:i13} and 
     298\autoref{eq:i11}, into triad components, a lateral tracer 
    298299flux 
    299300\begin{equation} 
    300   \label{eq:triad:latflux-triad} 
     301  \label{eq:latflux-triad} 
    301302  _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - \Alts_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p} 
    302303  \left( 
    … …  
    312313density flux associated with each triad separately disappears. 
    313314\begin{equation} 
    314   \label{eq:triad:latflux-rho} 
     315  \label{eq:latflux-rho} 
    315316  {\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 
    316317\end{equation} 
    … …  
    319320to $i+1,k$ must also vanish since it is a sum of four such triad fluxes. 
    320321 
    321 The squared slope $r_1^2$ in the expression \eqref{eq:triad:i33c} for the 
     322The squared slope $r_1^2$ in the expression \autoref{eq:i33c} for the 
    322323$_{33}$ component is also expressed in terms of area-weighted 
    323324squared triad slopes, so the area-integrated vertical flux from tracer 
    324325cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is 
    325326\begin{equation} 
    326   \label{eq:triad:i33} 
     327  \label{eq:i33} 
    327328  \left( F_w^{33} \right) _i^{k+\frac{1}{2}} = 
    328329    - \left( \Alts_i^{k+1} a_{1}' s_{1}'^2 
    … …  
    332333\end{equation} 
    333334where the areas $a'$ and slopes $s'$ are the same as in 
    334 \eqref{eq:triad:i31}. 
    335 Then, separating the total vertical flux, the sum of \eqref{eq:triad:i31} and 
    336 \eqref{eq:triad:i33}, into triad components,  a vertical flux 
     335\autoref{eq:i31}. 
     336Then, separating the total vertical flux, the sum of \autoref{eq:i31} and 
     337\autoref{eq:i33}, into triad components,  a vertical flux 
    337338\begin{align} 
    338   \label{eq:triad:vertflux-triad} 
     339  \label{eq:vertflux-triad} 
    339340  _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) 
    340341  &= \Alts_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p} 
    … …  
    345346  \right) \\ 
    346347  &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right) 
    347    {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:triad:vertflux-triad2} 
     348   {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2} 
    348349\end{align} 
    349350may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$ 
    … …  
    355356fluxes. 
    356357 
    357 We can explicitly identify (Fig.~\ref{fig:triad:qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of 
     358We 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 
    358359the $u$-fluxes and $w$-fluxes in 
    359 \eqref{eq:triad:i31}, \eqref{eq:triad:i13}, \eqref{eq:triad:i11} \eqref{eq:triad:i33} and 
    360 Fig.~\ref{fig:triad:ISO_triad} to  write out the iso-neutral fluxes at $u$- and 
     360\autoref{eq:i31}, \autoref{eq:i13}, \autoref{eq:i11} \autoref{eq:i33} and 
     361\autoref{fig:ISO_triad} to  write out the iso-neutral fluxes at $u$- and 
    361362$w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces: 
    362 %(Fig.~\ref{Fig_ISO_triad}): 
    363 \begin{flalign} \label{Eq_iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv 
     363%(\autoref{fig:ISO_triad}): 
     364\begin{flalign} \label{eq:iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv 
    364365  \sum_{\substack{i_p,\,k_p}} 
    365366  \begin{pmatrix} 
    … …  
    371372 
    372373\subsection{Ensuring the scheme does not increase tracer variance} 
    373 \label{sec:triad:variance} 
     374\label{subsec:variance} 
    374375 
    375376We now require that this operator should not increase the 
    … …  
    397398  &= -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 
    398399  {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\ 
    399   &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:triad:dvar_iso_i} 
     400  &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:dvar_iso_i} 
    400401 \end{aligned} 
    401402\end{multline} 
    … …  
    404405$i,k+k_p+\half$ (below) of 
    405406\begin{equation} 
    406 \label{eq:triad:dvar_iso_k} 
     407\label{eq:dvar_iso_k} 
    407408  _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 
    408409\end{equation} 
    409410The total variance tendency driven by the triad is the sum of these 
    410411two. Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and 
    411 $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \eqref{eq:triad:latflux-triad} and 
    412 \eqref{eq:triad:vertflux-triad}, it is 
     412$_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \autoref{eq:latflux-triad} and 
     413\autoref{eq:vertflux-triad}, it is 
    413414\begin{multline*} 
    414415-\Alts_i^k\left \{ 
    … …  
    430431to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by 
    431432\begin{equation} 
    432   \label{eq:triad:V-A} 
     433  \label{eq:V-A} 
    433434  _i^k\mathbb{V}_{i_p}^{k_p} 
    434435  ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k} 
    … …  
    437438the variance tendency reduces to the perfect square 
    438439\begin{equation} 
    439   \label{eq:triad:perfect-square} 
     440  \label{eq:perfect-square} 
    440441  -\Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 
    441442  \left( 
    … …  
    445446  \right)^2\leq 0. 
    446447\end{equation} 
    447 Thus, the constraint \eqref{eq:triad:V-A} ensures that the fluxes (\ref{eq:triad:latflux-triad}, \ref{eq:triad:vertflux-triad}) associated 
     448Thus, the constraint \autoref{eq:V-A} ensures that the fluxes (\autoref{eq:latflux-triad}, \autoref{eq:vertflux-triad}) associated 
    448449with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase 
    449450the net variance. Since the total fluxes are sums of such fluxes from 
    … …  
    452453increase. 
    453454 
    454 The expression \eqref{eq:triad:V-A} can be interpreted as a discretization 
     455The expression \autoref{eq:V-A} can be interpreted as a discretization 
    455456of the global integral 
    456457\begin{equation} 
    457   \label{eq:triad:cts-var} 
     458  \label{eq:cts-var} 
    458459  \frac{\partial}{\partial t}\int\!\half T^2\, dV = 
    459460  \int\!\mathbf{F}\cdot\nabla T\, dV, 
    … …  
    480481cells, defined in terms of the distances between $T$, $u$,$f$ and 
    481482$w$-points. This is the natural discretization of 
    482 \eqref{eq:triad:cts-var}. The \NEMO model, however, operates with scale 
     483\autoref{eq:cts-var}. The \NEMO model, however, operates with scale 
    483484factors instead of grid sizes, and scale factors for the quarter 
    484485cells are not defined. Instead, therefore we simply choose 
    485486\begin{equation} 
    486   \label{eq:triad:V-NEMO} 
     487  \label{eq:V-NEMO} 
    487488  _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k, 
    488489\end{equation} 
    … …  
    492493$i+1,k$ reduces to the classical form 
    493494\begin{equation} 
    494   \label{eq:triad:lat-normal} 
     495  \label{eq:lat-normal} 
    495496-\overline\Alts_{\,i+1/2}^k\; 
    496497\frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    … …  
    500501In fact if the diffusive coefficient is defined at $u$-points, so that 
    501502we employ $\Alts_{i+i_p}^k$ instead of  $\Alts_i^k$ in the definitions of the 
    502 triad fluxes \eqref{eq:triad:latflux-triad} and \eqref{eq:triad:vertflux-triad}, 
     503triad fluxes \autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad}, 
    503504we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above. 
    504505 
    … …  
    506507The iso-neutral fluxes at $u$- and 
    507508$w$-points are the sums of the triad fluxes that cross the $u$- and 
    508 $w$-faces \eqref{Eq_iso_flux}: 
    509 \begin{subequations}\label{eq:triad:alltriadflux} 
    510   \begin{flalign}\label{eq:triad:vect_isoflux} 
     509$w$-faces \autoref{eq:iso_flux}: 
     510\begin{subequations}\label{eq:alltriadflux} 
     511  \begin{flalign}\label{eq:vect_isoflux} 
    511512    \vect{F}_{\mathrm{iso}}(T) &\equiv 
    512513    \sum_{\substack{i_p,\,k_p}} 
    … …  
    517518    \end{pmatrix}, 
    518519  \end{flalign} 
    519   where \eqref{eq:triad:latflux-triad}: 
     520  where \autoref{eq:latflux-triad}: 
    520521  \begin{align} 
    521     \label{eq:triad:triadfluxu} 
     522    \label{eq:triadfluxu} 
    522523    _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - \Alts_i^k{ 
    523524      \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    … …  
    534535      -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ 
    535536      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
    536     \right),\label{eq:triad:triadfluxw} 
     537    \right),\label{eq:triadfluxw} 
    537538  \end{align} 
    538   with \eqref{eq:triad:V-NEMO} 
     539  with \autoref{eq:V-NEMO} 
    539540  \begin{equation} 
    540     \label{eq:triad:V-NEMO2} 
     541    \label{eq:V-NEMO2} 
    541542    _i^k{\mathbb{V}}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k. 
    542543  \end{equation} 
    543544\end{subequations} 
    544545 
    545  The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 
     546 The divergence of the expression \autoref{eq:iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 
    546547each tracer point: 
    547 \begin{equation} \label{eq:triad:iso_operator} D_l^T = \frac{1}{b_T} 
     548\begin{equation} \label{eq:iso_operator} D_l^T = \frac{1}{b_T} 
    548549  \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k 
    549550        {\mathbb{F}_u }_{i_p}^{k_p}} \right] + \delta_{k} \left[ 
    … …  
    554555\begin{description} 
    555556\item[$\bullet$ horizontal diffusion] The discretization of the 
    556   diffusion operator recovers \eqref{eq:triad:lat-normal} the traditional five-point Laplacian in 
     557  diffusion operator recovers \autoref{eq:lat-normal} the traditional five-point Laplacian in 
    557558  the limit of flat iso-neutral direction : 
    558   \begin{equation} \label{eq:triad:iso_property0} D_l^T = \frac{1}{b_T} \ 
     559  \begin{equation} \label{eq:iso_property0} D_l^T = \frac{1}{b_T} \ 
    559560    \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; 
    560561      \overline\Alts^{\,i} \; \delta_{i+1/2}[T] \right] \qquad 
    … …  
    564565\item[$\bullet$ implicit treatment in the vertical] Only tracer values 
    565566  associated with a single water column appear in the expression 
    566   \eqref{eq:triad:i33} for the $_{33}$ fluxes, vertical fluxes driven by 
     567  \autoref{eq:i33} for the $_{33}$ fluxes, vertical fluxes driven by 
    567568  vertical gradients. This is of paramount importance since it means 
    568569  that a time-implicit algorithm can be used to solve the vertical 
    … …  
    582583\item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of 
    583584  locally referenced potential density is zero. See 
    584   \eqref{eq:triad:latflux-rho} and \eqref{eq:triad:vertflux-triad2}. 
     585  \autoref{eq:latflux-rho} and \autoref{eq:vertflux-triad2}. 
    585586 
    586587\item[$\bullet$ conservation of tracer] The iso-neutral diffusion 
    587588  conserves tracer content, $i.e.$ 
    588   \begin{equation} \label{eq:triad:iso_property1} \sum_{i,j,k} \left\{ D_l^T \ 
     589  \begin{equation} \label{eq:iso_property1} \sum_{i,j,k} \left\{ D_l^T \ 
    589590      b_T \right\} = 0 
    590591  \end{equation} 
    … …  
    594595\item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion 
    595596  does not increase the tracer variance, $i.e.$ 
    596   \begin{equation} \label{eq:triad:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T 
     597  \begin{equation} \label{eq:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T 
    597598      \ b_T \right\} \leq 0 
    598599  \end{equation} 
    599600  The property is demonstrated in 
    600   \S\ref{sec:triad:variance} above. It is a key property for a diffusion 
     601  \autoref{subsec:variance} above. It is a key property for a diffusion 
    601602  term. It means that it is also a dissipation term, $i.e.$ it 
    602603  dissipates the square of the quantity on which it is applied.  It 
    … …  
    607608\item[$\bullet$ self-adjoint operator] The iso-neutral diffusion 
    608609  operator is self-adjoint, $i.e.$ 
    609   \begin{equation} \label{eq:triad:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T 
     610  \begin{equation} \label{eq:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T 
    610611      \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 
    611612  \end{equation} 
    … …  
    614615  routine. This property can be demonstrated similarly to the proof of 
    615616  the `no increase of tracer variance' property. The contribution by a 
    616   single triad towards the left hand side of \eqref{eq:triad:iso_property3}, can 
    617   be found by replacing $\delta[T]$ by $\delta[S]$ in \eqref{eq:triad:dvar_iso_i} 
    618   and \eqref{eq:triad:dvar_iso_k}. This results in a term similar to 
    619   \eqref{eq:triad:perfect-square}, 
    620 \begin{equation} 
    621   \label{eq:triad:TScovar} 
     617  single triad towards the left hand side of \autoref{eq:iso_property3}, can 
     618  be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:dvar_iso_i} 
     619  and \autoref{eq:dvar_iso_k}. This results in a term similar to 
     620  \autoref{eq:perfect-square}, 
     621\begin{equation} 
     622  \label{eq:TScovar} 
    622623  - \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 
    623624  \left( 
    … …  
    634635This is symmetrical in $T $ and $S$, so exactly the same term arises 
    635636from the discretization of this triad's contribution towards the 
    636 RHS of \eqref{eq:triad:iso_property3}. 
     637RHS of \autoref{eq:iso_property3}. 
    637638\end{description} 
    638639 
    639 \subsection{Treatment of the triads at the boundaries}\label{sec:triad:iso_bdry} 
     640\subsection{Treatment of the triads at the boundaries}\label{sec:iso_bdry} 
    640641The triad slope can only be defined where both the grid boxes centred at 
    641642the end of the arms exist. Triads that would poke up 
    642643through the upper ocean surface into the atmosphere, or down into the 
    643 ocean floor, must be masked out. See Fig.~\ref{fig:triad:bdry_triads}. Surface layer triads 
     644ocean floor, must be masked out. See \autoref{fig:bdry_triads}. Surface layer triads 
    644645$\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and 
    645646$\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be 
    646 specified above the ocean surface are masked (Fig.~\ref{fig:triad:bdry_triads}a): this ensures that lateral 
     647specified above the ocean surface are masked (\autoref{fig:bdry_triads}a): this ensures that lateral 
    647648tracer gradients produce no flux through the ocean surface. However, 
    648649to prevent surface noise, it is customary to retain the $_{11}$ contributions towards 
    … …  
    650651$\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer 
    651652fluxes. Similar comments apply to triads that would intersect the 
    652 ocean floor (Fig.~\ref{fig:triad:bdry_triads}b). Note that both near bottom 
     653ocean floor (\autoref{fig:bdry_triads}b). Note that both near bottom 
    653654triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and 
    654655$\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ 
    … …  
    665666\begin{figure}[h] \begin{center} 
    666667    \includegraphics[width=0.60\textwidth]{Fig_GRIFF_bdry_triads} 
    667     \caption{  \protect\label{fig:triad:bdry_triads} 
     668    \caption{  \protect\label{fig:bdry_triads} 
    668669      (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer 
    669670      points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad 
    … …  
    678679      or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point 
    679680      is masked. The associated lateral fluxes (grey-black dashed 
    680       line) are masked if \np{botmix\_triad}\forcode{ = .false.}, but left 
    681       unmasked, giving bottom mixing, if \np{botmix\_triad}\forcode{ = .true.}} 
     681      line) are masked if \protect\np{botmix\_triad}\forcode{ = .false.}, but left 
     682      unmasked, giving bottom mixing, if \protect\np{botmix\_triad}\forcode{ = .true.}} 
    682683 \end{center} \end{figure} 
    683684% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    684685 
    685 \subsection{ Limiting of the slopes within the interior}\label{sec:triad:limit} 
    686 As discussed in \S\ref{LDF_slp_iso}, iso-neutral slopes relative to 
     686\subsection{ Limiting of the slopes within the interior}\label{sec:limit} 
     687As discussed in \autoref{subsec:LDF_slp_iso}, iso-neutral slopes relative to 
    687688geopotentials must be bounded everywhere, both for consistency with the small-slope 
    688689approximation and for numerical stability \citep{Cox1987, 
    … …  
    692693It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to 
    693694geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to 
    694 geopotentials) \eqref{Eq_PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate 
     695geopotentials) \autoref{eq:PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate 
    695696surfaces, so we require 
    696697\begin{equation*} 
    … …  
    700701Each individual triad slope 
    701702 \begin{equation} 
    702    \label{eq:triad:Rtilde} 
     703   \label{eq:Rtilde} 
    703704_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}} 
    704705 \end{equation} 
    … …  
    709710is always downwards, and so acts to reduce gravitational potential energy. 
    710711 
    711 \subsection{Tapering within the surface mixed layer}\label{sec:triad:taper} 
     712\subsection{Tapering within the surface mixed layer}\label{sec:taper} 
    712713Additional tapering of the iso-neutral fluxes is necessary within the 
    713714surface mixed layer. When the Griffies triads are used, we offer two 
    714715options for this. 
    715716 
    716 \subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper} 
     717\subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:lintaper} 
    717718This is the option activated by the default choice 
    718719\np{ln\_triad\_iso}\forcode{ = .false.}. Slopes $\tilde{r}_i$ relative to 
    719720geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the 
    720 surface, as described in option (c) of Fig.~\ref{Fig_eiv_slp}, to values 
     721surface, as described in option (c) of \autoref{fig:eiv_slp}, to values 
    721722\begin{subequations} 
    722723  \begin{equation} 
    723    \label{eq:triad:rmtilde} 
     724   \label{eq:rmtilde} 
    724725     \rMLt = 
    725726  -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h, 
    … …  
    728729adjusted to 
    729730  \begin{equation} 
    730    \label{eq:triad:rm} 
     731   \label{eq:rm} 
    731732 \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h. 
    732733  \end{equation} 
    733734\end{subequations} 
    734735Thus the diffusion operator within the mixed layer is given by: 
    735 \begin{equation} \label{eq:triad:iso_tensor_ML} 
     736\begin{equation} \label{eq:iso_tensor_ML} 
    736737D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    737738\mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} 
    … …  
    745746mixed-layer and in isopycnal layers immediately below, in the 
    746747thermocline. It is consistent with the way the $\tilde{r}_i$ are 
    747 tapered within the mixed layer (see \S\ref{sec:triad:taperskew} below) 
     748tapered within the mixed layer (see \autoref{sec:taperskew} below) 
    748749so as to ensure a uniform GM eddy-induced velocity throughout the 
    749750mixed layer. However, it gives a downwards density flux and so acts so 
    750751as to reduce potential energy in the same way as does the slope 
    751 limiting discussed above in \S\ref{sec:triad:limit}. 
     752limiting discussed above in \autoref{sec:limit}. 
    752753  
    753 As in \S\ref{sec:triad:limit} above, the tapering 
    754 \eqref{eq:triad:rmtilde} is applied separately to each triad 
     754As in \autoref{sec:limit} above, the tapering 
     755\autoref{eq:rmtilde} is applied separately to each triad 
    755756$_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the 
    756757$_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume 
    757758$z$-coordinates in the following; the conversion from 
    758759$\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described 
    759 above by \eqref{eq:triad:Rtilde}. 
     760above by \autoref{eq:Rtilde}. 
    760761\begin{enumerate} 
    761762\item Mixed-layer depth is defined so as to avoid including regions of weak 
    762763vertical stratification in the slope definition. 
    763764 At each $i,j$ (simplified to $i$ in 
    764 Fig.~\ref{fig:triad:MLB_triad}), we define the mixed-layer by setting 
     765\autoref{fig:MLB_triad}), we define the mixed-layer by setting 
    765766the vertical index of the tracer point immediately below the mixed 
    766767layer, $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) 
    … …  
    768769${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, where $i,k_{10}$ is 
    769770the tracer gridbox within which the depth reaches 10~m. See the left 
    770 side of Fig.~\ref{fig:triad:MLB_triad}. We use the $k_{10}$-gridbox 
     771side of \autoref{fig:MLB_triad}. We use the $k_{10}$-gridbox 
    771772instead of the surface gridbox to avoid problems e.g.\ with thin 
    772773daytime mixed-layers. Currently we use the same 
    … …  
    784785${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are 
    785786representative of the thermocline. The four basal triads defined in the bottom part 
    786 of Fig.~\ref{fig:triad:MLB_triad} are then 
     787of \autoref{fig:MLB_triad} are then 
    787788\begin{align} 
    788789  {\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p} &= 
    789  {\:}^{k_{\mathrm{ML}}-k_p-1/2}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}, \label{eq:triad:Rbase} 
     790 {\:}^{k_{\mathrm{ML}}-k_p-1/2}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}, \label{eq:Rbase} 
    790791\\ 
    791792\intertext{with e.g.\ the green triad} 
    … …  
    797798so it is this depth 
    798799\begin{equation} 
    799   \label{eq:triad:zbase} 
     800  \label{eq:zbase} 
    800801  {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2} 
    801802\end{equation} 
    802803(one gridbox deeper than the 
    803804diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper 
    804 the slopes in \eqref{eq:triad:rmtilde}. 
     805the slopes in \autoref{eq:rmtilde}. 
    805806\item Finally, we calculate the adjusted triads 
    806807${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within the mixed 
    … …  
    815816\intertext{and more generally} 
    816817 {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p} &= 
    817 \frac{{z_w}_{k+k_p}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}.\label{eq:triad:RML} 
     818\frac{{z_w}_{k+k_p}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}.\label{eq:RML} 
    818819\end{align} 
    819820\end{enumerate} 
    … …  
    822823\begin{figure}[h] 
    823824%  \fcapside { 
    824     \caption{\protect\label{fig:triad:MLB_triad} Definition of 
     825    \caption{\protect\label{fig:MLB_triad} Definition of 
    825826      mixed-layer depth and calculation of linearly tapered 
    826827      triads. The figure shows a water column at a given $i,j$ 
    … …  
    846847 
    847848\subsubsection{Additional truncation of skew iso-neutral flux components} 
    848 \label{sec:triad:Gerdes-taper} 
     849\label{subsec:Gerdes-taper} 
    849850The alternative option is activated by setting \np{ln\_triad\_iso} = 
    850851  true. This retains the same tapered slope $\rML$  described above for the 
    … …  
    853854replaces the $\rML$ in the skew term by 
    854855\begin{equation} 
    855   \label{eq:triad:rm*} 
     856  \label{eq:rm*} 
    856857  \rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i, 
    857858\end{equation} 
    858859giving a ML diffusive operator 
    859 \begin{equation} \label{eq:triad:iso_tensor_ML2} 
     860\begin{equation} \label{eq:iso_tensor_ML2} 
    860861D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    861862\mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} 
    … …  
    887888% Skew flux formulation for Eddy Induced Velocity : 
    888889% ================================================================ 
    889 \section{Eddy induced advection formulated as a skew flux}\label{sec:triad:skew-flux} 
    890  
    891 \subsection{Continuous skew flux formulation}\label{sec:triad:continuous-skew-flux} 
     890\section{Eddy induced advection formulated as a skew flux}\label{sec:skew-flux} 
     891 
     892\subsection{Continuous skew flux formulation}\label{sec:continuous-skew-flux} 
    892893 
    893894 When Gent and McWilliams's [1990] diffusion is used, 
    … …  
    895896eddy induced velocity, the formulation of which depends on the slopes of iso- 
    896897neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used 
    897 here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo} 
    898 is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo} 
    899 + \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates. 
     898here are referenced to the geopotential surfaces, $i.e.$ \autoref{eq:ldfslp_geo} 
     899is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo} 
     900+ \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. 
    900901 
    901902The eddy induced velocity is given by: 
    902 \begin{subequations} \label{eq:triad:eiv} 
    903 \begin{equation}\label{eq:triad:eiv_v} 
     903\begin{subequations} \label{eq:eiv} 
     904\begin{equation}\label{eq:eiv_v} 
    904905\begin{split} 
    905906 u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\ 
    … …  
    910911\end{equation} 
    911912where the streamfunctions $\psi_i$ are given by 
    912 \begin{equation} \label{eq:triad:eiv_psi} 
     913\begin{equation} \label{eq:eiv_psi} 
    913914\begin{split} 
    914915\psi_1 & = A_{e} \; \tilde{r}_1,   \\ 
    … …  
    924925default implementation, where \np{ln\_traldf\_triad} is set 
    925926false. This allows us to take advantage of all the advection schemes 
    926 offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$ 
     927offered for the tracers (see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ 
    927928order advection scheme. This is particularly useful for passive 
    928929tracers where \emph{positivity} of the advection scheme is of 
    … …  
    962963and since the eddy induced velocity field is non-divergent, we end up with the skew 
    963964form of the eddy induced advective fluxes per unit area in $ijk$ space: 
    964 \begin{equation} \label{eq:triad:eiv_skew_ijk} 
     965\begin{equation} \label{eq:eiv_skew_ijk} 
    965966\textbf{F}_\mathrm{eiv}^T = \begin{pmatrix} 
    966967           {+ e_{2} \, \psi_1  \; \partial_k T}   \\ 
    … …  
    969970\end{equation} 
    970971The total fluxes per unit physical area are then 
    971 \begin{equation}\label{eq:triad:eiv_skew_physical} 
     972\begin{equation}\label{eq:eiv_skew_physical} 
    972973\begin{split} 
    973974 f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T   \\ 
    … …  
    977978\end{split} 
    978979\end{equation} 
    979 Note that Eq.~ \eqref{eq:triad:eiv_skew_physical} takes the same form whatever the 
     980Note that  \autoref{eq:eiv_skew_physical} takes the same form whatever the 
    980981vertical coordinate, though of course the slopes 
    981 $\tilde{r}_i$ which define the $\psi_i$ in \eqref{eq:triad:eiv_psi} are relative to geopotentials. 
     982$\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:eiv_psi} are relative to geopotentials. 
    982983The tendency associated with eddy induced velocity is then simply the convergence 
    983 of the fluxes (\ref{eq:triad:eiv_skew_ijk}, \ref{eq:triad:eiv_skew_physical}), so 
    984 \begin{equation} \label{eq:triad:skew_eiv_conv} 
     984of the fluxes (\autoref{eq:eiv_skew_ijk}, \autoref{eq:eiv_skew_physical}), so 
     985\begin{equation} \label{eq:skew_eiv_conv} 
    985986\frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[ 
    986987  \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right) 
    … …  
    995996 
    996997\subsection{Discrete skew flux formulation} 
    997 The skew fluxes in (\ref{eq:triad:eiv_skew_physical}, \ref{eq:triad:eiv_skew_ijk}), like the off-diagonal terms 
    998 (\ref{eq:triad:i13c}, \ref{eq:triad:i31c}) of the small angle diffusion tensor, are best 
    999 expressed in terms of the triad slopes, as in Fig.~\ref{fig:triad:ISO_triad} 
    1000 and Eqs~(\ref{eq:triad:i13}, \ref{eq:triad:i31}); but now in terms of the triad slopes 
     998The skew fluxes in (\autoref{eq:eiv_skew_physical}, \autoref{eq:eiv_skew_ijk}), like the off-diagonal terms 
     999(\autoref{eq:i13c}, \autoref{eq:i31c}) of the small angle diffusion tensor, are best 
     1000expressed in terms of the triad slopes, as in \autoref{fig:ISO_triad} 
     1001and (\autoref{eq:i13}, \autoref{eq:i31}); but now in terms of the triad slopes 
    10011002$\tilde{\mathbb{R}}$ relative to geopotentials instead of the 
    10021003$\mathbb{R}$ relative to coordinate surfaces. The discrete form of 
    1003 \eqref{eq:triad:eiv_skew_ijk} using the slopes \eqref{eq:triad:R} and 
     1004\autoref{eq:eiv_skew_ijk} using the slopes \autoref{eq:R} and 
    10041005defining $A_e$ at $T$-points is then given by: 
    10051006 
    10061007 
    1007 \begin{subequations}\label{eq:triad:allskewflux} 
    1008   \begin{flalign}\label{eq:triad:vect_skew_flux} 
     1008\begin{subequations}\label{eq:allskewflux} 
     1009  \begin{flalign}\label{eq:vect_skew_flux} 
    10091010    \vect{F}_{\mathrm{eiv}}(T) &\equiv 
    10101011    \sum_{\substack{i_p,\,k_p}} 
    … …  
    10161017  \end{flalign} 
    10171018  where the skew flux in the $i$-direction associated with a given 
    1018   triad is (\ref{eq:triad:latflux-triad}, \ref{eq:triad:triadfluxu}): 
     1019  triad is (\autoref{eq:latflux-triad}, \autoref{eq:triadfluxu}): 
    10191020  \begin{align} 
    1020     \label{eq:triad:skewfluxu} 
     1021    \label{eq:skewfluxu} 
    10211022    _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \quarter {A_e}_i^k{ 
    10221023      \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    … …  
    10241025      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }, 
    10251026   \\ 
    1026     \intertext{and \eqref{eq:triad:triadfluxw} in the $k$-direction, changing the sign 
    1027       to be consistent with \eqref{eq:triad:eiv_skew_ijk}:} 
     1027    \intertext{and \autoref{eq:triadfluxw} in the $k$-direction, changing the sign 
     1028      to be consistent with \autoref{eq:eiv_skew_ijk}:} 
    10281029    _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T) 
    10291030    &= -\quarter {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}} 
    1030      {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:triad:skewfluxw} 
     1031     {_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} 
    10311032  \end{align} 
    10321033\end{subequations} 
    … …  
    10401041include a diffusive component but is a `pure' advection term. This can 
    10411042be seen 
    1042 %either from Appendix \ref{Apdx_eiv_skew} or 
     1043%either from Appendix \autoref{apdx:eiv_skew} or 
    10431044by considering the 
    10441045fluxes associated with a given triad slope 
    10451046$_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. For, following 
    1046 \S\ref{sec:triad:variance} and \eqref{eq:triad:dvar_iso_i}, the 
     1047\autoref{subsec:variance} and \autoref{eq:dvar_iso_i}, the 
    10471048associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ 
    10481049drives a net rate of change of variance, summed over the two 
    10491050$T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of 
    10501051\begin{equation} 
    1051 \label{eq:triad:dvar_eiv_i} 
     1052\label{eq:dvar_eiv_i} 
    10521053  _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], 
    10531054\end{equation} 
    … …  
    10551056$T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of 
    10561057\begin{equation} 
    1057 \label{eq:triad:dvar_eiv_k} 
     1058\label{eq:dvar_eiv_k} 
    10581059  _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 
    10591060\end{equation} 
    1060 Inspection of the definitions (\ref{eq:triad:skewfluxu}, \ref{eq:triad:skewfluxw}) 
    1061 shows that these two variance changes (\ref{eq:triad:dvar_eiv_i}, \ref{eq:triad:dvar_eiv_k}) 
     1061Inspection of the definitions (\autoref{eq:skewfluxu}, \autoref{eq:skewfluxw}) 
     1062shows that these two variance changes (\autoref{eq:dvar_eiv_i}, \autoref{eq:dvar_eiv_k}) 
    10621063sum to zero. Hence the two fluxes associated with each triad make no 
    10631064net contribution to the variance budget. 
    … …  
    10721073For the change in gravitational PE driven by the $k$-flux is 
    10731074\begin{align} 
    1074   \label{eq:triad:vert_densityPE} 
     1075  \label{eq:vert_densityPE} 
    10751076  g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) 
    10761077  &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k 
    … …  
    10781079    {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\ 
    10791080\intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from 
    1080   \eqref{eq:triad:skewfluxw}, gives} 
     1081  \autoref{eq:skewfluxw}, gives} 
    10811082% and separating out 
    10821083% $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$, 
    … …  
    10901091\end{align} 
    10911092using the definition of the triad slope $\rtriad{R}$, 
    1092 \eqref{eq:triad:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+ 
     1093\autoref{eq:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+ 
    10931094\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of  $-\alpha_i^k \delta_{k+ 
    10941095  k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$. 
    … …  
    10961097Where the coordinates slope, the $i$-flux gives a PE change 
    10971098\begin{multline} 
    1098   \label{eq:triad:lat_densityPE} 
     1099  \label{eq:lat_densityPE} 
    10991100 g \delta_{i+i_p}[z_T^k] 
    11001101\left[ 
    … …  
    11061107\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}}, 
    11071108\end{multline} 
    1108 (using \eqref{eq:triad:skewfluxu}) and so the total PE change 
    1109 \eqref{eq:triad:vert_densityPE} + \eqref{eq:triad:lat_densityPE} associated with the triad fluxes is 
     1109(using \autoref{eq:skewfluxu}) and so the total PE change 
     1110\autoref{eq:vert_densityPE} + \autoref{eq:lat_densityPE} associated with the triad fluxes is 
    11101111\begin{multline} 
    1111   \label{eq:triad:tot_densityPE} 
     1112  \label{eq:tot_densityPE} 
    11121113  g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) + 
    11131114g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\ 
    … …  
    11191120\beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative. 
    11201121 
    1121 \subsection{Treatment of the triads at the boundaries}\label{sec:triad:skew_bdry} 
     1122\subsection{Treatment of the triads at the boundaries}\label{sec:skew_bdry} 
    11221123Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes 
    11231124are masked at the boundaries in exactly the same way as are the triad 
    11241125slopes \rtriad{R} used for the iso-neutral diffusive fluxes, as 
    1125 described in \S\ref{sec:triad:iso_bdry} and 
    1126 Fig.~\ref{fig:triad:bdry_triads}. Thus surface layer triads 
     1126described in \autoref{sec:iso_bdry} and 
     1127\autoref{fig:bdry_triads}. Thus surface layer triads 
    11271128$\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are 
    11281129masked, and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ 
    … …  
    11321133no effect on the eddy-induced skew-fluxes. 
    11331134 
    1134 \subsection{Limiting of the slopes within the interior}\label{sec:triad:limitskew} 
     1135\subsection{Limiting of the slopes within the interior}\label{sec:limitskew} 
    11351136Presently, the iso-neutral slopes $\tilde{r}_i$ relative 
    11361137to geopotentials are limited to be less than $1/100$, exactly as in 
    1137 calculating the iso-neutral diffusion, \S \ref{sec:triad:limit}. Each 
     1138calculating the iso-neutral diffusion, \S \autoref{sec:limit}. Each 
    11381139individual triad \rtriadt{R} is so limited. 
    11391140 
    1140 \subsection{Tapering within the surface mixed layer}\label{sec:triad:taperskew} 
     1141\subsection{Tapering within the surface mixed layer}\label{sec:taperskew} 
    11411142The slopes $\tilde{r}_i$ relative to 
    11421143geopotentials (and thus the individual triads \rtriadt{R}) are always tapered linearly from their value immediately below the mixed layer to zero at the 
    1143 surface \eqref{eq:triad:rmtilde}, as described in \S\ref{sec:triad:lintaper}. This is 
    1144 option (c) of Fig.~\ref{Fig_eiv_slp}. This linear tapering for the 
     1144surface \autoref{eq:rmtilde}, as described in \autoref{sec:lintaper}. This is 
     1145option (c) of \autoref{fig:eiv_slp}. This linear tapering for the 
    11451146slopes used to calculate the eddy-induced fluxes is 
    11461147unaffected by the value of \np{ln\_triad\_iso}. 
    … …  
    11481149The justification for this linear slope tapering is that, for $A_e$ 
    11491150that is constant or varies only in the horizontal (the most commonly 
    1150 used options in \NEMO: see \S\ref{LDF_coef}), it is 
     1151used options in \NEMO: see \autoref{sec:LDF_coef}), it is 
    11511152equivalent to a horizontal eiv (eddy-induced velocity) that is uniform 
    1152 within the mixed layer \eqref{eq:triad:eiv_v}. This ensures that the 
     1153within the mixed layer \autoref{eq:eiv_v}. This ensures that the 
    11531154eiv velocities do not restratify the mixed layer \citep{Treguier1997, 
    11541155  Danabasoglu_al_2008}. Equivantly, in terms 
    … …  
    11581159horizontal flux convergence is relatively insignificant within the mixed-layer). 
    11591160 
    1160 \subsection{Streamfunction diagnostics}\label{sec:triad:sfdiag} 
     1161\subsection{Streamfunction diagnostics}\label{sec:sfdiag} 
    11611162Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.}, diagnosed 
    11621163mean eddy-induced velocities are output. Each time step, 
    … …  
    11641165$uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ 
    11651166(integer $i$, integer +1/2 $j$, integer +1/2 $k$) points (see Table 
    1166 \ref{Tab_cell}) respectively. We follow \citep{Griffies_Bk04} and 
     1167\autoref{tab:cell}) respectively. We follow \citep{Griffies_Bk04} and 
    11671168calculate the streamfunction at a given $uw$-point from the 
    11681169surrounding four triads according to: 
    11691170\begin{equation} 
    1170   \label{eq:triad:sfdiagi} 
     1171  \label{eq:sfdiagi} 
    11711172  {\psi_1}_{i+1/2}^{k+1/2}={\quarter}\sum_{\substack{i_p,\,k_p}} 
    11721173  {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}. 
    … …  
    11741175The streamfunction $\psi_1$ is calculated similarly at $vw$ points. 
    11751176The eddy-induced velocities are then calculated from the 
    1176 straightforward discretisation of \eqref{eq:triad:eiv_v}: 
    1177 \begin{equation}\label{eq:triad:eiv_v_discrete} 
     1177straightforward discretisation of \autoref{eq:eiv_v}: 
     1178\begin{equation}\label{eq:eiv_v_discrete} 
    11781179\begin{split} 
    11791180 {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),   \\ 
Note: See TracChangeset for help on using the changeset viewer.