# source:NEMO/trunk/doc/latex/NEMO/subfiles/annex_E.tex@10442 Last change on this file since 10442 was 10442, checked in by nicolasmartin, 22 months ago

Front page edition, cleaning in custom LaTeX commands and add index for single subfile compilation

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