source: NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex @ 11571

Last change on this file since 11571 was 11571, checked in by nicolasmartin, 13 months ago

Cosmetic modifications on section titles

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