New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Annex_E.tex in branches/UKMO/r5518_INGV1_WAVE-coupling/DOC/TexFiles/Chapters – NEMO

source: branches/UKMO/r5518_INGV1_WAVE-coupling/DOC/TexFiles/Chapters/Annex_E.tex @ 7152

Last change on this file since 7152 was 7152, checked in by jcastill, 7 years ago

Initial implementation of wave coupling branch - INGV wave branch + UKMO wave coupling branch

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