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

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

Apply two thirds ratio on figures width as default

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