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

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

Cleaning of section headings, reinstating the index by mixing \np and \forcode macros, continued conversion of source code inclusions

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 
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{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 \eqref{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(\ref{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 \eqref{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\eqref{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 \eqref{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{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{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 \eqref{Eq_DOM_nxt} can be defined as:
205\begin{equation} \label{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 \eqref{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{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 Fig.~\ref{Fig_zgr_e3}).
255
256In the formulation \eqref{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 \S\ref{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{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{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 \eqref{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 \eqref{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 Appendix~\ref{Apdx_Gf_operator}).
299
300%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
301\begin{figure}[!ht] \label{Fig_ISO_triad}
302\begin{center}
303\includegraphics[width=0.70\textwidth]{Fig_ISO_triad}
304\caption{  \protect\label{Fig_ISO_triad}   
305Triads used in the Griffies's like iso-neutral diffision scheme for
306$u$-component (upper panel) and $w$-component (lower panel).}
307\end{center}
308\end{figure}
309%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
310
311The four iso-neutral fluxes associated with the triads are defined at $T$-point.
312They take the following expression :
313\begin{flalign} \label{Gf_fluxes}
314\begin{split}
315{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)
316   &= \ \; \qquad  \quad    { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+i_p}^{\,k}}    \\
317{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T)
318   &=  -\; { _i^k \mathbb{R}_{i_p}^{k_p} }
319             \ \; { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+k_p}}
320\end{split}
321\end{flalign}
322
323The resulting iso-neutral fluxes at $u$- and $w$-points are then given by the
324sum of the fluxes that cross the $u$- and $w$-face (Fig.~\ref{Fig_ISO_triad}):
325\begin{flalign} \label{Eq_iso_flux} 
326\textbf{F}_{iso}(T)
327&\equiv  \sum_{\substack{i_p,\,k_p}} 
328   \begin{pmatrix} 
329      {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\
330      \\
331      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)      \\   
332   \end{pmatrix}    \notag \\
333&  \notag \\
334&\equiv  \sum_{\substack{i_p,\,k_p}} 
335   \begin{pmatrix} 
336      && { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+1/2}^{\,k} }    \\
337      \\
338      & -\; { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }
339        & {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+1/2} }   \\   
340   \end{pmatrix}      % \\
341% &\\
342% &\equiv  \sum_{\substack{i_p,\,k_p}}
343%    \begin{pmatrix}
344%       \qquad  \qquad  \qquad
345%       \frac{1}{ {e_{1u}}_{\,i+1/2}^{\,k} }  \ \; 
346%        { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T)\\
347%       \\
348%       -\frac{1}{ {e_{3w}}_{\,i}^{\,k+1/2} }  \ \;
349%        { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \; 
350%                  {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T)\\   
351%    \end{pmatrix}     
352\end{flalign}
353resulting in a iso-neutral diffusion tendency on temperature given by the divergence
354of the sum of all the four triad fluxes :
355\begin{equation} \label{Gf_operator}
356D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{ 
357       \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] 
358        + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]   \right\}
359\end{equation}
360where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells.
361
362This expression of the iso-neutral diffusion has been chosen in order to satisfy
363the following six properties:
364\begin{description}
365\item[$\bullet$ horizontal diffusion] The discretization of the diffusion operator
366recovers the traditional five-point Laplacian in the limit of flat iso-neutral direction :
367\begin{equation} \label{Gf_property1a}
368D_l^T = \frac{1}{b_T}  \ \delta_{i} 
369   \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right]
370\qquad  \text{when} \quad 
371   { _i^k \mathbb{R}_{i_p}^{k_p} }=0
372\end{equation}
373
374\item[$\bullet$ implicit treatment in the vertical]  In the diagonal term associated
375with the vertical divergence of the iso-neutral fluxes (i.e. the term associated
376with a second order vertical derivative) appears only tracer values associated
377with a single water column. This is of paramount importance since it means
378that the implicit in time algorithm for solving the vertical diffusion equation can
379be used to evaluate this term. It is a necessity since the vertical eddy diffusivity
380associated with this term, 
381\begin{equation}
382    \sum_{\substack{i_p, \,k_p}} \left\{ 
383      A_i^k \; \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
384   \right\} 
385\end{equation}
386can be quite large.
387
388\item[$\bullet$ pure iso-neutral operator]  The iso-neutral flux of locally referenced
389potential density is zero, $i.e.$
390\begin{align} \label{Gf_property2}
391\begin{matrix}
392&{_i^k {\mathbb{F}_u}_{i_p}^{k_p} (\rho)} 
393   &=    &\alpha_i^k   &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)
394   &- \ \;  \beta _i^k    &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (S) & = \ 0   \\
395&{_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)} 
396   &=    &\alpha_i^k   &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T)
397   &- \  \; \beta _i^k    &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (S)  &= \ 0
398\end{matrix}
399\end{align}
400This result is trivially obtained using the \eqref{Gf_triads} applied to $T$ and $S$ 
401and the definition of the triads' slopes \eqref{Gf_slopes}.
402
403\item[$\bullet$ conservation of tracer] The iso-neutral diffusion term conserve the
404total tracer content, $i.e.$
405\begin{equation} \label{Gf_property1}
406\sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0
407\end{equation}
408This property is trivially satisfied since the iso-neutral diffusive operator
409is written in flux form.
410
411\item[$\bullet$ decrease of tracer variance] The iso-neutral diffusion term does
412not increase the total tracer variance, $i.e.$
413\begin{equation} \label{Gf_property1}
414\sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0
415\end{equation}
416The property is demonstrated in the Appendix~\ref{Apdx_Gf_operator}. It is a
417key property for a diffusion term. It means that the operator is also a dissipation
418term, $i.e.$ it is a sink term for the square of the quantity on which it is applied.
419It therfore ensure that, when the diffusivity coefficient is large enough, the field
420on which it is applied become free of grid-point noise.
421
422\item[$\bullet$ self-adjoint operator] The iso-neutral diffusion operator is self-adjoint,
423$i.e.$
424\begin{equation} \label{Gf_property1}
425\sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 
426\end{equation}
427In other word, there is no needs to develop a specific routine from the adjoint of this
428operator. We just have to apply the same routine. This properties can be demonstrated
429quite easily in a similar way the "non increase of tracer variance" property has been
430proved (see Appendix~\ref{Apdx_Gf_operator}).
431\end{description}
432
433
434$\ $\newline      %force an empty line
435% ================================================================
436% Skew flux formulation for Eddy Induced Velocity :
437% ================================================================
438\subsection{Eddy induced velocity and skew flux formulation}
439
440When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined),
441an additional advection term is added. The associated velocity is the so called
442eddy induced velocity, the formulation of which depends on the slopes of iso-
443neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used
444here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo} 
445is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo}
446+ \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates.
447
448The eddy induced velocity is given by:
449\begin{equation} \label{Eq_eiv_v}
450\begin{split}
451 u^* & = - \frac{1}{e_2\,e_{3}}          \;\partial_k \left( e_2 \, A_e \; r_\right)   
452          = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_\right)            \\
453 v^* & = - \frac{1}{e_1\,e_3}\;             \partial_k \left( e_1 \, A_e \; r_\right)   
454          = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_\right)             \\
455w^* & =    \frac{1}{e_1\,e_2}\; \left\{   \partial_\left( e_2 \, A_e \; r_\right)
456                             + \partial_\left( e_1 \, A_e \;r_j   \right) \right\}   \\
457\end{split}
458\end{equation}
459where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the
460slopes between the iso-neutral and the geopotential surfaces.
461%%gm wrong: to be modified with 2 2D streamfunctions
462 In other words,
463the eddy induced velocity can be derived from a vector streamfuntion, $\phi$, which
464is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^*  = \textbf{k} \times \nabla \phi$
465%%end gm
466
467A traditional way to implement this additional advection is to add it to the eulerian
468velocity prior to compute the tracer advection. This allows us to take advantage of
469all the advection schemes offered for the tracers (see \S\ref{TRA_adv}) and not just
470a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers
471where \emph{positivity} of the advection scheme is of paramount importance.
472% give here the expression using the triads. It is different from the one given in \eqref{Eq_ldfeiv}
473% see just below a copy of this equation:
474%\begin{equation} \label{Eq_ldfeiv}
475%\begin{split}
476% u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\
477% v^* & = \frac{1}{e_{1u}e_{3v}}\; \delta_k \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right]\\
478%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\} \\
479%\end{split}
480%\end{equation}
481\begin{equation} \label{Eq_eiv_vd} 
482\textbf{F}_{eiv}^T   \equiv   \left( \begin{aligned}                               
483 \sum_{\substack{i_p,\,k_p}} &
484 +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k} 
485\ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}]      \\
486    \\
487 \sum_{\substack{i_p,\,k_p}} &
488 - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p} 
489\ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]    \\   
490\end{aligned}   \right)
491\end{equation}
492
493\ref{Griffies_JPO98} introduces another way to implement the eddy induced advection,
494the so-called skew form. It is based on a transformation of the advective fluxes
495using the non-divergent nature of the eddy induced velocity.
496For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be
497transformed as follows:
498\begin{flalign*}
499\begin{split}
500\textbf{F}_{eiv}^T =
501\begin{pmatrix} 
502           {e_{2}\,e_{3}\;  u^*}       \\
503      {e_{1}\,e_{2}\; w^*}  \\
504\end{pmatrix}   \;   T
505&=
506\begin{pmatrix} 
507           { - \partial_k \left( e_{2} \, A_{e} \; r_i \right) \; T \;}       \\
508      {+ \partial_\left( e_{2} \, A_{e} \; r_i \right) \; T \;}    \\
509\end{pmatrix}        \\
510&=       
511\begin{pmatrix} 
512           { - \partial_k \left( e_{2} \, A_{e} \; r_\; T \right) \;}  \\
513      {+ \partial_\left( e_{2} \, A_{e} \; r_\; T \right) \;}   \\
514\end{pmatrix}       
515 +
516\begin{pmatrix} 
517           {+ e_{2} \, A_{e} \; r_\; \partial_k T}  \\
518      { - e_{2} \, A_{e} \; r_\; \partial_i  T}  \\
519\end{pmatrix}   
520\end{split}
521\end{flalign*}
522and since the eddy induces velocity field is no-divergent, we end up with the skew
523form of the eddy induced advective fluxes:
524\begin{equation} \label{Eq_eiv_skew_continuous}
525\textbf{F}_{eiv}^T = \begin{pmatrix} 
526           {+ e_{2} \, A_{e} \; r_\; \partial_k T}   \\
527      { - e_{2} \, A_{e} \; r_\; \partial_i  T}  \\
528                                 \end{pmatrix}
529\end{equation}
530The tendency associated with eddy induced velocity is then simply the divergence
531of the \eqref{Eq_eiv_skew_continuous} fluxes. It naturally conserves the tracer
532content, as it is expressed in flux form and, as the advective form, it preserve the
533tracer variance. Another interesting property of \eqref{Eq_eiv_skew_continuous} 
534form is that when $A=A_e$, a simplification occurs in the sum of the iso-neutral
535diffusion and eddy induced velocity terms:
536\begin{flalign} \label{Eq_eiv_skew+eiv_continuous}
537\textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &=
538\begin{pmatrix} 
539           + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T -  e_2 \, A \; r_i                              \;\partial_k T   \\
540      -  e_2 \, A_{e} \; r_i           \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T   \\
541\end{pmatrix}
542+
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&= \begin{pmatrix} 
548           + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T    \\
549      -  2\; e_2 \, A_{e} \; r_i      \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T   \\
550\end{pmatrix}
551\end{flalign}
552The horizontal component reduces to the one use for an horizontal laplacian
553operator and the vertical one keep the same complexity, but not more. This property
554has been used to reduce the computational time \citep{Griffies_JPO98}, but it is
555not of practical use as usually $A \neq A_e$. Nevertheless this property can be used to
556choose a discret form of  \eqref{Eq_eiv_skew_continuous} which is consistent with the
557iso-neutral operator \eqref{Gf_operator}. Using the slopes \eqref{Gf_slopes} 
558and defining $A_e$ at $T$-point($i.e.$ as $A$, the eddy diffusivity coefficient),
559the resulting discret form is given by:
560\begin{equation} \label{Eq_eiv_skew} 
561\textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left( \begin{aligned}                               
562 \sum_{\substack{i_p,\,k_p}} &
563 +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k} 
564\ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}]      \\
565    \\
566 \sum_{\substack{i_p,\,k_p}} &
567 - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p} 
568\ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]    \\   
569\end{aligned}   \right)
570\end{equation}
571Note that \eqref{Eq_eiv_skew} is valid in $z$-coordinate with or without partial cells.
572In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces
573must be added to $\mathbb{R}$ for the discret form to be exact.
574
575Such a choice of discretisation is consistent with the iso-neutral operator as it uses the
576same definition for the slopes. It also ensures the conservation of the tracer variance
577(see Appendix \ref{Apdx_eiv_skew}), $i.e.$ it does not include a diffusive component
578but is a "pure" advection term.
579
580
581
582
583$\ $\newpage      %force an empty line
584% ================================================================
585% Discrete Invariants of the iso-neutral diffrusion
586% ================================================================
587\subsection{Discrete invariants of the iso-neutral diffrusion}
588\label{Apdx_Gf_operator}
589
590Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane.
591
592This part will be moved in an Appendix.
593
594The continuous property to be demonstrated is :
595\begin{align*}
596\int_D  D_l^T \; T \;dv   \leq 0
597\end{align*}
598The discrete form of its left hand side is obtained using \eqref{Eq_iso_flux}
599
600\begin{align*}
601&\int_D  D_l^T \; T \;dv \equiv  \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\}    \\
602&\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ 
603      \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] 
604        + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]  \ T \right\}    \\
605&\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ 
606                {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T]
607             + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\
608&\equiv -\sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ 
609   \frac{ _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} (T) }{ {e_{1u}}_{\,i+1/2}^{\,k} }  \ \delta_{i+1/2} [T]
610 - { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \; 
611   \frac{ _i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} (T) }{ {e_{3w}}_{\,i}^{\,k+1/2}  } \ \delta_{k+1/2} [T]   
612                                                                       \right\}      \\
613%
614\allowdisplaybreaks
615\intertext{ Expending the summation on $i_p$ and $k_p$, it becomes:}
616%
617&\equiv -\sum_{i,k}
618\begin{Bmatrix} 
619&\ \ \Bigl{ _{i+1}^{k} \mathbb{T}_{-1/2}^{-1/2} (T) } 
620&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
621& -\ \ {_{i}^{k+1} \mathbb{R}_{-1/2}^{-1/2}} 
622&      {_{i}^{k+1} \mathbb{T}_{-1/2}^{-1/2} (T) }   
623&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
624& \\
625&+\Bigl( \ \;\; { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) } 
626&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
627& -\ \ {_i^{k+1} \mathbb{R}_{+1/2}^{-1/2}}
628& { _i^{k+1} \mathbb{T}_{+1/2}^{-1/2} (T) }   
629&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}      \Bigr)
630& \\
631&+\Bigl{ _{i+1}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) } 
632&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
633& -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{-1/2}^{+1/2}} 
634&      \ \;\;{_{i}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }   
635&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
636& \\
637&+\Bigl( \ \;\; { _{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) } 
638&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
639& -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{+1/2}^{+1/2}} 
640&      \ \;\;{_{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }   
641&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
642\end{Bmatrix}
643%
644\allowdisplaybreaks
645\intertext{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: }
646%
647&\equiv -\sum_{i,k}
648\begin{Bmatrix} 
649&\ \ \Bigl{_i^k \mathbb{T}_{-1/2}^{-1/2} (T) } 
650&\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} 
651& -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}} 
652&      {_i^k \mathbb{T}_{-1/2}^{-1/2} (T) }   
653&\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}     \Bigr)
654& \\
655&+\Bigl{ _i^k \mathbb{T}_{+1/2}^{-1/2} (T) } 
656&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
657& -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}}
658&      { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }   
659&\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}      \Bigr)
660& \\
661&+\Bigl{_i^k \mathbb{T}_{-1/2}^{+1/2} (T) } 
662&\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} 
663& -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}} 
664&      {_i^k \mathbb{T}_{-1/2}^{+1/2} (T) }   
665&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
666& \\
667&+\Bigl( { _i^k \mathbb{T}_{+1/2}^{+1/2} (T) } 
668&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
669& -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}} 
670&      {_i^k \mathbb{T}_{+1/2}^{+1/2} (T) }   
671&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
672\end{Bmatrix}   \\
673%
674\allowdisplaybreaks
675\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 \eqref{Gf_triads}. 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$, then the previous demonstration would have let to:
714\begin{align*}
715\int_D  S \; D_l^\;dv &\equiv  \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\}    \\
716&\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ 
717\left( \frac{ \delta_{i +i_p} [S] }{{e_{1u} }_{\,i+i_p}^{\,k}} 
718 - {_i^k \mathbb{R}_{i_p}^{k_p}} 
719\frac{ \delta_{k+k_p} [S] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right\right.   
720\\   & \qquad \qquad \qquad \ \left.
721\left( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}} 
722 - {_i^k \mathbb{R}_{i_p}^{k_p}} 
723\frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right)
724 \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \
725 \right\}   
726%
727\allowdisplaybreaks
728\intertext{which, by applying the same operation as before but in reverse order, leads to: }
729%
730&\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\}   
731\end{align*} 
732This means that the iso-neutral operator is self-adjoint. There is no need to develop a specific to obtain it.
733
734
735
736$\ $\newpage      %force an empty line
737% ================================================================
738% Discrete Invariants of the skew flux formulation
739% ================================================================
740\subsection{Discrete invariants of the skew flux formulation}
741\label{Apdx_eiv_skew}
742
743
744Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.
745
746This have to be moved in an Appendix.
747
748The continuous property to be demonstrated is :
749\begin{align*}
750\int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0
751\end{align*}
752The discrete form of its left hand side is obtained using \eqref{Eq_eiv_skew}
753\begin{align*}
754 \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
755 \delta_&\left[                                                   
756{e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k} 
757\ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]         
758   \right] \; T_i^k      \\
759- \delta_k &\left[
760{e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2} 
761\ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]   
762   \right] \; T_i^k      \         \Biggr\}   
763\end{align*}
764apply the adjoint of delta operator, it becomes
765\begin{align*}
766 \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
767  &\left(                                                   
768{e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k} 
769\ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]         
770   \right) \; \delta_{i+1/2}[T^{k}]      \\
771- &\left(
772{e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2} 
773\ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]   
774     \right) \; \delta_{k+1/2}[T_{i}]       \         \Biggr\}       
775\end{align*}
776Expending the summation on $i_p$ and $k_p$, it becomes:
777\begin{align*}
778 \begin{matrix} 
779&\sum\limits_{i,k}   \Bigl\{ 
780  &+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k} 
781  &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}]    &\delta_{i+1/2}[T^{k}]   &\\
782&&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}     
783  &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}}  &\delta_{k-1/2}[T_{i\ \ \ \;}&\delta_{i+1/2}[T^{k}] &\\
784&&+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k} 
785  &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}]     &\delta_{i+1/2}[T^{k}] &\\
786&&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}       
787    &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\
788%
789&&-{e_{2u}}_i^{k+1}                                &{A_{e}}_i^{k+1} 
790  &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}}   &\delta_{i-1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\   
791&&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:} 
792  &{\ \ \;_i^\mathbb{R}_{-1/2}^{+1/2}}   &\delta_{i-1/2}[T^{k\ \ \ \:}&\delta_{k+1/2}[T_{i}] &\\
793&&-{e_{2u}}_i^{k+1    }                             &{A_{e}}_i^{k+1} 
794  &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}}   &\delta_{i+1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\   
795&&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:} 
796  &{\ \ \;_i^\mathbb{R}_{+1/2}^{+1/2}}   &\delta_{i+1/2}[T^{k\ \ \ \:}&\delta_{k+1/2}[T_{i}]
797&\Bigr\}  \\
798\end{matrix}   
799\end{align*}
800The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the
801same but of opposite signs, they cancel out.
802Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$.
803The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the
804same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$ 
805they cancel out with the neighbouring grid points.
806Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the
807$i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the
808tracer is preserved by the discretisation of the skew fluxes.
809
810\end{document}
Note: See TracBrowser for help on using the repository browser.