New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
apdx_algos.tex in NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

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

Last change on this file since 11597 was 11597, checked in by nicolasmartin, 5 years ago

Continuation of coding rules application
Recovery of some sections deleted by the previous commit

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