source: NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex @ 11693

Last change on this file since 11693 was 11693, checked in by nicolasmartin, 13 months ago

Macros renaming

File size: 40.2 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4
5\chapter{Note on some algorithms}
6\label{apdx:ALGOS}
7
8\thispagestyle{plain}
9
10\chaptertoc
11
12\paragraph{Changes record} ~\\
13
14{\footnotesize
15  \begin{tabularx}{\textwidth}{l||X|X}
16    Release & Author(s) & Modifications \\
17    \hline
18    {\em   4.0} & {\em ...} & {\em ...} \\
19    {\em   3.6} & {\em ...} & {\em ...} \\
20    {\em   3.4} & {\em ...} & {\em ...} \\
21    {\em <=3.4} & {\em ...} & {\em ...}
22  \end{tabularx}
23}
24
25\clearpage
26
27This appendix some on going consideration on algorithms used or planned to be used in \NEMO.
28
29%% =================================================================================================
30\section{Upstream Biased Scheme (UBS) (\protect\np[=.true.]{ln_traadv_ubs}{ln\_traadv\_ubs})}
31\label{sec:ALGOS_tra_adv_ubs}
32
33The UBS advection scheme is an upstream biased third order scheme based on
34an upstream-biased parabolic interpolation.
35It is also known as Cell Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective Kinematics).
36For example, in the $i$-direction:
37\begin{equation}
38  \label{eq:ALGOS_tra_adv_ubs2}
39  \tau_u^{ubs} = \left\{
40    \begin{aligned}
41      & \tau_u^{cen4} + \frac{1}{12} \,\tau"_i     & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
42      & \tau_u^{cen4} - \frac{1}{12} \,\tau"_{i+1} & \quad \text{if }\ u_{i+1/2}       <       0
43    \end{aligned}
44  \right.
45\end{equation}
46or equivalently, the advective flux is
47\begin{equation}
48  \label{eq:ALGOS_tra_adv_ubs2}
49  U_{i+1/2} \ \tau_u^{ubs}
50  =U_{i+1/2} \ \overline{ T_i - \frac{1}{6}\,\tau"_i }^{\,i+1/2}
51  - \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
52\end{equation}
53where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and
54$\tau "_i =\delta_i \left[ {\delta_{i+1/2} \left[ \tau \right]} \right]$.
55By choosing this expression for $\tau "$ we consider a fourth order approximation of $\partial_i^2$ with
56a constant i-grid spacing ($\Delta i=1$).
57
58Alternative choice: introduce the scale factors:
59$\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta_i \left[ \frac{e_{2u} e_{3u} }{e_{1u} }\delta_{i+1/2}[\tau] \right]$.
60
61This results in a dissipatively dominant (\ie\ hyper-diffusive) truncation error
62\citep{shchepetkin.mcwilliams_OM05}.
63The overall performance of the advection scheme is similar to that reported in \cite{farrow.stevens_JPO95}.
64It is a relatively good compromise between accuracy and smoothness.
65It is not a \emph{positive} scheme meaning false extrema are permitted but
66the amplitude of such are significantly reduced over the centred second order method.
67Nevertheless it is not recommended to apply it to a passive tracer that requires positivity.
68
69The intrinsic diffusion of UBS makes its use risky in the vertical direction where
70the control of artificial diapycnal fluxes is of paramount importance.
71It has therefore been preferred to evaluate the vertical flux using the TVD scheme when
72\np[=.true.]{ln_traadv_ubs}{ln\_traadv\_ubs}.
73
74For stability reasons, in \autoref{eq:TRA_adv_ubs}, the first term which corresponds to
75a second order centred scheme is evaluated using the \textit{now} velocity (centred in time) while
76the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity
77(forward in time).
78This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme.
79UBS and QUICK schemes only differ by one coefficient.
80Substituting 1/6 with 1/8 in (\autoref{eq:TRA_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
81This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded.
82Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme.
83
84NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can be controlled by vertical advection
85(not vertical diffusion which is usually solved using an implicit scheme).
86Computer time can be saved by using a time-splitting technique on vertical advection.
87This possibility have been implemented and validated in ORCA05-L301.
88It is not currently offered in the current reference version.
89
90NB 2: In a forthcoming release four options will be proposed for the vertical component used in the UBS scheme.
91$\tau_w^{ubs}$ will be evaluated using either \textit{(a)} a centered $2^{nd}$ order scheme,
92or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative parabolic splines following
93\citet{shchepetkin.mcwilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS.
94The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme.
95
96NB 3: It is straight forward to rewrite \autoref{eq:TRA_adv_ubs} as follows:
97\begin{equation}
98  \label{eq:ALGOS_tra_adv_ubs2}
99  \tau_u^{ubs} = \left\{
100    \begin{aligned}
101      & \tau_u^{cen4} + \frac{1}{12} \tau"_i    & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
102      & \tau_u^{cen4} - \frac{1}{12} \tau"_{i+1}   & \quad \text{if }\ u_{i+1/2}       <       0
103    \end{aligned}
104  \right.
105\end{equation}
106or equivalently
107\begin{equation}
108  \label{eq:ALGOS_tra_adv_ubs2}
109  \begin{split}
110    e_{2u} e_{3u}\,u_{i+1/2} \ \tau_u^{ubs}
111    &= e_{2u} e_{3u}\,u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\tau"_i }^{\,i+1/2} \\
112    & - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
113  \end{split}
114\end{equation}
115\autoref{eq:TRA_adv_ubs2} has several advantages.
116First it clearly evidences that the UBS scheme is based on the fourth order scheme to which
117is added an upstream biased diffusive term.
118Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step,
119not only the $2^{th}$ order part as stated above using \autoref{eq:TRA_adv_ubs}.
120Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient which
121is simply proportional to the velocity.
122
123laplacian diffusion:
124\begin{equation}
125  \label{eq:ALGOS_tra_ldf_lap}
126  \begin{split}
127    D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\;  e_{3T} } &\left[ {\quad \delta_i
128        \left[ {A_u^{lT} \frac{e_{2u} e_{3u} }{e_{1u} }\;\delta_{i+1/2}
129            \left[ T \right]} \right]} \right. \\
130    &\ \left. {+\; \delta_j \left[
131          {A_v^{lT} \left( {\frac{e_{1v} e_{3v} }{e_{2v} }\;\delta_{j+1/2} \left[ T
132                \right]} \right)} \right]\quad } \right]
133  \end{split}
134\end{equation}
135
136bilaplacian:
137\begin{equation}
138  \label{eq:ALGOS_tra_ldf_lap}
139  \begin{split}
140    D_T^{lT} =&-\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
141    & \delta_i \left\sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta_{i+1/2}
142      \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}}
143        \delta_i \left[ \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta_{i+1/2}
144          [T] \right] \right] \right]
145  \end{split}
146\end{equation}
147with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$,
148\ie\ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$
149it comes:
150\begin{equation}
151  \label{eq:ALGOS_tra_ldf_lap}
152  \begin{split}
153    D_T^{lT} =&-\frac{1}{12}\,\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
154    & \delta_i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta_{i+1/2}
155      \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}}
156        \delta_i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta_{i+1/2}
157          [T] \right] \right] \right]
158  \end{split}
159\end{equation}
160if the velocity is uniform (\ie\ $|u|=cst$) then the diffusive flux is
161\begin{equation}
162  \label{eq:ALGOS_tra_ldf_lap}
163  \begin{split}
164    F_u^{lT} = - \frac{1}{12}
165    e_{2u}\,e_{3u}\,|u| \;\sqrt{ e_{1u}}\,\delta_{i+1/2}
166    \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}}
167      \delta_i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}}\:\delta_{i+1/2}
168        [T] \right] \right]
169  \end{split}
170\end{equation}
171beurk....  reverte the logic: starting from the diffusive part of the advective flux it comes:
172
173\begin{equation}
174  \label{eq:ALGOS_tra_adv_ubs2}
175  \begin{split}
176    F_u^{lT} &= - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
177  \end{split}
178\end{equation}
179if the velocity is uniform (\ie\ $|u|=cst$) and
180choosing $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta_i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta_{i+1/2}[\tau] \right]$
181
182sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$):
183\begin{equation}
184  \label{eq:ALGOS_tra_adv_ubs2}
185  \begin{split}
186    F_u^{lT} &= - \frac{1}{12} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{e_{1T}^3\,|u|}{e_{1T}e_{2T}\,e_{3T}}\,\delta_i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta_{i+1/2}[\tau] \right] \right]
187  \end{split}
188\end{equation}
189which leads to ${A_T^{lT}}^2 = \frac{1}{12} {e_{1T}}^3\ \overline{|u|}^{\,i+1/2}$
190
191sol 2 coefficient at u-point: split $|u|$ into $\sqrt{|u|}$ and $e_{1T}$ into $\sqrt{e_{1u}}$
192\begin{equation}
193  \label{eq:ALGOS_tra_adv_ubs2}
194  \begin{split}
195    F_u^{lT} &= - \frac{1}{12} {e_{1u}}^1 \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{1}{e_{2T}\,e_{3T}}\,\delta_i \left[ \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u} }{e_{1u} } \delta_{i+1/2}[\tau] \right] \right] \\
196    &= - \frac{1}{12} e_{1u} \sqrt{e_{1u}|u|\,} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{1}{e_{1T}\,e_{2T}\,e_{3T}}\,\delta_i \left[ e_{1u} \sqrt{e_{1u}|u|\,} \frac{e_{2u} e_{3u} }{e_{1u}} \delta_{i+1/2}[\tau] \right] \right]
197  \end{split}
198\end{equation}
199which leads to ${A_u^{lT}} = \frac{1}{12} {e_{1u}}^3\ |u|$
200
201%% =================================================================================================
202\section{Leapfrog energetic}
203\label{sec:ALGOS_LF}
204
205We adopt the following semi-discrete notation for time derivative.
206Given the values of a variable $q$ at successive time step,
207the time derivation and averaging operators at the mid time step are:
208\[
209  % \label{eq:ALGOS_dt_mt}
210  \begin{split}
211    \delta_{t+\rdt/2} [q]     &=  \  \ \,   q^{t+\rdt- q^{t}      \\
212    \overline q^{\,t+\rdt/2} &= \left\{ q^{t+\rdt} + q^{t} \right\} \; / \; 2
213  \end{split}
214\]
215As for space operator,
216the adjoint of the derivation and averaging time operators are $\delta_t^*=\delta_{t+\rdt/2}$ and
217$\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$, respectively.
218
219The Leap-frog time stepping given by \autoref{eq:DOM_nxt} can be defined as:
220\[
221  % \label{eq:ALGOS_LF}
222  \frac{\partial q}{\partial t}
223  \equiv \frac{1}{\rdt} \overline{ \delta_{t+\rdt/2}[q]}^{\,t}
224  =         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt}
225\]
226Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$,
227not $2\rdt$ as it can be found sometimes in literature.
228The leap-Frog time stepping is a second order centered scheme.
229As such it respects the quadratic invariant in integral forms, \ie\ the following continuous property,
230\[
231  % \label{eq:ALGOS_Energy}
232  \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}
233  =\int_{t_0}^{t_1} {\frac{1}{2}\, \frac{\partial q^2}{\partial t} \;dt}
234  =  \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) ,
235\]
236is satisfied in discrete form.
237Indeed,
238\[
239  \begin{split}
240    \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}
241    &\equiv \sum\limits_{0}^{N}
242    {\frac{1}{\rdt} q^t \ \overline{ \delta_{t+\rdt/2}[q]}^{\,t} \ \rdt}
243    \equiv \sum\limits_{0}^{N}  { q^t \ \overline{ \delta_{t+\rdt/2}[q]}^{\,t} } \\
244    &\equiv \sum\limits_{0}^{N}  { \overline{q}^{\,t+\Delta/2}{ \delta_{t+\rdt/2}[q]}}
245    \equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta_{t+\rdt/2}[q^2] }\\
246    &\equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta_{t+\rdt/2}[q^2] }
247    \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right)
248  \end{split}
249\]
250NB here pb of boundary condition when applying the adjoint!
251In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition
252(equivalently of the boundary value of the integration by part).
253In time this boundary condition is not physical and \textbf{add something here!!!}
254
255%% =================================================================================================
256\section{Lateral diffusion operator}
257
258%% =================================================================================================
259\subsection{Griffies iso-neutral diffusion operator}
260
261Let try to define a scheme that get its inspiration from the \citet{griffies.gnanadesikan.ea_JPO98} scheme,
262but is formulated within the \NEMO\ framework
263(\ie\ using scale factors rather than grid-size and having a position of $T$-points that
264is not necessary in the middle of vertical velocity points, see \autoref{fig:DOM_zgr_e3}).
265
266In the formulation \autoref{eq:TRA_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO,
267the off-diagonal terms of the small angle diffusion tensor contain several double spatial averages of a gradient,
268for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$.
269It is apparent that the combination of a $k$ average and a $k$ derivative of the tracer allows for
270the presence of grid point oscillation structures that will be invisible to the operator.
271These structures are \textit{computational modes}.
272They will not be damped by the iso-neutral operator, and even possibly amplified by it.
273In other word, the operator applied to a tracer does not warranties the decrease of its global average variance.
274To circumvent this, we have introduced a smoothing of the slopes of the iso-neutral surfaces
275(see \autoref{chap:LDF}).
276Nevertheless, this technique works fine for $T$ and $S$ as they are active tracers
277(\ie\ they enter the computation of density), but it does not work for a passive tracer.
278\citep{griffies.gnanadesikan.ea_JPO98} introduce a different way to discretise the off-diagonal terms that
279nicely solve the problem.
280The idea is to get rid of combinations of an averaged in one direction combined with
281a derivative in the same direction by considering triads.
282For example in the (\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows:
283\begin{equation}
284  \label{eq:ALGOS_Gf_triads}
285  _i^k \mathbb{T}_{i_p}^{k_p} (T)
286  = \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k     \left(
287    \frac{ \delta_{i + i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
288    -\ {_i^k \mathbb{R}_{i_p}^{k_p}} \ \frac{ \delta_{k+k_p} [T^i] }{ {e_{3w}}_{\,i}^{\,k+k_p} }
289  \right)
290\end{equation}
291where the indices $i_p$ and $k_p$ define the four triads and take the following value:
292$i_p = -1/2$ or $1/2$ and $k_p = -1/2$ or $1/2$,
293$b_u= e_{1u}\,e_{2u}\,e_{3u}$ is the volume of $u$-cells,
294$A_i^k$ is the lateral eddy diffusivity coefficient defined at $T$-point,
295and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad:
296\begin{equation}
297  \label{eq:ALGOS_Gf_slopes}
298  _i^k \mathbb{R}_{i_p}^{k_p}
299  =\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} \ \frac
300  {\left(\alpha / \beta \right)_i^\ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] }
301  {\left(\alpha / \beta \right)_i^\ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }
302\end{equation}
303Note that in \autoref{eq:ALGOS_Gf_slopes} we use the ratio $\alpha / \beta$ instead of
304multiplying the temperature derivative by $\alpha$ and the salinity derivative by $\beta$.
305This is more efficient as the ratio $\alpha / \beta$ can to be evaluated directly.
306
307Note that in \autoref{eq:ALGOS_Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of ${b_{uw}}_{\,i+i_p}^{\,k+k_p}$.
308This choice has been motivated by the decrease of tracer variance and
309the presence of partial cell at the ocean bottom (see \autoref{subsec:ALGOS_Gf_operator}).
310
311\begin{figure}[!ht]
312  \centering
313  %\includegraphics[width=0.66\textwidth]{ALGOS_ISO_triad}
314  \caption[Triads used in the Griffies's like iso-neutral diffision scheme for
315    $u$- and $w$-components)]{
316    Triads used in the Griffies's like iso-neutral diffision scheme for
317    $u$-component (upper panel) and $w$-component (lower panel).}
318  \label{fig:ALGOS_ISO_triad}
319\end{figure}
320
321The four iso-neutral fluxes associated with the triads are defined at $T$-point.
322They take the following expression:
323\begin{flalign*}
324  % \label{eq:ALGOS_Gf_fluxes}
325  \begin{split}
326    {_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)
327    &= \ \; \qquad  \quad    { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+i_p}^{\,k}}    \\
328    {_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T)
329    &=  -\; { _i^k \mathbb{R}_{i_p}^{k_p} }
330    \ \; { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+k_p}}
331  \end{split}
332\end{flalign*}
333
334The resulting iso-neutral fluxes at $u$- and $w$-points are then given by
335the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:TRIADS_ISO_triad}):
336\begin{flalign}
337  \label{eq:ALGOS_iso_flux}
338  \textbf{F}_{iso}(T)
339  &\equiv  \sum_{\substack{i_p,\,k_p}}
340  \begin{pmatrix}
341    {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) \\ \\
342    {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)
343  \end{pmatrix}
344  \notag \\
345  &  \notag \\
346  &\equiv  \sum_{\substack{i_p,\,k_p}}
347  \begin{pmatrix}
348    && { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+1/2}^{\,k} } \\ \\
349    & -\; { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }
350    & {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+1/2} }
351  \end{pmatrix}      % \\
352  % &\\
353  % &\equiv  \sum_{\substack{i_p,\,k_p}}
354  % \begin{pmatrix}
355  %   \qquad  \qquad  \qquad
356  %   \frac{1}{ {e_{1u}}_{\,i+1/2}^{\,k} }  \ \;
357  %   { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T)\\
358  %   \\
359  %   -\frac{1}{ {e_{3w}}_{\,i}^{\,k+1/2} }  \ \;
360  %   { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \;
361  %   {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T)\\
362  % \end{pmatrix}
363\end{flalign}
364resulting in a iso-neutral diffusion tendency on temperature given by
365the divergence of the sum of all the four triad fluxes:
366\begin{equation}
367  \label{eq:ALGOS_Gf_operator}
368  D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{
369    \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
370    + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]   \right\}
371\end{equation}
372where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells.
373
374This expression of the iso-neutral diffusion has been chosen in order to satisfy the following six properties:
375\begin{description}
376\item [Horizontal diffusion] The discretization of the diffusion operator recovers the traditional five-point Laplacian in the limit of flat iso-neutral direction:
377  \[
378    % \label{eq:ALGOS_Gf_property1a}
379    D_l^T = \frac{1}{b_T}  \ \delta_{i}
380    \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right]
381    \qquad  \text{when} \quad
382    { _i^k \mathbb{R}_{i_p}^{k_p} }=0
383  \]
384\item [Implicit treatment in the vertical] In the diagonal term associated with the vertical divergence of the iso-neutral fluxes
385  \ie\ the term associated with a second order vertical derivative)
386  appears only tracer values associated with a single water column.
387  This is of paramount importance since it means that
388  the implicit in time algorithm for solving the vertical diffusion equation can be used to evaluate this term.
389  It is a necessity since the vertical eddy diffusivity associated with this term,
390  \[
391    \sum_{\substack{i_p, \,k_p}} \left\{
392      A_i^k \; \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
393    \right\}
394  \]
395  can be quite large.
396\item [Pure iso-neutral operator] The iso-neutral flux of locally referenced potential density is zero, \ie
397  \begin{align*}
398    % \label{eq:ALGOS_Gf_property2}
399    \begin{matrix}
400      &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} (\rho)}
401      &=    &\alpha_i^k   &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)
402      &- \ \;  \beta _i^k    &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (S) & = \ 0   \\
403      &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)}
404      &=    &\alpha_i^k   &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T)
405      &- \  \; \beta _i^k    &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (S)  &= \ 0
406    \end{matrix}
407  \end{align*}
408  This result is trivially obtained using the \autoref{eq:ALGOS_Gf_triads} applied to $T$ and $S$ and
409  the definition of the triads' slopes \autoref{eq:ALGOS_Gf_slopes}.
410\item [Conservation of tracer] The iso-neutral diffusion term conserve the total tracer content, \ie
411  \[
412    % \label{eq:ALGOS_Gf_property1}
413    \sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0
414  \]
415This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form.
416\item [Decrease of tracer variance] The iso-neutral diffusion term does not increase the total tracer variance, \ie
417  \[
418    % \label{eq:ALGOS_Gf_property1}
419    \sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0
420  \]
421The property is demonstrated in the \autoref{subsec:ALGOS_Gf_operator}.
422It is a key property for a diffusion term.
423It means that the operator is also a dissipation term,
424\ie\ it is a sink term for the square of the quantity on which it is applied.
425It therfore ensures that, when the diffusivity coefficient is large enough,
426the field on which it is applied become free of grid-point noise.
427\item [Self-adjoint operator] The iso-neutral diffusion operator is self-adjoint, \ie
428  \[
429    % \label{eq:ALGOS_Gf_property1}
430    \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}
431  \]
432In other word, there is no needs to develop a specific routine from the adjoint of this operator.
433We just have to apply the same routine.
434This properties can be demonstrated quite easily in a similar way the "non increase of tracer variance" property
435has been proved (see \autoref{apdx:ALGOS_Gf_operator}).
436\end{description}
437
438%% =================================================================================================
439\subsection{Eddy induced velocity and skew flux formulation}
440
441When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined),
442an additional advection term is added.
443The associated velocity is the so called eddy induced velocity,
444the formulation of which depends on the slopes of iso-neutral surfaces.
445Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces,
446\ie\ \autoref{eq:LDF_slp_geo} is used in $z$-coordinate,
447and the sum \autoref{eq:LDF_slp_geo} + \autoref{eq:LDF_slp_iso} in $z^*$ or $s$-coordinates.
448
449The eddy induced velocity is given by:
450\begin{equation}
451  \label{eq:ALGOS_eiv_v}
452  \begin{split}
453    u^* & = - \frac{1}{e_2\,e_{3}}          \;\partial_k \left( e_2 \, A_e \; r_\right)
454    = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_\right)            \\
455    v^* & = - \frac{1}{e_1\,e_3}\;             \partial_k \left( e_1 \, A_e \; r_\right)
456    = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_\right)             \\
457    w^* & =    \frac{1}{e_1\,e_2}\; \left\{   \partial_\left( e_2 \, A_e \; r_\right)
458      + \partial_\left( e_1 \, A_e \;r_j   \right) \right\}
459  \end{split}
460\end{equation}
461where $A_{e}$ is the eddy induced velocity coefficient,
462and $r_i$ and $r_j$ the slopes between the iso-neutral and the geopotential surfaces.
463\cmtgm{Wrong: to be modified with 2 2D streamfunctions}
464In other words, the eddy induced velocity can be derived from a vector streamfuntion, $\phi$,
465which is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^*  = \textbf{k} \times \nabla \phi$.
466
467A traditional way to implement this additional advection is to add it to the eulerian velocity prior to
468compute the tracer advection.
469This allows us to take advantage of all the advection schemes offered for the tracers
470(see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme.
471This is particularly useful for passive tracers where
472\emph{positivity} of the advection scheme is of paramount importance.
473% give here the expression using the triads. It is different from the one given in \autoref{eq:LDF_eiv}
474% see just below a copy of this equation:
475%\begin{equation} \label{eq:ALGOS_ldfeiv}
476%\begin{split}
477% u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\
478% v^* & = \frac{1}{e_{1u}e_{3v}}\; \delta_k \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right]\\
479%w^* & = \frac{1}{e_{1w}e_{2w}}\; \left\{ \delta_i \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right] + %\delta_j \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right] \right\} \\
480%\end{split}
481%\end{equation}
482\[
483  % \label{eq:ALGOS_eiv_vd}
484  \textbf{F}_{eiv}^T   \equiv   \left(
485    \begin{aligned}
486      \sum_{\substack{i_p,\,k_p}} &
487      +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k}
488      \ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ \\
489      \sum_{\substack{i_p,\,k_p}} &
490      - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p}
491      \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]
492    \end{aligned}
493  \right)
494\]
495
496\citep{griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form.
497It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity.
498For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows:
499\begin{flalign*}
500  \begin{split}
501    \textbf{F}_{eiv}^T =
502    \begin{pmatrix}
503      {e_{2}\,e_{3}\;  u^*}      \\
504      {e_{1}\,e_{2}\; w^*}
505    \end{pmatrix}
506    \;   T
507    &=
508    \begin{pmatrix}
509      { - \partial_k \left( e_{2} \, A_{e} \; r_i \right) \; T \;}      \\
510      {+ \partial_\left( e_{2} \, A_{e} \; r_i \right) \; T \;}
511    \end{pmatrix}
512    \\
513    &=
514    \begin{pmatrix}
515      { - \partial_k \left( e_{2} \, A_{e} \; r_\; T \right) \;}  \\
516      {+ \partial_\left( e_{2} \, A_{e} \; r_\; T \right) \;}
517    \end{pmatrix}
518    +
519    \begin{pmatrix}
520      {+ e_{2} \, A_{e} \; r_\; \partial_k T}  \\
521      { - e_{2} \, A_{e} \; r_\; \partial_i  T}
522    \end{pmatrix}
523  \end{split}
524\end{flalign*}
525and since the eddy induces velocity field is no-divergent,
526we end up with the skew form of the eddy induced advective fluxes:
527\begin{equation}
528  \label{eq:ALGOS_eiv_skew_continuous}
529  \textbf{F}_{eiv}^T =
530  \begin{pmatrix}
531    {+ e_{2} \, A_{e} \; r_\; \partial_k T}   \\
532    { - e_{2} \, A_{e} \; r_\; \partial_i  T}
533  \end{pmatrix}
534\end{equation}
535The tendency associated with eddy induced velocity is then simply the divergence of
536the \autoref{eq:ALGOS_eiv_skew_continuous} fluxes.
537It naturally conserves the tracer content, as it is expressed in flux form and,
538as the advective form, it preserves the tracer variance.
539Another interesting property of \autoref{eq:ALGOS_eiv_skew_continuous} form is that when $A=A_e$,
540a simplification occurs in the sum of the iso-neutral diffusion and eddy induced velocity terms:
541\begin{flalign*}
542  % \label{eq:ALGOS_eiv_skew+eiv_continuous}
543  \textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &=
544  \begin{pmatrix}
545    + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T -  e_2 \, A \; r_i                              \;\partial_k T   \\
546    -  e_2 \, A_{e} \; r_i           \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T
547  \end{pmatrix}
548  +
549  \begin{pmatrix}
550    {+ e_{2} \, A_{e} \; r_\; \partial_k T}   \\
551    { - e_{2} \, A_{e} \; r_\; \partial_i  T}
552  \end{pmatrix}
553  \\
554  &=
555  \begin{pmatrix}
556    + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T    \\
557    -  2\; e_2 \, A_{e} \; r_i      \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T
558  \end{pmatrix}
559\end{flalign*}
560The horizontal component reduces to the one use for an horizontal laplacian operator and
561the vertical one keeps the same complexity, but not more.
562This property has been used to reduce the computational time \citep{griffies_JPO98},
563but it is not of practical use as usually $A \neq A_e$.
564Nevertheless this property can be used to choose a discret form of \autoref{eq:ALGOS_eiv_skew_continuous} which
565is consistent with the iso-neutral operator \autoref{eq:ALGOS_Gf_operator}.
566Using the slopes \autoref{eq:ALGOS_Gf_slopes} and defining $A_e$ at $T$-point(\ie\ as $A$,
567the eddy diffusivity coefficient), the resulting discret form is given by:
568\begin{equation}
569  \label{eq:ALGOS_eiv_skew}
570  \textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left(
571    \begin{aligned}
572      \sum_{\substack{i_p,\,k_p}} &
573      +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k}
574      \ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ \\
575      \sum_{\substack{i_p,\,k_p}} &
576      - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p}
577      \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]
578    \end{aligned}
579  \right)
580\end{equation}
581Note that \autoref{eq:ALGOS_eiv_skew} is valid in $z$-coordinate with or without partial cells.
582In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces must be added to
583$\mathbb{R}$ for the discret form to be exact.
584
585Such a choice of discretisation is consistent with the iso-neutral operator as
586it uses the same definition for the slopes.
587It also ensures the conservation of the tracer variance (see \autoref{subsec:ALGOS_eiv_skew}),
588\ie\ it does not include a diffusive component but is a "pure" advection term.
589
590%% =================================================================================================
591\subsection{Discrete invariants of the iso-neutral diffrusion}
592\label{subsec:ALGOS_Gf_operator}
593
594Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane.
595
596This part will be moved in an Appendix.
597
598The continuous property to be demonstrated is:
599\[
600  \int_D  D_l^T \; T \;dv   \leq 0
601\]
602The discrete form of its left hand side is obtained using \autoref{eq:TRIADS_iso_flux}
603
604\begin{align*}
605  &\int_D  D_l^T \; T \;dv \equiv  \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\}    \\
606  &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
607    \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
608    + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]  \ T \right\}    \\
609  &\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
610    {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T]
611    + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\
612  &\equiv -\sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
613    \frac{ _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} (T) }{ {e_{1u}}_{\,i+1/2}^{\,k} }  \ \delta_{i+1/2} [T]
614    - { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \;
615    \frac{ _i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} (T) }{ {e_{3w}}_{\,i}^{\,k+1/2}  } \ \delta_{k+1/2} [T]
616    \right\}      \\
617    %
618  \allowdisplaybreaks
619  \intertext{ Expending the summation on $i_p$ and $k_p$, it becomes:}
620  %
621  &\equiv -\sum_{i,k}
622    \begin{Bmatrix}
623      &\ \ \Bigl{ _{i+1}^{k} \mathbb{T}_{-1/2}^{-1/2} (T) }
624      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
625      & -\ \ {_{i}^{k+1} \mathbb{R}_{-1/2}^{-1/2}}
626      &      {_{i}^{k+1} \mathbb{T}_{-1/2}^{-1/2} (T) }
627      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
628      & \\
629      &+\Bigl( \ \;\; { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }
630      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
631      & -\ \ {_i^{k+1} \mathbb{R}_{+1/2}^{-1/2}}
632      & { _i^{k+1} \mathbb{T}_{+1/2}^{-1/2} (T) }
633      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}      \Bigr)
634      & \\
635      &+\Bigl{ _{i+1}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }
636      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
637      & -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{-1/2}^{+1/2}}
638      &      \ \;\;{_{i}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }
639      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
640      & \\
641      &+\Bigl( \ \;\; { _{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }
642      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
643      & -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{+1/2}^{+1/2}}
644      &      \ \;\;{_{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }
645      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
646    \end{Bmatrix}
647    %
648  \allowdisplaybreaks
649  \intertext{
650  The summation is done over all $i$ and $k$ indices,
651  it is therefore possible to introduce a shift of $-1$ either in $i$ or $k$ direction in order to
652  regroup all the terms of the summation by triad at a ($i$,$k$) point.
653  In other words, we regroup all the terms in the neighbourhood that contain a triad at the same ($i$,$k$) indices.
654  It becomes:
655  }
656  %
657  &\equiv -\sum_{i,k}
658    \begin{Bmatrix}
659      &\ \ \Bigl{_i^k \mathbb{T}_{-1/2}^{-1/2} (T) }
660      &\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
661      & -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}}
662      &      {_i^k \mathbb{T}_{-1/2}^{-1/2} (T) }
663      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}     \Bigr)
664      & \\
665      &+\Bigl{ _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }
666      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
667      & -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}}
668      &      { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }
669      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}      \Bigr)
670      & \\
671      &+\Bigl{_i^k \mathbb{T}_{-1/2}^{+1/2} (T) }
672      &\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
673      & -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}}
674      &      {_i^k \mathbb{T}_{-1/2}^{+1/2} (T) }
675      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
676      & \\
677      &+\Bigl( { _i^k \mathbb{T}_{+1/2}^{+1/2} (T) }
678      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
679      & -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}}
680      &      {_i^k \mathbb{T}_{+1/2}^{+1/2} (T) }
681      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
682    \end{Bmatrix}   \\
683    %
684  \allowdisplaybreaks
685  \intertext{
686  Then outing in factor the triad in each of the four terms of the summation and
687  substituting the triads by their expression given in \autoref{eq:ALGOS_Gf_triads}.
688  It becomes:
689  }
690  %
691  &\equiv -\sum_{i,k}
692    \begin{Bmatrix}
693      &\ \ \Bigl\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
694      & -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}}
695      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}     \Bigr)^2
696      & \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k}  \  A_i^k
697      & \\
698      &+\Bigl\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
699      & -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}}
700      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}      \Bigr)^2
701      & \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k}  \  A_i^k
702      & \\
703      &+\Bigl\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
704      & -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}}
705      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)^2
706      & \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k}  \  A_i^k
707      & \\
708      &+\Bigl( \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
709      & -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}}
710      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)^2
711      & \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k}  \  A_i^k      \\
712    \end{Bmatrix}
713  \\
714  & \\
715  %
716  &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
717    \begin{matrix}
718      &\Bigl( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}}
719      & -\ \ {_i^k \mathbb{R}_{i_p}^{k_p}}
720      &\frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \Bigr)^2
721      & \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \ \
722    \end{matrix}
723        \right\}
724        \quad   \leq 0
725\end{align*}
726The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities.
727
728Note that, if instead of multiplying $D_l^T$ by $T$, we were using another tracer field, let say $S$,
729then the previous demonstration would have let to:
730\begin{align*}
731  \int_D  S \; D_l^\;dv &\equiv  \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\}    \\
732                           &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
733                             \left( \frac{ \delta_{i +i_p} [S] }{{e_{1u} }_{\,i+i_p}^{\,k}}
734                             - {_i^k \mathbb{R}_{i_p}^{k_p}}
735                             \frac{ \delta_{k+k_p} [S] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right\right. \\
736                           & \qquad \qquad \qquad \ \left.
737                             \left( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}}
738                             - {_i^k \mathbb{R}_{i_p}^{k_p}}
739                             \frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right)
740                             \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \
741                             \right\}
742                             %
743                             \allowdisplaybreaks
744                             \intertext{
745                             which, by applying the same operation as before but in reverse order, leads to:
746                             }
747                             %
748                           &\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\}
749\end{align*}
750This means that the iso-neutral operator is self-adjoint.
751There is no need to develop a specific to obtain it.
752
753%% =================================================================================================
754\subsection{Discrete invariants of the skew flux formulation}
755\label{subsec:ALGOS_eiv_skew}
756
757Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.
758
759This have to be moved in an Appendix.
760
761The continuous property to be demonstrated is:
762\begin{align*}
763  \int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0
764\end{align*}
765The discrete form of its left hand side is obtained using \autoref{eq:ALGOS_eiv_skew}
766\begin{align*}
767  \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
768  \delta_&\left[
769              {e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}
770              \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]
771              \right] \; T_i^k      \\
772  - \delta_k &\left[
773               {e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}
774               \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]
775               \right] \; T_i^k      \         \Biggr\}
776\end{align*}
777apply the adjoint of delta operator, it becomes
778\begin{align*}
779  \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
780  &\left(
781    {e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}
782    \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]
783    \right) \; \delta_{i+1/2}[T^{k}]      \\
784  - &\left(
785      {e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}
786      \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]
787      \right) \; \delta_{k+1/2}[T_{i}]       \         \Biggr\}
788\end{align*}
789Expending the summation on $i_p$ and $k_p$, it becomes:
790\begin{align*}
791  \begin{matrix}
792    &\sum\limits_{i,k}   \Bigl\{
793    &+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}
794    &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}]    &\delta_{i+1/2}[T^{k}]   &\\
795    &&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}
796    &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}}  &\delta_{k-1/2}[T_{i\ \ \ \;}&\delta_{i+1/2}[T^{k}] &\\
797    &&+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}
798    &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}]     &\delta_{i+1/2}[T^{k}] &\\
799    &&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}
800    &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\
801    %
802    &&-{e_{2u}}_i^{k+1}                                &{A_{e}}_i^{k+1}
803    &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}}   &\delta_{i-1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\
804    &&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}
805    &{\ \ \;_i^\mathbb{R}_{-1/2}^{+1/2}}   &\delta_{i-1/2}[T^{k\ \ \ \:}&\delta_{k+1/2}[T_{i}] &\\
806    &&-{e_{2u}}_i^{k+1    }                             &{A_{e}}_i^{k+1}
807    &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}}   &\delta_{i+1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\
808    &&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}
809    &{\ \ \;_i^\mathbb{R}_{+1/2}^{+1/2}}   &\delta_{i+1/2}[T^{k\ \ \ \:}&\delta_{k+1/2}[T_{i}]
810    &\Bigr\}  \\
811  \end{matrix}
812\end{align*}
813The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs,
814they cancel out.
815Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$.
816The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the same but both of opposite signs and
817shifted by 1 in $k$ direction.
818When summing over $k$ they cancel out with the neighbouring grid points.
819Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the $i$ direction.
820Therefore the sum over the domain is zero,
821\ie\ the variance of the tracer is preserved by the discretisation of the skew fluxes.
822
823\subinc{\input{../../global/epilogue}}
824
825\end{document}
Note: See TracBrowser for help on using the repository browser.