# source:branches/2017/dev_merge_2017/DOC/TexFiles/Chapters/Annex_ISO.tex@9364

Last change on this file since 9364 was 9364, checked in by nicolasmartin, 3 years ago

Fix mathrm syntax, equations in multline math env is still not accepted

File size: 58.2 KB
Line
1\documentclass[NEMO_book]{subfiles}
2\begin{document}
3% ================================================================
4% Iso-neutral diffusion :
5% ================================================================
6\chapter[Iso-neutral diffusion and eddy advection using
9\minitoc
10\pagebreak
11\section{Choice of \protect\ngn{namtra\_ldf} namelist parameters}
12%-----------------------------------------nam_traldf------------------------------------------------------
13\namdisplay{namtra_ldf}
14%---------------------------------------------------------------------------------------------------------
15
16Two scheme are available to perform the iso-neutral diffusion.
17If the namelist logical \np{ln\_traldf\_triad} is set true,
18\NEMO updates both active and passive tracers using the Griffies triad representation
19of iso-neutral diffusion and the eddy-induced advective skew (GM) fluxes.
20If the namelist logical \np{ln\_traldf\_iso} is set true,
21the filtered version of Cox's original scheme (the Standard scheme) is employed (\S\ref{LDF_slp}).
22In the present implementation of the Griffies scheme,
23the advective skew fluxes are implemented even if \np{ln\_traldf\_eiv} is false.
24
25Values of iso-neutral diffusivity and GM coefficient are set as
26described in \S\ref{LDF_coef}. Note that when GM fluxes are used,
27the eddy-advective (GM) velocities are output for diagnostic purposes using xIOS,
28even though the eddy advection is accomplished by means of the skew fluxes.
29
30
31The options specific to the Griffies scheme include:
32\begin{description}[font=\normalfont]
34  iso-neutral' mixing is accomplished within the surface mixed-layer
35  along slopes linearly decreasing with depth from the value immediately below
36  the mixed-layer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}).
37  This is the same treatment as used in the default implementation \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}
38  Where \np{ln\_triad\_iso} is set true, the vertical skew flux is further reduced
39  to ensure no vertical buoyancy flux, giving an almost pure
40  horizontal diffusive tracer flux within the mixed layer. This is similar to
41  the tapering suggested by \citet{Gerdes1991}. See \S\ref{sec:triad:Gerdes-taper}
43  If this is set false (the default) then the lateral diffusive fluxes
45  If it is set true, however, then these lateral diffusive fluxes are applied,
46  giving smoother bottom tracer fields at the cost of introducing diapycnal mixing.
48\end{description}
49The options shared with the Standard scheme include:
50\begin{description}[font=\normalfont]
51\item[\np{ln\_traldf\_msc}]   blah blah to be added
52\item[\np{rn\_slpmax}]  blah blah to be added
53\end{description}
56We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98},
57but formulated within the \NEMO framework, using scale factors rather than grid-sizes.
58
59\subsection{The iso-neutral diffusion operator}
60The iso-neutral second order tracer diffusive operator for small
61angles between iso-neutral surfaces and geopotentials is given by
62\eqref{Eq_PE_iso_tensor}:
64  \begin{equation}
65    D^{lT}=-\Div\vect{f}^{lT}\equiv
66    -\frac{1}{e_1e_2e_3}\left[\pd{i}\left (f_1^{lT}e_2e_3\right) +
67      \pd{j}\left (f_2^{lT}e_2e_3\right) + \pd{k}\left (f_3^{lT}e_1e_2\right)\right],
68  \end{equation}
69  where the diffusive flux per unit area of physical space
70  \begin{equation}
72  \end{equation}
73  \begin{equation}
76    \begin{pmatrix}
77       1   &  0   & -r_1           \mystrut \\
78       0   &  1   & -r_2           \mystrut \\
79      -r_1 & -r_2 &  r_1 ^2+r_2 ^2 \mystrut
80    \end{pmatrix}
82    \begin{pmatrix}
83      \frac{1}{e_1} \pd[T]{i} \mystrut \\
84      \frac{1}{e_2} \pd[T]{j} \mystrut \\
85      \frac{1}{e_3} \pd[T]{k} \mystrut
86    \end{pmatrix}.
87  \end{equation}
88\end{subequations}
89% \left( {{\begin{array}{*{20}c}
90%  1 \hfill & 0 \hfill & {-r_1 } \hfill \\
91%  0 \hfill & 1 \hfill & {-r_2 } \hfill \\
92%  {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\
93% \end{array} }} \right)
94 Here \eqref{Eq_PE_iso_slopes}
95\begin{align*}
96  r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i}
97  \right)
98  \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \\
99  &=-\frac{e_3 }{e_1 } \left( -\alpha\frac{\partial T }{\partial i} +
100    \beta\frac{\partial S }{\partial i} \right) \left(
101    -\alpha\frac{\partial T }{\partial k} + \beta\frac{\partial S
102    }{\partial k} \right)^{-1}
103\end{align*}
104is the $i$-component of the slope of the iso-neutral surface relative to the computational
105surface, and $r_2$ is the $j$-component.
106
107We will find it useful to consider the fluxes per unit area in $i,j,k$
108space; we write
109\begin{equation}
111  \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right).
112\end{equation}
113Additionally, we will sometimes write the contributions towards the
114fluxes $\vect{f}$ and $\vect{F}_{\mathrm{iso}}$ from the component
115$R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$, with
116$f_{ij}=R_{ij}e_i^{-1}\partial T/\partial x_i$ (no summation) etc.
117
118The 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$i$- and $j$-directions resulting from the vertical tracer gradient:
121\begin{align}
123  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}\\
124\intertext{and in the k-direction resulting from the lateral tracer gradients}
126 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}
127\end{align}
128
129The vertical diffusive flux associated with the $_{33}$
130component of the small angle diffusion tensor is
131\begin{equation}
133  f_{33}=-\Alt(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}.
134\end{equation}
135
136Since there are no cross terms involving $r_1$ and $r_2$ in the above, we can
137consider the iso-neutral diffusive fluxes separately in the $i$-$k$ and $j$-$k$
138planes, just adding together the vertical components from each
139plane. The following description will describe the fluxes on the $i$-$k$
140plane.
141
142There is no natural discretization for the $i$-component of the
144although it must be evaluated at $u$-points, it involves vertical
145gradients (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 but involves horizontal gradients defined at $u$-points.
148
149\subsection{The standard discretization}
150The 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
152into OPA, \eqref{Eq_tra_ldf_iso}, is to calculate a mean vertical
153gradient at the $u$-point from the average of the four surrounding
154vertical tracer gradients, and multiply this by a mean slope at the
155$u$-point, calculated from the averaged surrounding vertical density
156gradients. The total area-integrated skew-flux (flux per unit area in
157$ijk$ space) from tracer cell $i,k$
158to $i+1,k$, noting that the $e_{{3}_{i+1/2}^k}$ in the area
159$e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with
160the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer
162\begin{equation*}
163  \left(F_u^{13} \right)_{i+\hhalf}^k = \Alts_{i+\hhalf}^k
164  {e_{2}}_{i+1/2}^k \overline{\overline
165    r_1} ^{\,i,k}\,\overline{\overline{\delta_k T}}^{\,i,k},
166\end{equation*}
167where
168\begin{equation*}
169  \overline{\overline
170   r_1} ^{\,i,k} = -\frac{{e_{3u}}_{i+1/2}^k}{{e_{1u}}_{i+1/2}^k}
171  \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}},
172\end{equation*}
173and here and in the following we drop the $^{lT}$ superscript from
174$\Alt$ for simplicity.
175Unfortunately the resulting combination $\overline{\overline{\delta_k 176 \bullet}}^{\,i,k}$ of a $k$ average and a $k$ difference %of the tracer
177reduces to $\bullet_{k+1}-\bullet_{k-1}$, so two-grid-point oscillations are
178invisible to this discretization of the iso-neutral operator. These
179\emph{computational modes} will not be damped by this operator, and
180may even possibly be amplified by it.  Consequently, applying this
181operator to a tracer does not guarantee the decrease of its
182global-average variance. To correct this, we introduced a smoothing of
183the slopes of the iso-neutral surfaces (see \S\ref{LDF}). This
184technique works for $T$ and $S$ in so far as they are active tracers
185($i.e.$ they enter the computation of density), but it does not work
186for a passive tracer.
187\subsection{Expression of the skew-flux in terms of triad slopes}
188\citep{Griffies_al_JPO98} introduce a different discretization of the
189off-diagonal terms that nicely solves the problem.
190% Instead of multiplying the mean slope calculated at the $u$-point by
191% the mean vertical gradient at the $u$-point,
192% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
193\begin{figure}[tb] \begin{center}
196      (a) Arrangement of triads $S_i$ and tracer gradients to
197           give lateral tracer flux from box $i,k$ to $i+1,k$
198      (b) Triads $S'_i$ and tracer gradients to give vertical tracer flux from
199            box $i,k$ to $i,k+1$.}
200 \end{center} \end{figure}
201% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
202They get the skew flux from the products of the vertical gradients at
203each $w$-point surrounding the $u$-point with the corresponding triad'
204slope calculated from the lateral density gradient across the $u$-point
205divided by the vertical density gradient at the same $w$-point as the
207denote the tracer gradients, and the thin lines the corresponding
208triads, with slopes $s_1, \dotsc s_4$. The total area-integrated
209skew-flux from tracer cell $i,k$ to $i+1,k$
210\begin{multline}
212  \left( F_u^{13}  \right)_{i+\frac{1}{2}}^k = \Alts_{i+1}^k a_1 s_1
213  \delta _{k+\frac{1}{2}} \left[ T^{i+1}
214  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}  + \Alts _i^k a_2 s_2 \delta
215  _{k+\frac{1}{2}} \left[ T^i
216  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} \\
217   +\Alts _{i+1}^k a_3 s_3 \delta _{k-\frac{1}{2}} \left[ T^{i+1}
218  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}  +\Alts _i^k a_4 s_4 \delta
219  _{k-\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}},
220\end{multline}
221where the contributions of the triad fluxes are weighted by areas
222$a_1, \dotsc a_4$, and $\Alts$ is now defined at the tracer points
223rather than the $u$-points. This discretization gives a much closer
224stencil, and disallows the two-point computational modes.
225
226 The vertical skew flux \eqref{eq:triad:i31c} from tracer cell $i,k$ to $i,k+1$ at the
227$w$-point $i,k+\hhalf$ is constructed similarly (Fig.~\ref{fig:triad:ISO_triad}b)
228by multiplying lateral tracer gradients from each of the four
229surrounding $u$-points by the appropriate triad slope:
230\begin{multline}
232  \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} =  \Alts_i^{k+1} a_{1}'
233  s_{1}' \delta _{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1}
234   +\Alts_i^{k+1} a_{2}' s_{2}' \delta _{i+\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i+\frac{1}{2}}^{k+1}\\
235  + \Alts_i^k a_{3}' s_{3}' \delta _{i-\frac{1}{2}} \left[ T^k\right]/{e_{3u}}_{i-\frac{1}{2}}^k
236  +\Alts_i^k a_{4}' s_{4}' \delta _{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k.
237\end{multline}
238
239We notate the triad slopes $s_i$ and $s'_i$ in terms of the anchor point' $i,k$
240(appearing in both the vertical and lateral gradient), and the $u$- and
241$w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the arms' of the
243\begin{equation}
245  _i^k \mathbb{R}_{i_p}^{k_p}
246  =-\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}}
247  \
248  \frac
249  { \alpha_i^\ \delta_{i+i_p}[T^k] - \beta_i^k \ \delta_{i+i_p}[S^k] }
250  { \alpha_i^\ \delta_{k+k_p}[T^i] - \beta_i^k \ \delta_{k+k_p}[S^i] }.
251\end{equation}
252In calculating the slopes of the local neutral surfaces,
253the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad,
254while the metrics are calculated at the $u$- and $w$-points on the arms.
255
256% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
257\begin{figure}[tb] \begin{center}
258    \includegraphics[width=0.80\textwidth]{Fig_GRIFF_qcells}
260    Triad notation for quarter cells. $T$-cells are inside
261      boxes, while the  $i+\half,k$ $u$-cell is shaded in green and the
262      $i,k+\half$ $w$-cell is shaded in pink.}
263  \end{center} \end{figure}
264% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
265
266Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter
267cell 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.
268Expressing the slopes $s_i$ and $s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation,
269we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$.
270Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$)
271to calculate the lateral flux along its $u$-arm, at $(i+i_p,k)$,
272and then again as an $s'$ to calculate the vertical flux along its $w$-arm at $(i,k+k_p)$.
273Each vertical area $a_i$ used to calculate the lateral flux and horizontal area $a'_i$ used
274to calculate the vertical flux can also be identified as the area across the $u$- and $w$-arms
275of a unique triad, and we notate these areas, similarly to the triad slopes,
276as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$,
277where $e.g.$ in \eqref{eq:triad:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,
278and in \eqref{eq:triad:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.
279
281A key property of iso-neutral diffusion is that it should not affect
282the (locally referenced) density. In particular there should be no
283lateral or vertical density flux. The lateral density flux disappears so long as the
284area-integrated lateral diffusive flux from tracer cell $i,k$ to
285$i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the
286form
287\begin{equation}
289  \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} =
290  - \left( \Alts_i^{k+1} a_{1} + \Alts_i^{k+1} a_{2} + \Alts_i^k
291    a_{3} + \Alts_i^k a_{4} \right)
292  \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}},
293\end{equation}
294where the areas $a_i$ are as in \eqref{eq:triad:i13}. In this case,
295separating the total lateral flux, the sum of \eqref{eq:triad:i13} and
297flux
298\begin{equation}
300  _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - \Alts_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p}
301  \left(
302    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
303    -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \
304    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
305  \right)
306\end{equation}
307can be identified with each triad. Then, because the
308same metric factors ${e_{3w}}_{\,i}^{\,k+k_p}$ and
309${e_{1u}}_{\,i+i_p}^{\,k}$ are employed for both the density gradients
310in $_i^k \mathbb{R}_{i_p}^{k_p}$ and the tracer gradients, the lateral
311density flux associated with each triad separately disappears.
312\begin{equation}
314  {\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
315\end{equation}
316Thus the total flux $\left( F_u^{31} \right) ^i _{i,k+\frac{1}{2}} + 317\left( F_u^{11} \right) ^i _{i,k+\frac{1}{2}}$ from tracer cell $i,k$
318to $i+1,k$ must also vanish since it is a sum of four such triad fluxes.
319
320The squared slope $r_1^2$ in the expression \eqref{eq:triad:i33c} for the
321$_{33}$ component is also expressed in terms of area-weighted
322squared triad slopes, so the area-integrated vertical flux from tracer
323cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is
324\begin{equation}
326  \left( F_w^{33} \right) _i^{k+\frac{1}{2}} =
327    - \left( \Alts_i^{k+1} a_{1}' s_{1}'^2
328    + \Alts_i^{k+1} a_{2}' s_{2}'^2
329    + \Alts_i^k a_{3}' s_{3}'^2
330    + \Alts_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right],
331\end{equation}
332where the areas $a'$ and slopes $s'$ are the same as in
334Then, separating the total vertical flux, the sum of \eqref{eq:triad:i31} and
336\begin{align}
338  _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T)
339  &= \Alts_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p}
340  \left(
341    {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
342    -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \
343    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
344  \right) \\
345  &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right)
347\end{align}
348may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$
349associated with a triad then separately disappears (because the
350lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (\rho)$
351disappears). Consequently the total vertical density flux $\left( F_w^{31} \right)_i ^{k+\frac{1}{2}} + 352\left( F_w^{33} \right)_i^{k+\frac{1}{2}}$ from tracer cell $i,k$
353to $i,k+1$ must also vanish since it is a sum of four such triad
354fluxes.
355
356We 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
357the $u$-fluxes and $w$-fluxes in
359Fig.~\ref{fig:triad:ISO_triad} to  write out the iso-neutral fluxes at $u$- and
360$w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces:
362\begin{flalign} \label{Eq_iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv
363  \sum_{\substack{i_p,\,k_p}}
364  \begin{pmatrix}
365    {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\
366    \\
367    {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)      \\
368  \end{pmatrix}.
369\end{flalign}
370\subsection{Ensuring the scheme does not increase tracer variance}
372
373We now require that this operator should not increase the
374globally-integrated tracer variance.
375%This changes according to
376% \begin{align*}
377% &\int_D  D_l^T \; T \;dv \equiv  \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\}    \\
378% &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
379%     \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
380%       + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]  \ T \right\}    \\
381% &\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
382%                 {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T]
383%              + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\
384% \end{align*}
385Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux
386$_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the $u$-point $i+i_p,k$ and
387a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the
388$w$-point $i,k+k_p$.  The lateral flux drives a net rate of change of
389variance, summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
390\begin{multline}
391  {b_T}_{i+i_p-1/2}^k\left(\frac{\partial T}{\partial t}T\right)_{i+i_p-1/2}^k+
393      t}T\right)_{i+i_p+1/2}^k \\
394 \begin{split}
396  {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\
397  &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:triad:dvar_iso_i}
398 \end{split}
399\end{multline}
400while the vertical flux similarly drives a net rate of change of
401variance summed over the $T$-points $i,k+k_p-\half$ (above) and
402$i,k+k_p+\half$ (below) of
403\begin{equation}
405  _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
406\end{equation}
407The total variance tendency driven by the triad is the sum of these
408two. Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and
409$_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \eqref{eq:triad:latflux-triad} and
411\begin{multline*}
412-\Alts_i^k\left \{
413{ } _i^k{\mathbb{A}_u}_{i_p}^{k_p}
414  \left(
415    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
416    - {_i^k\mathbb{R}_{i_p}^{k_p}} \
417    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }\right)\,\delta_{i+ i_p}[T^k] \right.\\
418- \left. { } _i^k{\mathbb{A}_w}_{i_p}^{k_p}
419  \left(
420    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
421    -{\:}_i^k\mathbb{R}_{i_p}^{k_p}
422    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
423  \right) {\,}_i^k\mathbb{R}_{i_p}^{k_p}\delta_{k+ k_p}[T^i]
424\right \}.
425\end{multline*}
426The key point is then that if we require
427$_i^k{\mathbb{A}_u}_{i_p}^{k_p}$ and $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$
428to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by
429\begin{equation}
431  _i^k\mathbb{V}_{i_p}^{k_p}
432  ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k}
433  ={\;}_i^k{\mathbb{A}_w}_{i_p}^{k_p}\,{e_{3w}}_{\,i}^{\,k + k_p},
434\end{equation}
435the variance tendency reduces to the perfect square
436\begin{equation}
438  -\Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
439  \left(
440    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
441    -{\:}_i^k\mathbb{R}_{i_p}^{k_p}
442    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
443  \right)^2\leq 0.
444\end{equation}
446with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase
447the net variance. Since the total fluxes are sums of such fluxes from
449sufficient to ensure that the globally integrated variance does not
450increase.
451
452The expression \eqref{eq:triad:V-A} can be interpreted as a discretization
453of the global integral
454\begin{equation}
456  \frac{\partial}{\partial t}\int\!\half T^2\, dV =
457  \int\!\mathbf{F}\cdot\nabla T\, dV,
458\end{equation}
459where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the
460lateral and vertical fluxes/unit area
461$462\mathbf{F}=\left( 463\left.{}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_u}_{i_p}^{k_p}, 464\left.{\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_w}_{i_p}^{k_p} 465 \right) 466$
468 $\nabla T = \left( 469\left.\delta_{i+ i_p}[T^k] \right/ {e_{1u}}_{\,i + i_p}^{\,k}, 470\left.\delta_{k+ k_p}[T^i] \right/ {e_{3w}}_{\,i}^{\,k + k_p} 471\right) 472$
473\subsection{Triad volumes in Griffes's scheme and in \NEMO}
474To complete the discretization we now need only specify the triad
475volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. \citet{Griffies_al_JPO98} identify
476these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter
477cells, defined in terms of the distances between $T$, $u$,$f$ and
478$w$-points. This is the natural discretization of
479\eqref{eq:triad:cts-var}. The \NEMO model, however, operates with scale
480factors instead of grid sizes, and scale factors for the quarter
481cells are not defined. Instead, therefore we simply choose
482\begin{equation}
484  _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k,
485\end{equation}
486as a quarter of the volume of the $u$-cell inside which the triad
487quarter-cell lies. This has the nice property that when the slopes
488$\mathbb{R}$ vanish, the lateral flux from tracer cell $i,k$ to
489$i+1,k$ reduces to the classical form
490\begin{equation}
492-\overline\Alts_{\,i+1/2}^k\;
493\frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}}
494\;\frac{\delta_{i+ 1/2}[T^k] }{{e_{1u}}_{\,i + i_p}^{\,k}}
495 = -\overline\Alts_{\,i+1/2}^k\;\frac{{e_{1w}}_{\,i + 1/2}^{\,k}\:{e_{1v}}_{\,i + 1/2}^{\,k}\;\delta_{i+ 1/2}[T^k]}{{e_{1u}}_{\,i + 1/2}^{\,k}}.
496\end{equation}
497In fact if the diffusive coefficient is defined at $u$-points, so that
498we employ $\Alts_{i+i_p}^k$ instead of  $\Alts_i^k$ in the definitions of the
500we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above.
501
502\subsection{Summary of the scheme}
503The iso-neutral fluxes at $u$- and
504$w$-points are the sums of the triad fluxes that cross the $u$- and
505$w$-faces \eqref{Eq_iso_flux}:
508    \vect{F}_{\mathrm{iso}}(T) &\equiv
509    \sum_{\substack{i_p,\,k_p}}
510    \begin{pmatrix}
511      {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\
512      \\
513      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)
514    \end{pmatrix},
515  \end{flalign}
517  \begin{align}
519    _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - \Alts_i^k{
520      \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}}
521    \left(
522      \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
523      -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \
524      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
525    \right),\\
526    \intertext{and}
527    _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T)
528    &= \Alts_i^k{\: }\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{3w}}_{\,i}^{\,k+k_p}}
529    \left(
530      {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
531      -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \
532      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
534  \end{align}
536  \begin{equation}
538    _i^k{\mathbb{V}}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k.
539  \end{equation}
540\end{subequations}
541
542 The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at
543each tracer point:
545  \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k
546        {\mathbb{F}_u }_{i_p}^{k_p}} \right] + \delta_{k} \left[
547      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \right\}
548\end{equation}
549where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells.
550The diffusion scheme satisfies the following six properties:
551\begin{description}
552\item[$\bullet$ horizontal diffusion] The discretization of the
554  the limit of flat iso-neutral direction :
555  \begin{equation} \label{eq:triad:iso_property0} D_l^T = \frac{1}{b_T} \
556    \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \;
557      \overline\Alts^{\,i} \; \delta_{i+1/2}[T] \right] \qquad
558    \text{when} \quad { _i^k \mathbb{R}_{i_p}^{k_p} }=0
559  \end{equation}
560
561\item[$\bullet$ implicit treatment in the vertical] Only tracer values
562  associated with a single water column appear in the expression
563  \eqref{eq:triad:i33} for the $_{33}$ fluxes, vertical fluxes driven by
564  vertical gradients. This is of paramount importance since it means
565  that a time-implicit algorithm can be used to solve the vertical
566  diffusion equation. This is necessary
567 since the vertical eddy
568  diffusivity associated with this term,
569  \begin{equation}
570    \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{
571      {\:}_i^k\mathbb{V}_{i_p}^{k_p} \: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
572    \right\}  =
573    \frac{1}{4b_w}\sum_{\substack{i_p, \,k_p}} \left\{
574      {b_u}_{i+i_p}^k\: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
575    \right\},
576 \end{equation}
577  (where $b_w= e_{1w}\,e_{2w}\,e_{3w}$ is the volume of $w$-cells) can be quite large.
578
579\item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of
580  locally referenced potential density is zero. See
582
583\item[$\bullet$ conservation of tracer] The iso-neutral diffusion
584  conserves tracer content, $i.e.$
585  \begin{equation} \label{eq:triad:iso_property1} \sum_{i,j,k} \left\{ D_l^T \
586      b_T \right\} = 0
587  \end{equation}
588  This property is trivially satisfied since the iso-neutral diffusive
589  operator is written in flux form.
590
591\item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion
592  does not increase the tracer variance, $i.e.$
593  \begin{equation} \label{eq:triad:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T
594      \ b_T \right\} \leq 0
595  \end{equation}
596  The property is demonstrated in
597  \S\ref{sec:triad:variance} above. It is a key property for a diffusion
598  term. It means that it is also a dissipation term, $i.e.$ it
599  dissipates the square of the quantity on which it is applied.  It
600  therefore ensures that, when the diffusivity coefficient is large
601  enough, the field on which it is applied becomes free of grid-point
602  noise.
603
604\item[$\bullet$ self-adjoint operator] The iso-neutral diffusion
605  operator is self-adjoint, $i.e.$
606  \begin{equation} \label{eq:triad:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T
607      \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}
608  \end{equation}
609  In other word, there is no need to develop a specific routine from
610  the adjoint of this operator. We just have to apply the same
611  routine. This property can be demonstrated similarly to the proof of
612  the no increase of tracer variance' property. The contribution by a
614  be found by replacing $\delta[T]$ by $\delta[S]$ in \eqref{eq:triad:dvar_iso_i}
615  and \eqref{eq:triad:dvar_iso_k}. This results in a term similar to
617\begin{equation}
619  - \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
620  \left(
621    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
622    -{\:}_i^k\mathbb{R}_{i_p}^{k_p}
623    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
624  \right)
625  \left(
626    \frac{ \delta_{i+ i_p}[S^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
627    -{\:}_i^k\mathbb{R}_{i_p}^{k_p}
628    \frac{ \delta_{k+k_p} [S^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
629  \right).
630\end{equation}
631This is symmetrical in $T$ and $S$, so exactly the same term arises
632from the discretization of this triad's contribution towards the
634\end{description}
636The triad slope can only be defined where both the grid boxes centred at
637the end of the arms exist. Triads that would poke up
638through the upper ocean surface into the atmosphere, or down into the
640$\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and
641$\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be
643tracer gradients produce no flux through the ocean surface. However,
644to prevent surface noise, it is customary to retain the $_{11}$ contributions towards
645the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and
646$\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer
649triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and
650$\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$
651or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is
652masked. The associated lateral fluxes (grey-black dashed line) are
655
656The default option \np{ln\_botmix\_triad}=false is suitable when the
657bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}=1),
658or  for simple idealized  problems. For setups with topography without
659bbl mixing, \np{ln\_botmix\_triad}=true may be necessary.
660% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
661\begin{figure}[h] \begin{center}
664      (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer
665      points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad
666      slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$
668      figure). However, the lateral $_{11}$ contributions towards
669      $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$
670      (yellow line) are still applied, giving diapycnal diffusive
671      fluxes.\newline
672      (b) Both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and
673      $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$
674      or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point
675      is masked. The associated lateral fluxes (grey-black dashed
678 \end{center} \end{figure}
679% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
680\subsection{ Limiting of the slopes within the interior}\label{sec:triad:limit}
681As discussed in \S\ref{LDF_slp_iso}, iso-neutral slopes relative to
682geopotentials must be bounded everywhere, both for consistency with the small-slope
683approximation and for numerical stability \citep{Cox1987,
684  Griffies_Bk04}. The bound chosen in \NEMO is applied to each
685component of the slope separately and has a value of $1/100$ in the ocean interior.
686%, ramping linearly down above 70~m depth to zero at the surface
687It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to
688geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to
689geopotentials) \eqref{Eq_PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate
690surfaces, so we require
691\begin{equation*}
692  |\tilde{r}_i|\leq \tilde{r}_\mathrm{max}=0.01.
693\end{equation*}
694and then recalculate the slopes $r_i$ relative to coordinates.
696 \begin{equation}
698_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}}
699 \end{equation}
700is limited like this and then the corresponding
701$_i^k\mathbb{R}_{i_p}^{k_p}$ are recalculated and combined to form the fluxes.
702Note that where the slopes have been limited, there is now a non-zero
703iso-neutral density flux that drives dianeutral mixing.  In particular this iso-neutral density flux
704is always downwards, and so acts to reduce gravitational potential energy.
705\subsection{Tapering within the surface mixed layer}\label{sec:triad:taper}
706
707Additional tapering of the iso-neutral fluxes is necessary within the
708surface mixed layer. When the Griffies triads are used, we offer two
709options for this.
710\subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper}
711This is the option activated by the default choice
712\np{ln\_triad\_iso}=false. Slopes $\tilde{r}_i$ relative to
713geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the
714surface, as described in option (c) of Fig.~\ref{Fig_eiv_slp}, to values
715\begin{subequations}
716  \begin{equation}
718     \rMLt =
719  -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h,
720  \end{equation}
721and then the $r_i$ relative to vertical coordinate surfaces are appropriately
723  \begin{equation}
725 \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h.
726  \end{equation}
727\end{subequations}
728Thus the diffusion operator within the mixed layer is given by:
730D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
732 1 \hfill & 0 \hfill & {-\rML[1]}\hfill \\
733 0 \hfill & 1 \hfill & {-\rML[2]} \hfill \\
734 {-\rML[1]}\hfill &   {-\rML[2]} \hfill & {\rML[1]^2+\rML[2]^2} \hfill
735\end{array} }} \right)
736\end{equation}
737
738This slope tapering gives a natural connection between tracer in the
739mixed-layer and in isopycnal layers immediately below, in the
740thermocline. It is consistent with the way the $\tilde{r}_i$ are
741tapered within the mixed layer (see \S\ref{sec:triad:taperskew} below)
742so as to ensure a uniform GM eddy-induced velocity throughout the
743mixed layer. However, it gives a downwards density flux and so acts so
744as to reduce potential energy in the same way as does the slope
746
747As in \S\ref{sec:triad:limit} above, the tapering
749$_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the
750$_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume
751$z$-coordinates in the following; the conversion from
752$\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described
754\begin{enumerate}
755\item Mixed-layer depth is defined so as to avoid including regions of weak
756vertical stratification in the slope definition.
757 At each $i,j$ (simplified to $i$ in
759the vertical index of the tracer point immediately below the mixed
760layer, $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point)
761such that the potential density
762${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, where $i,k_{10}$ is
763the tracer gridbox within which the depth reaches 10~m. See the left
764side of Fig.~\ref{fig:triad:MLB_triad}. We use the $k_{10}$-gridbox
765instead of the surface gridbox to avoid problems e.g.\ with thin
766daytime mixed-layers. Currently we use the same
767$\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is
768used to output the diagnosed mixed-layer depth
769$h_{\mathrm{ML}}=|z_{W}|_{k_{\mathrm{ML}}+1/2}$, the depth of the $w$-point
770above the $i,k_{\mathrm{ML}}$ tracer point.
771
772\item We define basal' triad slopes
773${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ as the slopes
774of those triads whose vertical arms' go down from the
775$i,k_{\mathrm{ML}}$ tracer point to the $i,k_{\mathrm{ML}}-1$ tracer point
776below. This is to ensure that the vertical density gradients
777associated with these basal triad slopes
778${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are
779representative of the thermocline. The four basal triads defined in the bottom part
781\begin{align}
782  {\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p} &=
784\\
786{\:}_i{\mathbb{R}_{\mathrm{base}}}_{1/2}^{-1/2}&=
787{\:}^{k_{\mathrm{ML}}}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2}. \notag
788\end{align}
789The vertical flux associated with each of these triads passes through the $w$-point
790$i,k_{\mathrm{ML}}-1/2$ lying \emph{below} the $i,k_{\mathrm{ML}}$ tracer point,
791so it is this depth
792\begin{equation}
794  {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2}
795\end{equation}
796(one gridbox deeper than the
797diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper
800${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within the mixed
801layer, by multiplying the appropriate
802${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ by the ratio of
803the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_{\mathrm{base}}}_{\,i}$. For
804instance the green triad centred on $i,k$
805\begin{align}
806  {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,1/2}^{-1/2} &=
807\frac{{z_w}_{k-1/2}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2}
808\notag \\
809\intertext{and more generally}
810 {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p} &=
812\end{align}
813\end{enumerate}
814
815% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
816\begin{figure}[h]
817%  \fcapside {
819      mixed-layer depth and calculation of linearly tapered
820      triads. The figure shows a water column at a given $i,j$
821      (simplified to $i$), with the ocean surface at the top. Tracer points are denoted by
822      bullets, and black lines the edges of the tracer cells; $k$
823      increases upwards. \newline
824      \hspace{5 em}We define the mixed-layer by setting the vertical index
825      of the tracer point immediately below the mixed layer,
826      $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point)
827      such that ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$,
828      where $i,k_{10}$ is the tracer gridbox within which the depth
829      reaches 10~m. We calculate the triad slopes within the mixed
830      layer by linearly tapering them from zero (at the surface) to
831      the basal' slopes, the slopes of the four triads passing through the
832      $w$-point $i,k_{\mathrm{ML}}-1/2$ (blue square),
833      ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$. Triads with
834    different $i_p,k_p$, denoted by different colours, (e.g. the green
835    triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.}
836%}
838\end{figure}
839% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
840
841\subsubsection{Additional truncation of skew iso-neutral flux
842  components}
844The alternative option is activated by setting \np{ln\_triad\_iso} =
845  true. This retains the same tapered slope $\rML$  described above for the
846calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the
847vertical tracer flux driven by vertical tracer gradients), but
848replaces the $\rML$ in the skew term by
849\begin{equation}
851  \rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i,
852\end{equation}
853giving a ML diffusive operator
855D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
857 1 \hfill & 0 \hfill & {-\rML[1]^*}\hfill \\
858 0 \hfill & 1 \hfill & {-\rML[2]^*} \hfill \\
859 {-\rML[1]^*}\hfill &   {-\rML[2]^*} \hfill & {\rML[1]^2+\rML[2]^2} \hfill \\
860\end{array} }} \right).
861\end{equation}
862This operator
863\footnote{To ensure good behaviour where horizontal density
865$\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.}
866then has the property it gives no vertical density flux, and so does
867not change the potential energy.
868This approach is similar to multiplying the iso-neutral  diffusion
869coefficient by $\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep
871Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$
872
873In practice, this approach gives weak vertical tracer fluxes through
874the mixed-layer, as well as vanishing density fluxes. While it is
875theoretically advantageous that it does not change the potential
876energy, it may give a discontinuity between the
877fluxes within the mixed-layer (purely horizontal) and just below (along
878iso-neutral surfaces).
879% This may give strange looking results,
880% particularly where the mixed-layer depth varies strongly laterally.
881% ================================================================
882% Skew flux formulation for Eddy Induced Velocity :
883% ================================================================
885
887
888 When Gent and McWilliams's [1990] diffusion is used,
890eddy induced velocity, the formulation of which depends on the slopes of iso-
891neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used
892here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo}
893is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo}
894+ \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates.
895
896The eddy induced velocity is given by:
899\begin{split}
900 u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\
901 v^* & = - \frac{1}{e_{3}}\;          \partial_j\psi_2,    \\
902w^* & =    \frac{1}{e_{1}e_{2}}\; \left\{ \partial_\left( e_{2} \, \psi_1\right)
903                         + \partial_\left( e_{1} \, \psi_2\right) \right\},
904\end{split}
905\end{equation}
906where the streamfunctions $\psi_i$ are given by
908\begin{split}
909\psi_1 & = A_{e} \; \tilde{r}_1,   \\
910\psi_2 & = A_{e} \; \tilde{r}_2,
911\end{split}
912\end{equation}
913\end{subequations}
914with $A_{e}$ the eddy induced velocity coefficient, and $\tilde{r}_1$ and $\tilde{r}_2$ the slopes between the iso-neutral and the geopotential surfaces.
915
917it to the Eulerian velocity prior to computing the tracer
918advection. This is implemented if \key{traldf\_eiv} is set in the
919default implementation, where \np{ln\_traldf\_triad} is set
920false. This allows us to take advantage of all the advection schemes
921offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$
922order advection scheme. This is particularly useful for passive
923tracers where \emph{positivity} of the advection scheme is of
924paramount importance.
925
927implements eddy induced advection according to the so-called skew form
928\citep{Griffies_JPO98}. It is based on a transformation of the advective fluxes
929using the non-divergent nature of the eddy induced velocity.
930For example in the (\textbf{i},\textbf{k}) plane, the tracer advective
931fluxes per unit area in $ijk$ space can be
932transformed as follows:
933\begin{flalign*}
934\begin{split}
935\textbf{F}_{\mathrm{eiv}}^T =
936\begin{pmatrix}
937           {e_{2}\,e_{3}\;  u^*}       \\
938      {e_{1}\,e_{2}\; w^*}  \\
939\end{pmatrix}   \;   T
940&=
941\begin{pmatrix}
942           { - \partial_k \left( e_{2} \,\psi_1 \right) \; T \;}     \\
943      {+ \partial_\left( e_{2} \, \psi_1 \right) \; T \;}    \\
944\end{pmatrix}        \\
945&=
946\begin{pmatrix}
947           { - \partial_k \left( e_{2} \, \psi_\; T \right) \;}  \\
948      {+ \partial_\left( e_{2} \,\psi_1 \; T \right) \;}  \\
949\end{pmatrix}
950 +
951\begin{pmatrix}
952           {+ e_{2} \, \psi_\; \partial_k T}  \\
953      { - e_{2} \, \psi_\; \partial_i  T}  \\
954\end{pmatrix}
955\end{split}
956\end{flalign*}
957and since the eddy induced velocity field is non-divergent, we end up with the skew
958form of the eddy induced advective fluxes per unit area in $ijk$ space:
960\textbf{F}_\mathrm{eiv}^T = \begin{pmatrix}
961           {+ e_{2} \, \psi_\; \partial_k T}   \\
962      { - e_{2} \, \psi_\; \partial_i  T}  \\
963                                 \end{pmatrix}
964\end{equation}
965The total fluxes per unit physical area are then
967\begin{split}
968 f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T   \\
969 f^*_2 & = \frac{1}{e_{3}}\; \psi_2 \partial_k T   \\
970 f^*_3 & =  -\frac{1}{e_{1}e_{2}}\; \left\{ e_{2} \psi_1 \partial_i T
971   + e_{1} \psi_2 \partial_j T \right\}. \\
972\end{split}
973\end{equation}
974Note that Eq.~ \eqref{eq:triad:eiv_skew_physical} takes the same form whatever the
975vertical coordinate, though of course the slopes
976$\tilde{r}_i$ which define the $\psi_i$ in \eqref{eq:triad:eiv_psi} are relative to geopotentials.
977The tendency associated with eddy induced velocity is then simply the convergence
980\frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[
981  \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right)
982  + \frac{\partial}{\partial j} \left( e_\;
983    \psi_2 \partial_k T\right)
984 -  \frac{\partial}{\partial k} \left( e_{2} \psi_1 \partial_i T
985   + e_{1} \psi_2 \partial_j T \right\right]
986\end{equation}
987 It naturally conserves the tracer content, as it is expressed in flux
988 form. Since it has the same divergence as the advective form it also
989 preserves the tracer variance.
990
991\subsection{The discrete skew flux formulation}
996$\tilde{\mathbb{R}}$ relative to geopotentials instead of the
997$\mathbb{R}$ relative to coordinate surfaces. The discrete form of
999defining $A_e$ at $T$-points is then given by:
1000
1001
1004    \vect{F}_{\mathrm{eiv}}(T) &\equiv
1005    \sum_{\substack{i_p,\,k_p}}
1006    \begin{pmatrix}
1007      {_{i+1/2-i_p}^k {\mathbb{S}_u}_{i_p}^{k_p} } (T)      \\
1008      \\
1009      {_i^{k+1/2-k_p} {\mathbb{S}_w}_{i_p}^{k_p} } (T)      \\
1010    \end{pmatrix},
1011  \end{flalign}
1012  where the skew flux in the $i$-direction associated with a given
1014  \begin{align}
1016    _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \quarter {A_e}_i^k{
1017      \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}}
1018     \ {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} \
1019      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} },
1020   \\
1021    \intertext{and \eqref{eq:triad:triadfluxw} in the $k$-direction, changing the sign
1022      to be consistent with \eqref{eq:triad:eiv_skew_ijk}:}
1023    _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T)
1024    &= -\quarter {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}}
1025     {_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}
1026  \end{align}
1027\end{subequations}
1028
1029Such a discretisation is consistent with the iso-neutral
1030operator as it uses the same definition for the slopes.  It also
1031ensures the following two key properties.
1032\subsubsection{No change in tracer variance}
1033The discretization conserves tracer variance, $i.e.$ it does not
1034include a diffusive component but is a `pure' advection term. This can
1035be seen
1036%either from Appendix \ref{Apdx_eiv_skew} or
1037by considering the
1038fluxes associated with a given triad slope
1039$_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. For, following
1041associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$
1042drives a net rate of change of variance, summed over the two
1043$T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
1044\begin{equation}
1046  _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k],
1047\end{equation}
1048while the associated vertical skew-flux gives a variance change summed over the
1049$T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of
1050\begin{equation}
1052  _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
1053\end{equation}
1056sum to zero. Hence the two fluxes associated with each triad make no
1057net contribution to the variance budget.
1058
1059\subsubsection{Reduction in gravitational PE}
1060The vertical density flux associated with the vertical skew-flux
1061always has the same sign as the vertical density gradient; thus, so
1062long as the fluid is stable (the vertical density gradient is
1063negative) the vertical density flux is negative (downward) and hence
1064reduces the gravitational PE.
1065
1066For the change in gravitational PE driven by the $k$-flux is
1067\begin{align}
1069  g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho)
1070  &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k
1071    {\mathbb{S}_w}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k
1072    {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\
1073\intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from
1075% and separating out
1076% $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$,
1077% gives two terms. The
1078% first $\rtriad{R}$ term (the only term for $z$-coordinates) is:
1079 &=-\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}
1080\frac{ -\alpha _i^k\delta_{i+ i_p}[T^k]+ \beta_i^k\delta_{i+ i_p}[S^k]} { {e_{1u}}_{\,i + i_p}^{\,k} } \notag \\
1081 &=+\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
1082     \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right) {_i^k\mathbb{R}_{i_p}^{k_p}}
1083\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}},
1084\end{align}
1085using the definition of the triad slope $\rtriad{R}$,
1086\eqref{eq:triad:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+ 1087\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of  $-\alpha_i^k \delta_{k+ 1088 k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$.
1089
1090Where the coordinates slope, the $i$-flux gives a PE change
1091\begin{multline}
1093 g \delta_{i+i_p}[z_T^k]
1094\left[
1095-\alpha _i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (S)
1096\right] \\
1097= +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
1098     \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}
1099\left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)
1100\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}},
1101\end{multline}
1102(using \eqref{eq:triad:skewfluxu}) and so the total PE change
1104\begin{multline}
1106  g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) +
1107g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\
1108= +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
1109     \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)^2
1110\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}}.
1111\end{multline}
1112Where the fluid is stable, with $-\alpha_i^k \delta_{k+ k_p}[T^i]+ 1113\beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative.
1114
1117are masked at the boundaries in exactly the same way as are the triad
1118slopes \rtriad{R} used for the iso-neutral diffusive fluxes, as
1121$\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are
1122masked, and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$
1123and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the
1124$i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$
1125$u$-point is masked. The namelist parameter \np{ln\_botmix\_triad} has
1126no effect on the eddy-induced skew-fluxes.
1127
1128\subsection{ Limiting of the slopes within the interior}\label{sec:triad:limitskew}
1129Presently, the iso-neutral slopes $\tilde{r}_i$ relative
1130to geopotentials are limited to be less than $1/100$, exactly as in
1131calculating the iso-neutral diffusion, \S \ref{sec:triad:limit}. Each
1133
1134\subsection{Tapering within the surface mixed layer}\label{sec:triad:taperskew}
1135The slopes $\tilde{r}_i$ relative to
1136geopotentials (and thus the individual triads \rtriadt{R}) are always tapered linearly from their value immediately below the mixed layer to zero at the
1138option (c) of Fig.~\ref{Fig_eiv_slp}. This linear tapering for the
1139slopes used to calculate the eddy-induced fluxes is
1140unaffected by the value of \np{ln\_triad\_iso}.
1141
1142The justification for this linear slope tapering is that, for $A_e$
1143that is constant or varies only in the horizontal (the most commonly
1144used options in \NEMO: see \S\ref{LDF_coef}), it is
1145equivalent to a horizontal eiv (eddy-induced velocity) that is uniform
1146within the mixed layer \eqref{eq:triad:eiv_v}. This ensures that the
1147eiv velocities do not restratify the mixed layer \citep{Treguier1997,
1148  Danabasoglu_al_2008}. Equivantly, in terms
1149of the skew-flux formulation we use here, the
1150linear slope tapering within the mixed-layer gives a linearly varying
1151vertical flux, and so a tracer convergence uniform in depth (the
1152horizontal flux convergence is relatively insignificant within the mixed-layer).
1153
1155Where the namelist parameter \np{ln\_traldf\_gdia}=true, diagnosed
1156mean eddy-induced velocities are output. Each time step,
1157streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
1158$uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$
1159(integer $i$, integer +1/2 $j$, integer +1/2 $k$) points (see Table
1160\ref{Tab_cell}) respectively. We follow \citep{Griffies_Bk04} and
1161calculate the streamfunction at a given $uw$-point from the
1163\begin{equation}
1165  {\psi_1}_{i+1/2}^{k+1/2}={\quarter}\sum_{\substack{i_p,\,k_p}}
1167\end{equation}
1168The streamfunction $\psi_1$ is calculated similarly at $vw$ points.
1169The eddy-induced velocities are then calculated from the