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

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

agn: changes to documentation in dev_NEMO_MERGE_2011

—Reordered NEMO_book.tex, changed from 11pt→12pt,

only needs pstricks for logos

—more abbreviations in mathabbrev.sty
—Corrected bug in Abstracts_Foreword.tex that put whole

text into small

—Clarified discussion of LDF operator rotation in
Chap_Model_Basics.tex, Chap_LDF.tex, Annex_B.tex
—Discuss surface tapering of Griffies triads, how they

work in s-coordinates, eddy-induced velocities in
Annex_ISO.tex. Extra ref in Biblio.bib
Extra figures Fig_bdry_triads.pdf, Fig_triad_MLB.pdf
Modified namelist section namtra_ldf

—Minor corrections to Chap_DIA.tex, Chap_CFG.tex,

Chap_SBC.tex, Annex_A.tex, Chap_DOM.tex, Chap_ASM.tex

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