source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/annex_iso.tex @ 10368

Last change on this file since 10368 was 10368, checked in by smasson, 2 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10365, see #2133

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