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.
annex_iso.tex in NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

source: NEMO/trunk/doc/latex/NEMO/subfiles/annex_iso.tex @ 11435

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

Various corrections on chapters

Cleaning the indexes by fixing/removing wrong entries (or appending a ? to unknown items) and
improve the classification with new index definitions for CPP keys and namelist blocks:

  • from \key{...} cmd, key_ prefix no longer precedes the index entry
  • namelist block declaration moves from \ngn{nam...} to \nam{...} (i.e. \ngn{namtra\_ldf} -> \nam{tra\_ldf}) The expected prefix nam is added to the printed word but not the index entry.

Now we have indexes with a better sorting instead of all CPP keys under 'K' and namelists blocks under 'N'.

Fix missing space issues with alias commands by adding a trailing backslash (\NEMO\, \ie\, \eg\, ...).
There is no perfect solution for this, and I prefer not using a particular package to solve it.

Review the initial LaTeX code snippet for the historic changes in chapters

Finally, for readability and future diff visualisations, please avoid writing paragraphs with continuous lines.
Break the lines around 80 to 100 characters long

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