New URL for NEMO forge!

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
annex_iso.tex in branches/2017/dev_merge_2017/DOC/tex_sub – NEMO

source: branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex @ 9407

Last change on this file since 9407 was 9407, checked in by nicolasmartin, 6 years ago

Complete refactoring of cross-referencing

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