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

Last change on this file since 10368 was 10368, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10365, see #2133

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