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 trunk/DOC/TexFiles/Chapters – NEMO

source: trunk/DOC/TexFiles/Chapters/Annex_ISO.tex @ 6997

Last change on this file since 6997 was 6997, checked in by nicolasmartin, 7 years ago

Duplication of changes in DOC directory for the trunk

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