source: branches/2017/dev_merge_2017/DOC/tex_sub/annex_E.tex @ 9407

Last change on this file since 9407 was 9407, checked in by nicolasmartin, 3 years ago

Complete refactoring of cross-referencing

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