source: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Annex_ISO.tex @ 3267

Last change on this file since 3267 was 3267, checked in by agn, 9 years ago

dev_NEMO_MERGE_2011:(DOC) Finish griffies description, remove need for pstricks

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