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

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

Update chapters according to previous commit

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%%gm 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%%end gm
467
468A traditional way to implement this additional advection is to add it to the eulerian velocity prior to
469compute the tracer advection.
470This allows us to take advantage of all the advection schemes offered for the tracers
471(see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme.
472This is particularly useful for passive tracers where
473\emph{positivity} of the advection scheme is of paramount importance.
474% give here the expression using the triads. It is different from the one given in \autoref{eq:LDF_eiv}
475% see just below a copy of this equation:
476%\begin{equation} \label{eq:ALGOS_ldfeiv}
477%\begin{split}
478% u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\
479% v^* & = \frac{1}{e_{1u}e_{3v}}\; \delta_k \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right]\\
480%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\} \\
481%\end{split}
482%\end{equation}
483\[
484  % \label{eq:ALGOS_eiv_vd}
485  \textbf{F}_{eiv}^T   \equiv   \left(
486    \begin{aligned}
487      \sum_{\substack{i_p,\,k_p}} &
488      +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k}
489      \ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ \\
490      \sum_{\substack{i_p,\,k_p}} &
491      - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p}
492      \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]
493    \end{aligned}
494  \right)
495\]
496
497\citep{griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form.
498It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity.
499For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows:
500\begin{flalign*}
501  \begin{split}
502    \textbf{F}_{eiv}^T =
503    \begin{pmatrix}
504      {e_{2}\,e_{3}\;  u^*}      \\
505      {e_{1}\,e_{2}\; w^*}
506    \end{pmatrix}
507    \;   T
508    &=
509    \begin{pmatrix}
510      { - \partial_k \left( e_{2} \, A_{e} \; r_i \right) \; T \;}      \\
511      {+ \partial_\left( e_{2} \, A_{e} \; r_i \right) \; T \;}
512    \end{pmatrix}
513    \\
514    &=
515    \begin{pmatrix}
516      { - \partial_k \left( e_{2} \, A_{e} \; r_\; T \right) \;}  \\
517      {+ \partial_\left( e_{2} \, A_{e} \; r_\; T \right) \;}
518    \end{pmatrix}
519    +
520    \begin{pmatrix}
521      {+ e_{2} \, A_{e} \; r_\; \partial_k T}  \\
522      { - e_{2} \, A_{e} \; r_\; \partial_i  T}
523    \end{pmatrix}
524  \end{split}
525\end{flalign*}
526and since the eddy induces velocity field is no-divergent,
527we end up with the skew form of the eddy induced advective fluxes:
528\begin{equation}
529  \label{eq:ALGOS_eiv_skew_continuous}
530  \textbf{F}_{eiv}^T =
531  \begin{pmatrix}
532    {+ e_{2} \, A_{e} \; r_\; \partial_k T}   \\
533    { - e_{2} \, A_{e} \; r_\; \partial_i  T}
534  \end{pmatrix}
535\end{equation}
536The tendency associated with eddy induced velocity is then simply the divergence of
537the \autoref{eq:ALGOS_eiv_skew_continuous} fluxes.
538It naturally conserves the tracer content, as it is expressed in flux form and,
539as the advective form, it preserves the tracer variance.
540Another interesting property of \autoref{eq:ALGOS_eiv_skew_continuous} form is that when $A=A_e$,
541a simplification occurs in the sum of the iso-neutral diffusion and eddy induced velocity terms:
542\begin{flalign*}
543  % \label{eq:ALGOS_eiv_skew+eiv_continuous}
544  \textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &=
545  \begin{pmatrix}
546    + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T -  e_2 \, A \; r_i                              \;\partial_k T   \\
547    -  e_2 \, A_{e} \; r_i           \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T
548  \end{pmatrix}
549  +
550  \begin{pmatrix}
551    {+ e_{2} \, A_{e} \; r_\; \partial_k T}   \\
552    { - e_{2} \, A_{e} \; r_\; \partial_i  T}
553  \end{pmatrix}
554  \\
555  &=
556  \begin{pmatrix}
557    + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T    \\
558    -  2\; e_2 \, A_{e} \; r_i      \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T
559  \end{pmatrix}
560\end{flalign*}
561The horizontal component reduces to the one use for an horizontal laplacian operator and
562the vertical one keeps the same complexity, but not more.
563This property has been used to reduce the computational time \citep{griffies_JPO98},
564but it is not of practical use as usually $A \neq A_e$.
565Nevertheless this property can be used to choose a discret form of \autoref{eq:ALGOS_eiv_skew_continuous} which
566is consistent with the iso-neutral operator \autoref{eq:ALGOS_Gf_operator}.
567Using the slopes \autoref{eq:ALGOS_Gf_slopes} and defining $A_e$ at $T$-point(\ie\ as $A$,
568the eddy diffusivity coefficient), the resulting discret form is given by:
569\begin{equation}
570  \label{eq:ALGOS_eiv_skew}
571  \textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left(
572    \begin{aligned}
573      \sum_{\substack{i_p,\,k_p}} &
574      +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k}
575      \ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ \\
576      \sum_{\substack{i_p,\,k_p}} &
577      - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p}
578      \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]
579    \end{aligned}
580  \right)
581\end{equation}
582Note that \autoref{eq:ALGOS_eiv_skew} is valid in $z$-coordinate with or without partial cells.
583In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces must be added to
584$\mathbb{R}$ for the discret form to be exact.
585
586Such a choice of discretisation is consistent with the iso-neutral operator as
587it uses the same definition for the slopes.
588It also ensures the conservation of the tracer variance (see \autoref{subsec:ALGOS_eiv_skew}),
589\ie\ it does not include a diffusive component but is a "pure" advection term.
590
591%% =================================================================================================
592\subsection{Discrete invariants of the iso-neutral diffrusion}
593\label{subsec:ALGOS_Gf_operator}
594
595Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane.
596
597This part will be moved in an Appendix.
598
599The continuous property to be demonstrated is:
600\[
601  \int_D  D_l^T \; T \;dv   \leq 0
602\]
603The discrete form of its left hand side is obtained using \autoref{eq:TRIADS_iso_flux}
604
605\begin{align*}
606  &\int_D  D_l^T \; T \;dv \equiv  \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\}    \\
607  &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
608    \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
609    + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]  \ T \right\}    \\
610  &\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
611    {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T]
612    + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\
613  &\equiv -\sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
614    \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]
615    - { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \;
616    \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]
617    \right\}      \\
618    %
619  \allowdisplaybreaks
620  \intertext{ Expending the summation on $i_p$ and $k_p$, it becomes:}
621  %
622  &\equiv -\sum_{i,k}
623    \begin{Bmatrix}
624      &\ \ \Bigl{ _{i+1}^{k} \mathbb{T}_{-1/2}^{-1/2} (T) }
625      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
626      & -\ \ {_{i}^{k+1} \mathbb{R}_{-1/2}^{-1/2}}
627      &      {_{i}^{k+1} \mathbb{T}_{-1/2}^{-1/2} (T) }
628      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
629      & \\
630      &+\Bigl( \ \;\; { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }
631      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
632      & -\ \ {_i^{k+1} \mathbb{R}_{+1/2}^{-1/2}}
633      & { _i^{k+1} \mathbb{T}_{+1/2}^{-1/2} (T) }
634      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}      \Bigr)
635      & \\
636      &+\Bigl{ _{i+1}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }
637      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
638      & -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{-1/2}^{+1/2}}
639      &      \ \;\;{_{i}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }
640      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
641      & \\
642      &+\Bigl( \ \;\; { _{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }
643      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
644      & -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{+1/2}^{+1/2}}
645      &      \ \;\;{_{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }
646      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
647    \end{Bmatrix}
648    %
649  \allowdisplaybreaks
650  \intertext{
651  The summation is done over all $i$ and $k$ indices,
652  it is therefore possible to introduce a shift of $-1$ either in $i$ or $k$ direction in order to
653  regroup all the terms of the summation by triad at a ($i$,$k$) point.
654  In other words, we regroup all the terms in the neighbourhood that contain a triad at the same ($i$,$k$) indices.
655  It becomes:
656  }
657  %
658  &\equiv -\sum_{i,k}
659    \begin{Bmatrix}
660      &\ \ \Bigl{_i^k \mathbb{T}_{-1/2}^{-1/2} (T) }
661      &\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
662      & -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}}
663      &      {_i^k \mathbb{T}_{-1/2}^{-1/2} (T) }
664      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}     \Bigr)
665      & \\
666      &+\Bigl{ _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }
667      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
668      & -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}}
669      &      { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }
670      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}      \Bigr)
671      & \\
672      &+\Bigl{_i^k \mathbb{T}_{-1/2}^{+1/2} (T) }
673      &\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
674      & -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}}
675      &      {_i^k \mathbb{T}_{-1/2}^{+1/2} (T) }
676      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
677      & \\
678      &+\Bigl( { _i^k \mathbb{T}_{+1/2}^{+1/2} (T) }
679      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
680      & -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}}
681      &      {_i^k \mathbb{T}_{+1/2}^{+1/2} (T) }
682      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
683    \end{Bmatrix}   \\
684    %
685  \allowdisplaybreaks
686  \intertext{
687  Then outing in factor the triad in each of the four terms of the summation and
688  substituting the triads by their expression given in \autoref{eq:ALGOS_Gf_triads}.
689  It becomes:
690  }
691  %
692  &\equiv -\sum_{i,k}
693    \begin{Bmatrix}
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      & \\
709      &+\Bigl( \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
710      & -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}}
711      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)^2
712      & \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k}  \  A_i^k      \\
713    \end{Bmatrix}
714  \\
715  & \\
716  %
717  &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
718    \begin{matrix}
719      &\Bigl( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}}
720      & -\ \ {_i^k \mathbb{R}_{i_p}^{k_p}}
721      &\frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \Bigr)^2
722      & \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \ \
723    \end{matrix}
724        \right\}
725        \quad   \leq 0
726\end{align*}
727The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities.
728
729Note that, if instead of multiplying $D_l^T$ by $T$, we were using another tracer field, let say $S$,
730then the previous demonstration would have let to:
731\begin{align*}
732  \int_D  S \; D_l^\;dv &\equiv  \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\}    \\
733                           &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
734                             \left( \frac{ \delta_{i +i_p} [S] }{{e_{1u} }_{\,i+i_p}^{\,k}}
735                             - {_i^k \mathbb{R}_{i_p}^{k_p}}
736                             \frac{ \delta_{k+k_p} [S] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right\right. \\
737                           & \qquad \qquad \qquad \ \left.
738                             \left( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}}
739                             - {_i^k \mathbb{R}_{i_p}^{k_p}}
740                             \frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right)
741                             \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \
742                             \right\}
743                             %
744                             \allowdisplaybreaks
745                             \intertext{
746                             which, by applying the same operation as before but in reverse order, leads to:
747                             }
748                             %
749                           &\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\}
750\end{align*}
751This means that the iso-neutral operator is self-adjoint.
752There is no need to develop a specific to obtain it.
753
754%% =================================================================================================
755\subsection{Discrete invariants of the skew flux formulation}
756\label{subsec:ALGOS_eiv_skew}
757
758Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.
759
760This have to be moved in an Appendix.
761
762The continuous property to be demonstrated is:
763\begin{align*}
764  \int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0
765\end{align*}
766The discrete form of its left hand side is obtained using \autoref{eq:ALGOS_eiv_skew}
767\begin{align*}
768  \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
769  \delta_&\left[
770              {e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}
771              \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]
772              \right] \; T_i^k      \\
773  - \delta_k &\left[
774               {e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}
775               \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]
776               \right] \; T_i^k      \         \Biggr\}
777\end{align*}
778apply the adjoint of delta operator, it becomes
779\begin{align*}
780  \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
781  &\left(
782    {e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}
783    \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]
784    \right) \; \delta_{i+1/2}[T^{k}]      \\
785  - &\left(
786      {e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}
787      \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]
788      \right) \; \delta_{k+1/2}[T_{i}]       \         \Biggr\}
789\end{align*}
790Expending the summation on $i_p$ and $k_p$, it becomes:
791\begin{align*}
792  \begin{matrix}
793    &\sum\limits_{i,k}   \Bigl\{
794    &+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}
795    &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}]    &\delta_{i+1/2}[T^{k}]   &\\
796    &&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}
797    &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}}  &\delta_{k-1/2}[T_{i\ \ \ \;}&\delta_{i+1/2}[T^{k}] &\\
798    &&+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}
799    &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}]     &\delta_{i+1/2}[T^{k}] &\\
800    &&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}
801    &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\
802    %
803    &&-{e_{2u}}_i^{k+1}                                &{A_{e}}_i^{k+1}
804    &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}}   &\delta_{i-1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\
805    &&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}
806    &{\ \ \;_i^\mathbb{R}_{-1/2}^{+1/2}}   &\delta_{i-1/2}[T^{k\ \ \ \:}&\delta_{k+1/2}[T_{i}] &\\
807    &&-{e_{2u}}_i^{k+1    }                             &{A_{e}}_i^{k+1}
808    &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}}   &\delta_{i+1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\
809    &&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}
810    &{\ \ \;_i^\mathbb{R}_{+1/2}^{+1/2}}   &\delta_{i+1/2}[T^{k\ \ \ \:}&\delta_{k+1/2}[T_{i}]
811    &\Bigr\}  \\
812  \end{matrix}
813\end{align*}
814The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs,
815they cancel out.
816Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$.
817The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the same but both of opposite signs and
818shifted by 1 in $k$ direction.
819When summing over $k$ they cancel out with the neighbouring grid points.
820Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the $i$ direction.
821Therefore the sum over the domain is zero,
822\ie\ the variance of the tracer is preserved by the discretisation of the skew fluxes.
823
824\onlyinsubfile{\input{../../global/epilogue}}
825
826\end{document}
Note: See TracBrowser for help on using the repository browser.