source: branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex @ 9393

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

Cleaning of section headings, reinstating the index by mixing \np and \forcode macros, continued conversion of source code inclusions

File size: 58.3 KB
Line 
1\documentclass[../tex_main/NEMO_manual]{subfiles}
2\begin{document}
3% ================================================================
4% Iso-neutral diffusion :
5% ================================================================
6\chapter{Iso-neutral diffusion and eddy advection using triads}
7\label{sec:triad}
8\minitoc
9\pagebreak
10\section{Choice of \protect\ngn{namtra\_ldf} namelist parameters}
11%-----------------------------------------nam_traldf------------------------------------------------------
12\forfile{../namelists/namtra_ldf}
13%---------------------------------------------------------------------------------------------------------
14
15Two scheme are available to perform the iso-neutral diffusion.
16If the namelist logical \np{ln\_traldf\_triad} is set true,
17\NEMO updates both active and passive tracers using the Griffies triad representation
18of iso-neutral diffusion and the eddy-induced advective skew (GM) fluxes.
19If the namelist logical \np{ln\_traldf\_iso} is set true,
20the filtered version of Cox's original scheme (the Standard scheme) is employed (\S\ref{LDF_slp}).
21In the present implementation of the Griffies scheme,
22the advective skew fluxes are implemented even if \np{ln\_traldf\_eiv} is false.
23
24Values of iso-neutral diffusivity and GM coefficient are set as
25described in \S\ref{LDF_coef}. Note that when GM fluxes are used,
26the eddy-advective (GM) velocities are output for diagnostic purposes using xIOS,
27even though the eddy advection is accomplished by means of the skew fluxes.
28
29
30The options specific to the Griffies scheme include:
31\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  `iso-neutral' mixing is accomplished within the surface mixed-layer
34  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}
37  Where \np{ln\_triad\_iso} is set true, the vertical skew flux is further reduced
38  to ensure no vertical buoyancy flux, giving an almost pure
39  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}.
42  If this is set false (the default) then the lateral diffusive fluxes
43  associated with triads partly masked by topography are neglected.
44  If it is set true, however, then these lateral diffusive fluxes are applied,
45  giving smoother bottom tracer fields at the cost of introducing diapycnal mixing.
46\item[\np{rn\_sw\_triad}]  blah blah to be added....
47\end{description}
48The options shared with the Standard scheme include:
49\begin{description}[font=\normalfont]
50\item[\np{ln\_traldf\_msc}]   blah blah to be added
51\item[\np{rn\_slpmax}]  blah blah to be added
52\end{description}
53
54\section{Triad formulation of iso-neutral diffusion}
55\label{sec:triad:iso}
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{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}:
63\begin{subequations} \label{eq:triad: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}
71    \vect{f}^{lT}=-\Alt\Re\cdot\grad T,
72  \end{equation}
73  \begin{equation}
74    \label{eq:triad:PE_iso_tensor:c}
75    \mbox{with}\quad \;\;\Re =
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}
81    \quad \text{and} \quad\grad T=
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}
110  \label{eq:triad:Fijk}
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}
122  \label{eq:triad:i13c}
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}
125  \label{eq:triad:i31c}
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}
132  \label{eq:triad:i33c}
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
143skew-flux, \eqref{eq:triad:i13c}, as
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{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
161gradient, is then \eqref{Eq_tra_ldf_iso}
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
188\subsection{Expression of the skew-flux in terms of triad slopes}
189\citep{Griffies_al_JPO98} introduce a different discretization of the
190off-diagonal terms that nicely solves the problem.
191% Instead of multiplying the mean slope calculated at the $u$-point by
192% the mean vertical gradient at the $u$-point,
193% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
194\begin{figure}[tb] \begin{center}
195    \includegraphics[width=1.05\textwidth]{Fig_GRIFF_triad_fluxes}
196    \caption{ \protect\label{fig:triad:ISO_triad}
197      (a) Arrangement of triads $S_i$ and tracer gradients to
198           give lateral tracer flux from box $i,k$ to $i+1,k$
199      (b) Triads $S'_i$ and tracer gradients to give vertical tracer flux from
200            box $i,k$ to $i,k+1$.}
201 \end{center} \end{figure}
202% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
203They get the skew flux from the products of the vertical gradients at
204each $w$-point surrounding the $u$-point with the corresponding `triad'
205slope calculated from the lateral density gradient across the $u$-point
206divided by the vertical density gradient at the same $w$-point as the
207tracer gradient. See Fig.~\ref{fig:triad:ISO_triad}a, where the thick lines
208denote the tracer gradients, and the thin lines the corresponding
209triads, with slopes $s_1, \dotsc s_4$. The total area-integrated
210skew-flux from tracer cell $i,k$ to $i+1,k$
211\begin{multline}
212  \label{eq:triad:i13}
213  \left( F_u^{13}  \right)_{i+\frac{1}{2}}^k = \Alts_{i+1}^k a_1 s_1
214  \delta _{k+\frac{1}{2}} \left[ T^{i+1}
215  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}  + \Alts _i^k a_2 s_2 \delta
216  _{k+\frac{1}{2}} \left[ T^i
217  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} \\
218   +\Alts _{i+1}^k a_3 s_3 \delta _{k-\frac{1}{2}} \left[ T^{i+1}
219  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}  +\Alts _i^k a_4 s_4 \delta
220  _{k-\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}},
221\end{multline}
222where the contributions of the triad fluxes are weighted by areas
223$a_1, \dotsc a_4$, and $\Alts$ is now defined at the tracer points
224rather than the $u$-points. This discretization gives a much closer
225stencil, and disallows the two-point computational modes.
226
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)
229by multiplying lateral tracer gradients from each of the four
230surrounding $u$-points by the appropriate triad slope:
231\begin{multline}
232  \label{eq:triad:i31}
233  \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} =  \Alts_i^{k+1} a_{1}'
234  s_{1}' \delta _{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1}
235   +\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}\\
236  + \Alts_i^k a_{3}' s_{3}' \delta _{i-\frac{1}{2}} \left[ T^k\right]/{e_{3u}}_{i-\frac{1}{2}}^k
237  +\Alts_i^k a_{4}' s_{4}' \delta _{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k.
238\end{multline}
239
240We notate the triad slopes $s_i$ and $s'_i$ in terms of the `anchor point' $i,k$
241(appearing in both the vertical and lateral gradient), and the $u$- and
242$w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the
243triad as follows (see also Fig.~\ref{fig:triad:ISO_triad}):
244\begin{equation}
245  \label{eq:triad:R}
246  _i^k \mathbb{R}_{i_p}^{k_p}
247  =-\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}}
248  \
249  \frac
250  { \alpha_i^\ \delta_{i+i_p}[T^k] - \beta_i^k \ \delta_{i+i_p}[S^k] }
251  { \alpha_i^\ \delta_{k+k_p}[T^i] - \beta_i^k \ \delta_{k+k_p}[S^i] }.
252\end{equation}
253In calculating the slopes of the local neutral surfaces,
254the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad,
255while the metrics are calculated at the $u$- and $w$-points on the arms.
256
257% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
258\begin{figure}[tb] \begin{center}
259    \includegraphics[width=0.80\textwidth]{Fig_GRIFF_qcells}
260    \caption{   \protect\label{fig:triad:qcells}
261    Triad notation for quarter cells. $T$-cells are inside
262      boxes, while the  $i+\half,k$ $u$-cell is shaded in green and the
263      $i,k+\half$ $w$-cell is shaded in pink.}
264  \end{center} \end{figure}
265% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
266
267Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter
268cell 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.
269Expressing the slopes $s_i$ and $s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation,
270we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$.
271Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$)
272to calculate the lateral flux along its $u$-arm, at $(i+i_p,k)$,
273and then again as an $s'$ to calculate the vertical flux along its $w$-arm at $(i,k+k_p)$.
274Each vertical area $a_i$ used to calculate the lateral flux and horizontal area $a'_i$ used
275to calculate the vertical flux can also be identified as the area across the $u$- and $w$-arms
276of a unique triad, and we notate these areas, similarly to the triad slopes,
277as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$,
278where $e.g.$ in \eqref{eq:triad:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,
279and in \eqref{eq:triad:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.
280
281\subsection{Full triad fluxes}
282A key property of iso-neutral diffusion is that it should not affect
283the (locally referenced) density. In particular there should be no
284lateral or vertical density flux. The lateral density flux disappears so long as the
285area-integrated lateral diffusive flux from tracer cell $i,k$ to
286$i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the
287form
288\begin{equation}
289  \label{eq:triad:i11}
290  \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} =
291  - \left( \Alts_i^{k+1} a_{1} + \Alts_i^{k+1} a_{2} + \Alts_i^k
292    a_{3} + \Alts_i^k a_{4} \right)
293  \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}},
294\end{equation}
295where the areas $a_i$ are as in \eqref{eq:triad:i13}. In this case,
296separating the total lateral flux, the sum of \eqref{eq:triad:i13} and
297\eqref{eq:triad:i11}, into triad components, a lateral tracer
298flux
299\begin{equation}
300  \label{eq:triad:latflux-triad}
301  _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - \Alts_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p}
302  \left(
303    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
304    -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \
305    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
306  \right)
307\end{equation}
308can be identified with each triad. Then, because the
309same metric factors ${e_{3w}}_{\,i}^{\,k+k_p}$ and
310${e_{1u}}_{\,i+i_p}^{\,k}$ are employed for both the density gradients
311in $ _i^k \mathbb{R}_{i_p}^{k_p}$ and the tracer gradients, the lateral
312density flux associated with each triad separately disappears.
313\begin{equation}
314  \label{eq:triad:latflux-rho}
315  {\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
316\end{equation}
317Thus the total flux $\left( F_u^{31} \right) ^i _{i,k+\frac{1}{2}} +
318\left( F_u^{11} \right) ^i _{i,k+\frac{1}{2}}$ from tracer cell $i,k$
319to $i+1,k$ must also vanish since it is a sum of four such triad fluxes.
320
321The squared slope $r_1^2$ in the expression \eqref{eq:triad:i33c} for the
322$_{33}$ component is also expressed in terms of area-weighted
323squared triad slopes, so the area-integrated vertical flux from tracer
324cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is
325\begin{equation}
326  \label{eq:triad:i33}
327  \left( F_w^{33} \right) _i^{k+\frac{1}{2}} =
328    - \left( \Alts_i^{k+1} a_{1}' s_{1}'^2
329    + \Alts_i^{k+1} a_{2}' s_{2}'^2
330    + \Alts_i^k a_{3}' s_{3}'^2
331    + \Alts_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right],
332\end{equation}
333where the areas $a'$ and slopes $s'$ are the same as in
334\eqref{eq:triad:i31}.
335Then, separating the total vertical flux, the sum of \eqref{eq:triad:i31} and
336\eqref{eq:triad:i33}, into triad components,  a vertical flux
337\begin{align}
338  \label{eq:triad:vertflux-triad}
339  _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T)
340  &= \Alts_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p}
341  \left(
342    {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
343    -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \
344    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
345  \right) \\
346  &= - \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\end{align}
349may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$
350associated with a triad then separately disappears (because the
351lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (\rho)$
352disappears). Consequently the total vertical density flux $\left( F_w^{31} \right)_i ^{k+\frac{1}{2}} +
353\left( F_w^{33} \right)_i^{k+\frac{1}{2}}$ from tracer cell $i,k$
354to $i,k+1$ must also vanish since it is a sum of four such triad
355fluxes.
356
357We 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
358the $u$-fluxes and $w$-fluxes in
359\eqref{eq:triad:i31}, \eqref{eq:triad:i13}, \eqref{eq:triad:i11} \eqref{eq:triad:i33} and
360Fig.~\ref{fig:triad:ISO_triad} to  write out the iso-neutral fluxes at $u$- and
361$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
364  \sum_{\substack{i_p,\,k_p}}
365  \begin{pmatrix}
366    {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\
367    \\
368    {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)      \\
369  \end{pmatrix}.
370\end{flalign}
371
372\subsection{Ensuring the scheme does not increase tracer variance}
373\label{sec:triad:variance}
374
375We now require that this operator should not increase the
376globally-integrated tracer variance.
377%This changes according to
378% \begin{align*}
379% &\int_D  D_l^T \; T \;dv \equiv  \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\}    \\
380% &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
381%     \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
382%       + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]  \ T \right\}    \\
383% &\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
384%                 {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T]
385%              + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\
386% \end{align*}
387Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux
388$_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the $u$-point $i+i_p,k$ and
389a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the
390$w$-point $i,k+k_p$.  The lateral flux drives a net rate of change of
391variance, summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
392\begin{multline}
393  {b_T}_{i+i_p-1/2}^k\left(\frac{\partial T}{\partial t}T\right)_{i+i_p-1/2}^k+
394  \quad {b_T}_{i+i_p+1/2}^k\left(\frac{\partial T}{\partial
395      t}T\right)_{i+i_p+1/2}^k \\
396 \begin{aligned}
397  &= -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
398  {\;}_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 \end{aligned}
401\end{multline}
402while the vertical flux similarly drives a net rate of change of
403variance summed over the $T$-points $i,k+k_p-\half$ (above) and
404$i,k+k_p+\half$ (below) of
405\begin{equation}
406\label{eq:triad:dvar_iso_k}
407  _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
408\end{equation}
409The total variance tendency driven by the triad is the sum of these
410two. 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
413\begin{multline*}
414-\Alts_i^k\left \{
415{ } _i^k{\mathbb{A}_u}_{i_p}^{k_p}
416  \left(
417    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
418    - {_i^k\mathbb{R}_{i_p}^{k_p}} \
419    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }\right)\,\delta_{i+ i_p}[T^k] \right.\\
420- \left. { } _i^k{\mathbb{A}_w}_{i_p}^{k_p}
421  \left(
422    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
423    -{\:}_i^k\mathbb{R}_{i_p}^{k_p}
424    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
425  \right) {\,}_i^k\mathbb{R}_{i_p}^{k_p}\delta_{k+ k_p}[T^i]
426\right \}.
427\end{multline*}
428The key point is then that if we require
429$_i^k{\mathbb{A}_u}_{i_p}^{k_p}$ and $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$
430to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by
431\begin{equation}
432  \label{eq:triad:V-A}
433  _i^k\mathbb{V}_{i_p}^{k_p}
434  ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k}
435  ={\;}_i^k{\mathbb{A}_w}_{i_p}^{k_p}\,{e_{3w}}_{\,i}^{\,k + k_p},
436\end{equation}
437the variance tendency reduces to the perfect square
438\begin{equation}
439  \label{eq:triad:perfect-square}
440  -\Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
441  \left(
442    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
443    -{\:}_i^k\mathbb{R}_{i_p}^{k_p}
444    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
445  \right)^2\leq 0.
446\end{equation}
447Thus, the constraint \eqref{eq:triad:V-A} ensures that the fluxes (\ref{eq:triad:latflux-triad}, \ref{eq:triad:vertflux-triad}) associated
448with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase
449the net variance. Since the total fluxes are sums of such fluxes from
450the various triads, this constraint, applied to all triads, is
451sufficient to ensure that the globally integrated variance does not
452increase.
453
454The expression \eqref{eq:triad:V-A} can be interpreted as a discretization
455of the global integral
456\begin{equation}
457  \label{eq:triad:cts-var}
458  \frac{\partial}{\partial t}\int\!\half T^2\, dV =
459  \int\!\mathbf{F}\cdot\nabla T\, dV,
460\end{equation}
461where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the
462lateral and vertical fluxes/unit area
463\[
464\mathbf{F}=\left(
465\left.{}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_u}_{i_p}^{k_p},
466\left.{\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_w}_{i_p}^{k_p}
467 \right)
468\]
469and the gradient
470 \[\nabla T = \left(
471\left.\delta_{i+ i_p}[T^k] \right/ {e_{1u}}_{\,i + i_p}^{\,k},
472\left.\delta_{k+ k_p}[T^i] \right/ {e_{3w}}_{\,i}^{\,k + k_p}
473\right)
474\]
475
476\subsection{Triad volumes in Griffes's scheme and in \NEMO}
477To complete the discretization we now need only specify the triad
478volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. \citet{Griffies_al_JPO98} identify
479these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter
480cells, defined in terms of the distances between $T$, $u$,$f$ and
481$w$-points. This is the natural discretization of
482\eqref{eq:triad:cts-var}. The \NEMO model, however, operates with scale
483factors instead of grid sizes, and scale factors for the quarter
484cells are not defined. Instead, therefore we simply choose
485\begin{equation}
486  \label{eq:triad:V-NEMO}
487  _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k,
488\end{equation}
489as a quarter of the volume of the $u$-cell inside which the triad
490quarter-cell lies. This has the nice property that when the slopes
491$\mathbb{R}$ vanish, the lateral flux from tracer cell $i,k$ to
492$i+1,k$ reduces to the classical form
493\begin{equation}
494  \label{eq:triad:lat-normal}
495-\overline\Alts_{\,i+1/2}^k\;
496\frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}}
497\;\frac{\delta_{i+ 1/2}[T^k] }{{e_{1u}}_{\,i + i_p}^{\,k}}
498 = -\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}}.
499\end{equation}
500In fact if the diffusive coefficient is defined at $u$-points, so that
501we employ $\Alts_{i+i_p}^k$ instead of  $\Alts_i^k$ in the definitions of the
502triad fluxes \eqref{eq:triad:latflux-triad} and \eqref{eq:triad:vertflux-triad},
503we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above.
504
505\subsection{Summary of the scheme}
506The iso-neutral fluxes at $u$- and
507$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}
511    \vect{F}_{\mathrm{iso}}(T) &\equiv
512    \sum_{\substack{i_p,\,k_p}}
513    \begin{pmatrix}
514      {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\
515      \\
516      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)
517    \end{pmatrix},
518  \end{flalign}
519  where \eqref{eq:triad:latflux-triad}:
520  \begin{align}
521    \label{eq:triad:triadfluxu}
522    _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - \Alts_i^k{
523      \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}}
524    \left(
525      \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
526      -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \
527      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
528    \right),\\
529    \intertext{and}
530    _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T)
531    &= \Alts_i^k{\: }\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{3w}}_{\,i}^{\,k+k_p}}
532    \left(
533      {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
534      -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \
535      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
536    \right),\label{eq:triad:triadfluxw}
537  \end{align}
538  with \eqref{eq:triad:V-NEMO}
539  \begin{equation}
540    \label{eq:triad:V-NEMO2}
541    _i^k{\mathbb{V}}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k.
542  \end{equation}
543\end{subequations}
544
545 The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at
546each tracer point:
547\begin{equation} \label{eq:triad:iso_operator} D_l^T = \frac{1}{b_T}
548  \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k
549        {\mathbb{F}_u }_{i_p}^{k_p}} \right] + \delta_{k} \left[
550      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \right\}
551\end{equation}
552where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells.
553The diffusion scheme satisfies the following six properties:
554\begin{description}
555\item[$\bullet$ horizontal diffusion] The discretization of the
556  diffusion operator recovers \eqref{eq:triad:lat-normal} the traditional five-point Laplacian in
557  the limit of flat iso-neutral direction :
558  \begin{equation} \label{eq:triad:iso_property0} D_l^T = \frac{1}{b_T} \
559    \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \;
560      \overline\Alts^{\,i} \; \delta_{i+1/2}[T] \right] \qquad
561    \text{when} \quad { _i^k \mathbb{R}_{i_p}^{k_p} }=0
562  \end{equation}
563
564\item[$\bullet$ implicit treatment in the vertical] Only tracer values
565  associated with a single water column appear in the expression
566  \eqref{eq:triad:i33} for the $_{33}$ fluxes, vertical fluxes driven by
567  vertical gradients. This is of paramount importance since it means
568  that a time-implicit algorithm can be used to solve the vertical
569  diffusion equation. This is necessary
570 since the vertical eddy
571  diffusivity associated with this term,
572  \begin{equation}
573    \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{
574      {\:}_i^k\mathbb{V}_{i_p}^{k_p} \: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
575    \right\}  =
576    \frac{1}{4b_w}\sum_{\substack{i_p, \,k_p}} \left\{
577      {b_u}_{i+i_p}^k\: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
578    \right\},
579 \end{equation}
580  (where $b_w= e_{1w}\,e_{2w}\,e_{3w}$ is the volume of $w$-cells) can be quite large.
581
582\item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of
583  locally referenced potential density is zero. See
584  \eqref{eq:triad:latflux-rho} and \eqref{eq:triad:vertflux-triad2}.
585
586\item[$\bullet$ conservation of tracer] The iso-neutral diffusion
587  conserves tracer content, $i.e.$
588  \begin{equation} \label{eq:triad:iso_property1} \sum_{i,j,k} \left\{ D_l^T \
589      b_T \right\} = 0
590  \end{equation}
591  This property is trivially satisfied since the iso-neutral diffusive
592  operator is written in flux form.
593
594\item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion
595  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      \ b_T \right\} \leq 0
598  \end{equation}
599  The property is demonstrated in
600  \S\ref{sec:triad:variance} above. It is a key property for a diffusion
601  term. It means that it is also a dissipation term, $i.e.$ it
602  dissipates the square of the quantity on which it is applied.  It
603  therefore ensures that, when the diffusivity coefficient is large
604  enough, the field on which it is applied becomes free of grid-point
605  noise.
606
607\item[$\bullet$ self-adjoint operator] The iso-neutral diffusion
608  operator is self-adjoint, $i.e.$
609  \begin{equation} \label{eq:triad:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T
610      \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}
611  \end{equation}
612  In other word, there is no need to develop a specific routine from
613  the adjoint of this operator. We just have to apply the same
614  routine. This property can be demonstrated similarly to the proof of
615  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}
622  - \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
623  \left(
624    \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
625    -{\:}_i^k\mathbb{R}_{i_p}^{k_p}
626    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
627  \right)
628  \left(
629    \frac{ \delta_{i+ i_p}[S^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
630    -{\:}_i^k\mathbb{R}_{i_p}^{k_p}
631    \frac{ \delta_{k+k_p} [S^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
632  \right).
633\end{equation}
634This is symmetrical in $T $ and $S$, so exactly the same term arises
635from the discretization of this triad's contribution towards the
636RHS of \eqref{eq:triad:iso_property3}.
637\end{description}
638
639\subsection{Treatment of the triads at the boundaries}\label{sec:triad:iso_bdry}
640The triad slope can only be defined where both the grid boxes centred at
641the end of the arms exist. Triads that would poke up
642through the upper ocean surface into the atmosphere, or down into the
643ocean floor, must be masked out. See Fig.~\ref{fig:triad:bdry_triads}. Surface layer triads
644$\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and
645$\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be
646specified above the ocean surface are masked (Fig.~\ref{fig:triad:bdry_triads}a): this ensures that lateral
647tracer gradients produce no flux through the ocean surface. However,
648to prevent surface noise, it is customary to retain the $_{11}$ contributions towards
649the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and
650$\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer
651fluxes. Similar comments apply to triads that would intersect the
652ocean floor (Fig.~\ref{fig:triad:bdry_triads}b). Note that both near bottom
653triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and
654$\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$
655or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is
656masked. The associated lateral fluxes (grey-black dashed line) are
657masked if \np{ln\_botmix\_triad}\forcode{ = .false.}, but left unmasked,
658giving bottom mixing, if \np{ln\_botmix\_triad}\forcode{ = .true.}.
659
660The default option \np{ln\_botmix\_triad}\forcode{ = .false.} is suitable when the
661bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}\forcode{ = 1}),
662or  for simple idealized  problems. For setups with topography without
663bbl mixing, \np{ln\_botmix\_triad}\forcode{ = .true.} may be necessary.
664% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
665\begin{figure}[h] \begin{center}
666    \includegraphics[width=0.60\textwidth]{Fig_GRIFF_bdry_triads}
667    \caption{  \protect\label{fig:triad:bdry_triads}
668      (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer
669      points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad
670      slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$
671      (blue) poking through the ocean surface are masked (faded in
672      figure). However, the lateral $_{11}$ contributions towards
673      $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$
674      (yellow line) are still applied, giving diapycnal diffusive
675      fluxes.\newline
676      (b) Both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and
677      $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$
678      or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point
679      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.}}
682 \end{center} \end{figure}
683% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
684
685\subsection{ Limiting of the slopes within the interior}\label{sec:triad:limit}
686As discussed in \S\ref{LDF_slp_iso}, iso-neutral slopes relative to
687geopotentials must be bounded everywhere, both for consistency with the small-slope
688approximation and for numerical stability \citep{Cox1987,
689  Griffies_Bk04}. The bound chosen in \NEMO is applied to each
690component of the slope separately and has a value of $1/100$ in the ocean interior.
691%, ramping linearly down above 70~m depth to zero at the surface
692It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to
693geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to
694geopotentials) \eqref{Eq_PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate
695surfaces, so we require
696\begin{equation*}
697  |\tilde{r}_i|\leq \tilde{r}_\mathrm{max}=0.01.
698\end{equation*}
699and then recalculate the slopes $r_i$ relative to coordinates.
700Each individual triad slope
701 \begin{equation}
702   \label{eq:triad:Rtilde}
703_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}}
704 \end{equation}
705is limited like this and then the corresponding
706$_i^k\mathbb{R}_{i_p}^{k_p} $ are recalculated and combined to form the fluxes.
707Note that where the slopes have been limited, there is now a non-zero
708iso-neutral density flux that drives dianeutral mixing.  In particular this iso-neutral density flux
709is always downwards, and so acts to reduce gravitational potential energy.
710
711\subsection{Tapering within the surface mixed layer}\label{sec:triad:taper}
712Additional tapering of the iso-neutral fluxes is necessary within the
713surface mixed layer. When the Griffies triads are used, we offer two
714options for this.
715
716\subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper}
717This is the option activated by the default choice
718\np{ln\_triad\_iso}\forcode{ = .false.}. Slopes $\tilde{r}_i$ relative to
719geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the
720surface, as described in option (c) of Fig.~\ref{Fig_eiv_slp}, to values
721\begin{subequations}
722  \begin{equation}
723   \label{eq:triad:rmtilde}
724     \rMLt =
725  -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h,
726  \end{equation}
727and then the $r_i$ relative to vertical coordinate surfaces are appropriately
728adjusted to
729  \begin{equation}
730   \label{eq:triad:rm}
731 \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h.
732  \end{equation}
733\end{subequations}
734Thus the diffusion operator within the mixed layer is given by:
735\begin{equation} \label{eq:triad:iso_tensor_ML}
736D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
737\mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c}
738 1 \hfill & 0 \hfill & {-\rML[1]}\hfill \\
739 0 \hfill & 1 \hfill & {-\rML[2]} \hfill \\
740 {-\rML[1]}\hfill &   {-\rML[2]} \hfill & {\rML[1]^2+\rML[2]^2} \hfill
741\end{array} }} \right)
742\end{equation}
743
744This slope tapering gives a natural connection between tracer in the
745mixed-layer and in isopycnal layers immediately below, in the
746thermocline. It is consistent with the way the $\tilde{r}_i$ are
747tapered within the mixed layer (see \S\ref{sec:triad:taperskew} below)
748so as to ensure a uniform GM eddy-induced velocity throughout the
749mixed layer. However, it gives a downwards density flux and so acts so
750as to reduce potential energy in the same way as does the slope
751limiting discussed above in \S\ref{sec:triad:limit}.
752 
753As in \S\ref{sec:triad:limit} above, the tapering
754\eqref{eq:triad:rmtilde} is applied separately to each triad
755$_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the
756$_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume
757$z$-coordinates in the following; the conversion from
758$\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described
759above by \eqref{eq:triad:Rtilde}.
760\begin{enumerate}
761\item Mixed-layer depth is defined so as to avoid including regions of weak
762vertical stratification in the slope definition.
763 At each $i,j$ (simplified to $i$ in
764Fig.~\ref{fig:triad:MLB_triad}), we define the mixed-layer by setting
765the vertical index of the tracer point immediately below the mixed
766layer, $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point)
767such that the potential density
768${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, where $i,k_{10}$ is
769the tracer gridbox within which the depth reaches 10~m. See the left
770side of Fig.~\ref{fig:triad:MLB_triad}. We use the $k_{10}$-gridbox
771instead of the surface gridbox to avoid problems e.g.\ with thin
772daytime mixed-layers. Currently we use the same
773$\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is
774used to output the diagnosed mixed-layer depth
775$h_{\mathrm{ML}}=|z_{W}|_{k_{\mathrm{ML}}+1/2}$, the depth of the $w$-point
776above the $i,k_{\mathrm{ML}}$ tracer point.
777
778\item We define `basal' triad slopes
779${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ as the slopes
780of those triads whose vertical `arms' go down from the
781$i,k_{\mathrm{ML}}$ tracer point to the $i,k_{\mathrm{ML}}-1$ tracer point
782below. This is to ensure that the vertical density gradients
783associated with these basal triad slopes
784${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are
785representative of the thermocline. The four basal triads defined in the bottom part
786of Fig.~\ref{fig:triad:MLB_triad} are then
787\begin{align}
788  {\:}_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\\
791\intertext{with e.g.\ the green triad}
792{\:}_i{\mathbb{R}_{\mathrm{base}}}_{1/2}^{-1/2}&=
793{\:}^{k_{\mathrm{ML}}}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2}. \notag
794\end{align}
795The vertical flux associated with each of these triads passes through the $w$-point
796$i,k_{\mathrm{ML}}-1/2$ lying \emph{below} the $i,k_{\mathrm{ML}}$ tracer point,
797so it is this depth
798\begin{equation}
799  \label{eq:triad:zbase}
800  {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2}
801\end{equation}
802(one gridbox deeper than the
803diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper
804the slopes in \eqref{eq:triad:rmtilde}.
805\item Finally, we calculate the adjusted triads
806${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within the mixed
807layer, by multiplying the appropriate
808${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ by the ratio of
809the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_{\mathrm{base}}}_{\,i}$. For
810instance the green triad centred on $i,k$
811\begin{align}
812  {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,1/2}^{-1/2} &=
813\frac{{z_w}_{k-1/2}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2}
814\notag \\
815\intertext{and more generally}
816 {\:}_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\end{align}
819\end{enumerate}
820
821% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
822\begin{figure}[h]
823%  \fcapside {
824    \caption{\protect\label{fig:triad:MLB_triad} Definition of
825      mixed-layer depth and calculation of linearly tapered
826      triads. The figure shows a water column at a given $i,j$
827      (simplified to $i$), with the ocean surface at the top. Tracer points are denoted by
828      bullets, and black lines the edges of the tracer cells; $k$
829      increases upwards. \newline
830      \hspace{5 em}We define the mixed-layer by setting the vertical index
831      of the tracer point immediately below the mixed layer,
832      $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point)
833      such that ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$,
834      where $i,k_{10}$ is the tracer gridbox within which the depth
835      reaches 10~m. We calculate the triad slopes within the mixed
836      layer by linearly tapering them from zero (at the surface) to
837      the `basal' slopes, the slopes of the four triads passing through the
838      $w$-point $i,k_{\mathrm{ML}}-1/2$ (blue square),
839      ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$. Triads with
840    different $i_p,k_p$, denoted by different colours, (e.g. the green
841    triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.}
842%}
843  {\includegraphics[width=0.60\textwidth]{Fig_GRIFF_MLB_triads}}
844\end{figure}
845% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
846
847\subsubsection{Additional truncation of skew iso-neutral flux components}
848\label{sec:triad:Gerdes-taper}
849The alternative option is activated by setting \np{ln\_triad\_iso} =
850  true. This retains the same tapered slope $\rML$  described above for the
851calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the
852vertical tracer flux driven by vertical tracer gradients), but
853replaces the $\rML$ in the skew term by
854\begin{equation}
855  \label{eq:triad:rm*}
856  \rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i,
857\end{equation}
858giving a ML diffusive operator
859\begin{equation} \label{eq:triad:iso_tensor_ML2}
860D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
861\mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c}
862 1 \hfill & 0 \hfill & {-\rML[1]^*}\hfill \\
863 0 \hfill & 1 \hfill & {-\rML[2]^*} \hfill \\
864 {-\rML[1]^*}\hfill &   {-\rML[2]^*} \hfill & {\rML[1]^2+\rML[2]^2} \hfill \\
865\end{array} }} \right).
866\end{equation}
867This operator
868\footnote{To ensure good behaviour where horizontal density
869  gradients are weak, we in fact follow \citet{Gerdes1991} and set
870$\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.}
871then has the property it gives no vertical density flux, and so does
872not change the potential energy.
873This approach is similar to multiplying the iso-neutral  diffusion
874coefficient by $\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep
875slopes, as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}).
876Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$
877
878In practice, this approach gives weak vertical tracer fluxes through
879the mixed-layer, as well as vanishing density fluxes. While it is
880theoretically advantageous that it does not change the potential
881energy, it may give a discontinuity between the
882fluxes within the mixed-layer (purely horizontal) and just below (along
883iso-neutral surfaces).
884% This may give strange looking results,
885% particularly where the mixed-layer depth varies strongly laterally.
886% ================================================================
887% Skew flux formulation for Eddy Induced Velocity :
888% ================================================================
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}
892
893 When Gent and McWilliams's [1990] diffusion is used,
894an additional advection term is added. The associated velocity is the so called
895eddy induced velocity, the formulation of which depends on the slopes of iso-
896neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used
897here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo}
898is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo}
899+ \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates.
900
901The eddy induced velocity is given by:
902\begin{subequations} \label{eq:triad:eiv}
903\begin{equation}\label{eq:triad:eiv_v}
904\begin{split}
905 u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\
906 v^* & = - \frac{1}{e_{3}}\;          \partial_j\psi_2,    \\
907w^* & =    \frac{1}{e_{1}e_{2}}\; \left\{ \partial_\left( e_{2} \, \psi_1\right)
908                         + \partial_\left( e_{1} \, \psi_2\right) \right\},
909\end{split}
910\end{equation}
911where the streamfunctions $\psi_i$ are given by
912\begin{equation} \label{eq:triad:eiv_psi}
913\begin{split}
914\psi_1 & = A_{e} \; \tilde{r}_1,   \\
915\psi_2 & = A_{e} \; \tilde{r}_2,
916\end{split}
917\end{equation}
918\end{subequations}
919with $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.
920
921The traditional way to implement this additional advection is to add
922it to the Eulerian velocity prior to computing the tracer
923advection. This is implemented if \key{traldf\_eiv} is set in the
924default implementation, where \np{ln\_traldf\_triad} is set
925false. This allows us to take advantage of all the advection schemes
926offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$
927order advection scheme. This is particularly useful for passive
928tracers where \emph{positivity} of the advection scheme is of
929paramount importance.
930
931However, when \np{ln\_traldf\_triad} is set true, \NEMO instead
932implements eddy induced advection according to the so-called skew form
933\citep{Griffies_JPO98}. It is based on a transformation of the advective fluxes
934using the non-divergent nature of the eddy induced velocity.
935For example in the (\textbf{i},\textbf{k}) plane, the tracer advective
936fluxes per unit area in $ijk$ space can be
937transformed as follows:
938\begin{flalign*}
939\begin{split}
940\textbf{F}_{\mathrm{eiv}}^T =
941\begin{pmatrix}
942           {e_{2}\,e_{3}\;  u^*}       \\
943      {e_{1}\,e_{2}\; w^*}  \\
944\end{pmatrix}   \;   T
945&=
946\begin{pmatrix}
947           { - \partial_k \left( e_{2} \,\psi_1 \right) \; T \;}     \\
948      {+ \partial_\left( e_{2} \, \psi_1 \right) \; T \;}    \\
949\end{pmatrix}        \\
950&=
951\begin{pmatrix}
952           { - \partial_k \left( e_{2} \, \psi_\; T \right) \;}  \\
953      {+ \partial_\left( e_{2} \,\psi_1 \; T \right) \;}  \\
954\end{pmatrix}
955 +
956\begin{pmatrix}
957           {+ e_{2} \, \psi_\; \partial_k T}  \\
958      { - e_{2} \, \psi_\; \partial_i  T}  \\
959\end{pmatrix}
960\end{split}
961\end{flalign*}
962and since the eddy induced velocity field is non-divergent, we end up with the skew
963form of the eddy induced advective fluxes per unit area in $ijk$ space:
964\begin{equation} \label{eq:triad:eiv_skew_ijk}
965\textbf{F}_\mathrm{eiv}^T = \begin{pmatrix}
966           {+ e_{2} \, \psi_\; \partial_k T}   \\
967      { - e_{2} \, \psi_\; \partial_i  T}  \\
968                                 \end{pmatrix}
969\end{equation}
970The total fluxes per unit physical area are then
971\begin{equation}\label{eq:triad:eiv_skew_physical}
972\begin{split}
973 f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T   \\
974 f^*_2 & = \frac{1}{e_{3}}\; \psi_2 \partial_k T   \\
975 f^*_3 & =  -\frac{1}{e_{1}e_{2}}\; \left\{ e_{2} \psi_1 \partial_i T
976   + e_{1} \psi_2 \partial_j T \right\}. \\
977\end{split}
978\end{equation}
979Note that Eq.~ \eqref{eq:triad:eiv_skew_physical} takes the same form whatever the
980vertical 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.
982The tendency associated with eddy induced velocity is then simply the convergence
983of the fluxes (\ref{eq:triad:eiv_skew_ijk}, \ref{eq:triad:eiv_skew_physical}), so
984\begin{equation} \label{eq:triad:skew_eiv_conv}
985\frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[
986  \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right)
987  + \frac{\partial}{\partial j} \left( e_\;
988    \psi_2 \partial_k T\right)
989 -  \frac{\partial}{\partial k} \left( e_{2} \psi_1 \partial_i T
990   + e_{1} \psi_2 \partial_j T \right\right]
991\end{equation}
992 It naturally conserves the tracer content, as it is expressed in flux
993 form. Since it has the same divergence as the advective form it also
994 preserves the tracer variance.
995
996\subsection{Discrete skew flux formulation}
997The 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
999expressed in terms of the triad slopes, as in Fig.~\ref{fig:triad:ISO_triad}
1000and Eqs~(\ref{eq:triad:i13}, \ref{eq:triad:i31}); but now in terms of the triad slopes
1001$\tilde{\mathbb{R}}$ relative to geopotentials instead of the
1002$\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
1004defining $A_e$ at $T$-points is then given by:
1005
1006
1007\begin{subequations}\label{eq:triad:allskewflux}
1008  \begin{flalign}\label{eq:triad:vect_skew_flux}
1009    \vect{F}_{\mathrm{eiv}}(T) &\equiv
1010    \sum_{\substack{i_p,\,k_p}}
1011    \begin{pmatrix}
1012      {_{i+1/2-i_p}^k {\mathbb{S}_u}_{i_p}^{k_p} } (T)      \\
1013      \\
1014      {_i^{k+1/2-k_p} {\mathbb{S}_w}_{i_p}^{k_p} } (T)      \\
1015    \end{pmatrix},
1016  \end{flalign}
1017  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  \begin{align}
1020    \label{eq:triad:skewfluxu}
1021    _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \quarter {A_e}_i^k{
1022      \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}}
1023     \ {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} \
1024      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} },
1025   \\
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}:}
1028    _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T)
1029    &= -\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  \end{align}
1032\end{subequations}
1033
1034Such a discretisation is consistent with the iso-neutral
1035operator as it uses the same definition for the slopes.  It also
1036ensures the following two key properties.
1037
1038\subsubsection{No change in tracer variance}
1039The discretization conserves tracer variance, $i.e.$ it does not
1040include a diffusive component but is a `pure' advection term. This can
1041be seen
1042%either from Appendix \ref{Apdx_eiv_skew} or
1043by considering the
1044fluxes associated with a given triad slope
1045$_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
1047associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$
1048drives a net rate of change of variance, summed over the two
1049$T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
1050\begin{equation}
1051\label{eq:triad:dvar_eiv_i}
1052  _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k],
1053\end{equation}
1054while the associated vertical skew-flux gives a variance change summed over the
1055$T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of
1056\begin{equation}
1057\label{eq:triad:dvar_eiv_k}
1058  _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
1059\end{equation}
1060Inspection of the definitions (\ref{eq:triad:skewfluxu}, \ref{eq:triad:skewfluxw})
1061shows that these two variance changes (\ref{eq:triad:dvar_eiv_i}, \ref{eq:triad:dvar_eiv_k})
1062sum to zero. Hence the two fluxes associated with each triad make no
1063net contribution to the variance budget.
1064
1065\subsubsection{Reduction in gravitational PE}
1066The vertical density flux associated with the vertical skew-flux
1067always has the same sign as the vertical density gradient; thus, so
1068long as the fluid is stable (the vertical density gradient is
1069negative) the vertical density flux is negative (downward) and hence
1070reduces the gravitational PE.
1071
1072For the change in gravitational PE driven by the $k$-flux is
1073\begin{align}
1074  \label{eq:triad:vert_densityPE}
1075  g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho)
1076  &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k
1077    {\mathbb{S}_w}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k
1078    {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\
1079\intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from
1080  \eqref{eq:triad:skewfluxw}, gives}
1081% and separating out
1082% $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$,
1083% gives two terms. The
1084% first $\rtriad{R}$ term (the only term for $z$-coordinates) is:
1085 &=-\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}
1086\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 \\
1087 &=+\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
1088     \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}}
1089\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}},
1090\end{align}
1091using 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\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of  $-\alpha_i^k \delta_{k+
1094  k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$.
1095
1096Where the coordinates slope, the $i$-flux gives a PE change
1097\begin{multline}
1098  \label{eq:triad:lat_densityPE}
1099 g \delta_{i+i_p}[z_T^k]
1100\left[
1101-\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)
1102\right] \\
1103= +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
1104     \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}
1105\left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)
1106\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}},
1107\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
1110\begin{multline}
1111  \label{eq:triad:tot_densityPE}
1112  g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) +
1113g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\
1114= +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
1115     \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
1116\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}}.
1117\end{multline}
1118Where the fluid is stable, with $-\alpha_i^k \delta_{k+ k_p}[T^i]+
1119\beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative.
1120
1121\subsection{Treatment of the triads at the boundaries}\label{sec:triad:skew_bdry}
1122Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes
1123are masked at the boundaries in exactly the same way as are the triad
1124slopes \rtriad{R} used for the iso-neutral diffusive fluxes, as
1125described in \S\ref{sec:triad:iso_bdry} and
1126Fig.~\ref{fig:triad:bdry_triads}. Thus surface layer triads
1127$\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are
1128masked, and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$
1129and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the
1130$i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$
1131$u$-point is masked. The namelist parameter \np{ln\_botmix\_triad} has
1132no effect on the eddy-induced skew-fluxes.
1133
1134\subsection{Limiting of the slopes within the interior}\label{sec:triad:limitskew}
1135Presently, the iso-neutral slopes $\tilde{r}_i$ relative
1136to geopotentials are limited to be less than $1/100$, exactly as in
1137calculating the iso-neutral diffusion, \S \ref{sec:triad:limit}. Each
1138individual triad \rtriadt{R} is so limited.
1139
1140\subsection{Tapering within the surface mixed layer}\label{sec:triad:taperskew}
1141The slopes $\tilde{r}_i$ relative to
1142geopotentials (and thus the individual triads \rtriadt{R}) are always tapered linearly from their value immediately below the mixed layer to zero at the
1143surface \eqref{eq:triad:rmtilde}, as described in \S\ref{sec:triad:lintaper}. This is
1144option (c) of Fig.~\ref{Fig_eiv_slp}. This linear tapering for the
1145slopes used to calculate the eddy-induced fluxes is
1146unaffected by the value of \np{ln\_triad\_iso}.
1147
1148The justification for this linear slope tapering is that, for $A_e$
1149that is constant or varies only in the horizontal (the most commonly
1150used options in \NEMO: see \S\ref{LDF_coef}), it is
1151equivalent to a horizontal eiv (eddy-induced velocity) that is uniform
1152within the mixed layer \eqref{eq:triad:eiv_v}. This ensures that the
1153eiv velocities do not restratify the mixed layer \citep{Treguier1997,
1154  Danabasoglu_al_2008}. Equivantly, in terms
1155of the skew-flux formulation we use here, the
1156linear slope tapering within the mixed-layer gives a linearly varying
1157vertical flux, and so a tracer convergence uniform in depth (the
1158horizontal flux convergence is relatively insignificant within the mixed-layer).
1159
1160\subsection{Streamfunction diagnostics}\label{sec:triad:sfdiag}
1161Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.}, diagnosed
1162mean eddy-induced velocities are output. Each time step,
1163streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
1164$uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$
1165(integer $i$, integer +1/2 $j$, integer +1/2 $k$) points (see Table
1166\ref{Tab_cell}) respectively. We follow \citep{Griffies_Bk04} and
1167calculate the streamfunction at a given $uw$-point from the
1168surrounding four triads according to:
1169\begin{equation}
1170  \label{eq:triad:sfdiagi}
1171  {\psi_1}_{i+1/2}^{k+1/2}={\quarter}\sum_{\substack{i_p,\,k_p}}
1172  {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}.
1173\end{equation}
1174The streamfunction $\psi_1$ is calculated similarly at $vw$ points.
1175The eddy-induced velocities are then calculated from the
1176straightforward discretisation of \eqref{eq:triad:eiv_v}:
1177\begin{equation}\label{eq:triad:eiv_v_discrete}
1178\begin{split}
1179 {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),   \\
1180 {v^*}_{j+1/2}^{k} & = - \frac{1}{{e_{3v}}_{j}^{k}}\left({\psi_2}_{j+1/2}^{k+1/2}-{\psi_2}_{j+1/2}^{k+1/2}\right),   \\
1181 {w^*}_{i,j}^{k+1/2} & =    \frac{1}{e_{1t}e_{2t}}\; \left\{
1182 {e_{2u}}_{i+1/2}^{k+1/2} \,{\psi_1}_{i+1/2}^{k+1/2} -
1183 {e_{2u}}_{i-1/2}^{k+1/2} \,{\psi_1}_{i-1/2}^{k+1/2} \right. + \\
1184\phantom{=} & \qquad\qquad\left. {e_{2v}}_{j+1/2}^{k+1/2} \,{\psi_2}_{j+1/2}^{k+1/2} - {e_{2v}}_{j-1/2}^{k+1/2} \,{\psi_2}_{j-1/2}^{k+1/2} \right\},
1185\end{split}
1186\end{equation}
1187\end{document}
Note: See TracBrowser for help on using the repository browser.