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

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

Reorganisation for future addition of .rst files from users wiki extraction

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