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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
apdx_triads.tex in NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

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

Last change on this file since 11561 was 11561, checked in by nicolasmartin, 4 years ago

Apply two thirds ratio on figures width as default

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