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.
Changeset 2282 for branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Annex_E.tex – NEMO

Ignore:
Timestamp:
2010-10-15T16:42:00+02:00 (14 years ago)
Author:
gm
Message:

ticket:#658 merge DOC of all the branches that form the v3.3 beta

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Annex_E.tex

    r1223 r2282  
    11% ================================================================ 
    2 % Appendix E : UBS scheme 
    3 % ================================================================ 
    4 \chapter{UBS scheme} 
    5 \label{AN_E} 
     2% Appendix E : Note on some algorithms 
     3% ================================================================ 
     4\chapter{Note on some algorithms} 
     5\label{Apdx_E} 
    66\minitoc 
    77 
     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 
    815 
    916% ------------------------------------------------------------------------------------------------------------- 
    … …  
    3946 
    4047This results in a dissipatively dominant (i.e. hyper-diffusive) truncation  
    41 error \citep{Sacha2005}. The overall performance of the  
     48error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of the  
    4249advection scheme is similar to that reported in \cite{Farrow1995}.  
    4350It is a relatively good compromise between accuracy and smoothness. It is  
    … …  
    5663(centred in time) while the second term which is the diffusive part of the scheme,  
    5764is evaluated using the \textit{before} velocity (forward in time. This is discussed  
    58 by \citet{Webb1998} in the context of the Quick advection scheme. UBS and QUICK  
     65by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. UBS and QUICK  
    5966schemes only differ by one coefficient. Substituting 1/6 with 1/8 in  
    60 (\ref{Eq_tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb1998}.  
     67(\ref{Eq_tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.  
    6168This option is not available through a namelist parameter, since the 1/6  
    6269coefficient is hard coded. Nevertheless it is quite easy to make the  
    … …  
    7481evaluated using either \textit{(a)} a centred $2^{nd}$ order scheme ,  
    7582or  \textit{(b)} a TVD scheme, or  \textit{(c)} an interpolation based on conservative  
    76 parabolic splines following \citet{Sacha2005} implementation of UBS in ROMS,  
     83parabolic splines following \citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS,  
    7784or  \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an  
    7885eight-order accurate conventional scheme. 
    … …  
    185192\begin{subequations} \label{dt_mt} 
    186193\begin{align} 
    187  \delta _{t+\Delta t/2} [q]     &=  \  \ \,   q^{t+\Delta t}  - q^{t}      \\ 
    188  \overline q^{\,t+\Delta t/2} &= \left\{ q^{t+\Delta t} + q^{t} \right\} \; / \; 2 
     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 
    189196\end{align} 
    190197\end{subequations} 
    191198As for space operator, the adjoint of the derivation and averaging time operators are  
    192 $\delta_t^*=\delta_{t+\Delta t/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$ 
     199$\delta_t^*=\delta_{t+\rdt/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$ 
    193200, respectively.  
    194201 
    … …  
    196203\begin{equation} \label{LF} 
    197204   \frac{\partial q}{\partial t}  
    198          \equiv \frac{1}{\Delta t} \overline{ \delta _{t+\Delta t/2}[q]}^{\,t}  
    199       =         \frac{q^{t+\Delta t}-q^{t-\Delta t}}{2\Delta t} 
     205         \equiv \frac{1}{\rdt} \overline{ \delta _{t+\rdt/2}[q]}^{\,t}  
     206      =         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt} 
    200207\end{equation}  
    201 Note that \eqref{LF} shows that the leapfrog time step is $\Delta t$, not $2\Delta t$  
     208Note that \eqref{LF} shows that the leapfrog time step is $\rdt$, not $2\rdt$  
    202209as it can be found sometime in literature.  
    203210The leap-Frog time stepping is a second order centered scheme. As such it respects  
    … …  
    212219\int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}  
    213220   &\equiv \sum\limits_{0}^{N}  
    214          {\frac{1}{\Delta t} q^t \ \overline{ \delta _{t+\Delta t/2}[q]}^{\,t} \ \Delta t}  
    215       \equiv \sum\limits_{0}^{N}  { q^t \ \overline{ \delta _{t+\Delta t/2}[q]}^{\,t} } \\ 
    216    &\equiv \sum\limits_{0}^{N}  { \overline{q}^{\,t+\Delta/2}{ \delta _{t+\Delta t/2}[q]}} 
    217       \equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta _{t+\Delta t/2}[q^2] }\\ 
    218    &\equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta _{t+\Delta t/2}[q^2] } 
     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] } 
    219226      \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) 
    220227\end{split} \end{equation} 
    … …  
    226233 
    227234 
     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^k  \ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] } 
     286{\left(\alpha / \beta \right)_i^k  \ \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{Triads used in the Griffies's like iso-neutral diffision scheme for  
     303$u$-component (upper panel) and $w$-component (lower panel).} 
     304\end{center} 
     305\end{figure} 
     306%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     307 
     308The four iso-neutral fluxes associated with the triads are defined at $T$-point.  
     309They take the following expression : 
     310\begin{flalign} \label{Gf_fluxes} 
     311\begin{split} 
     312{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)  
     313   &= \ \; \qquad  \quad    { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+i_p}^{\,k}}    \\ 
     314{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T) 
     315   &=  -\; { _i^k \mathbb{R}_{i_p}^{k_p} } 
     316             \ \; { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+k_p}} 
     317\end{split} 
     318\end{flalign} 
     319 
     320The resulting iso-neutral fluxes at $u$- and $w$-points are then given by the  
     321sum of the fluxes that cross the $u$- and $w$-face (Fig.~\ref{Fig_ISO_triad}): 
     322\begin{flalign} \label{Eq_iso_flux}  
     323\textbf{F}_{iso}(T)  
     324&\equiv  \sum_{\substack{i_p,\,k_p}}  
     325   \begin{pmatrix}  
     326      {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\ 
     327      \\ 
     328      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)      \\    
     329   \end{pmatrix}    \notag \\ 
     330&  \notag \\ 
     331&\equiv  \sum_{\substack{i_p,\,k_p}}  
     332   \begin{pmatrix}  
     333      && { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+1/2}^{\,k} }    \\ 
     334      \\ 
     335      & -\; { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } 
     336        & {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+1/2} }   \\    
     337   \end{pmatrix}      % \\ 
     338% &\\ 
     339% &\equiv  \sum_{\substack{i_p,\,k_p}}  
     340%    \begin{pmatrix}  
     341%       \qquad  \qquad  \qquad  
     342%       \frac{1}{ {e_{1u}}_{\,i+1/2}^{\,k} }  \ \;   
     343%        { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T)\\ 
     344%       \\ 
     345%       -\frac{1}{ {e_{3w}}_{\,i}^{\,k+1/2} }  \ \;  
     346%        { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \;   
     347%                  {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T)\\    
     348%    \end{pmatrix}       
     349\end{flalign} 
     350resulting in a iso-neutral diffusion tendency on temperature given by the divergence  
     351of the sum of all the four triad fluxes : 
     352\begin{equation} \label{Gf_operator} 
     353D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{   
     354       \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]  
     355        + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]   \right\} 
     356\end{equation} 
     357where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells.  
     358 
     359This expression of the iso-neutral diffusion has been chosen in order to satisfy  
     360the following six properties: 
     361\begin{description} 
     362\item[$\bullet$ horizontal diffusion] The discretization of the diffusion operator  
     363recovers the traditional five-point Laplacian in the limit of flat iso-neutral direction : 
     364\begin{equation} \label{Gf_property1a} 
     365D_l^T = \frac{1}{b_T}  \ \delta_{i}  
     366   \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right]  
     367\qquad  \text{when} \quad  
     368   { _i^k \mathbb{R}_{i_p}^{k_p} }=0 
     369\end{equation} 
     370 
     371\item[$\bullet$ implicit treatment in the vertical]  In the diagonal term associated  
     372with the vertical divergence of the iso-neutral fluxes (i.e. the term associated  
     373with a second order vertical derivative) appears only tracer values associated  
     374with a single water column. This is of paramount importance since it means 
     375that the implicit in time algorithm for solving the vertical diffusion equation can  
     376be used to evaluate this term. It is a necessity since the vertical eddy diffusivity  
     377associated with this term,   
     378\begin{equation} 
     379    \sum_{\substack{i_p, \,k_p}} \left\{   
     380      A_i^k \; \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 
     381   \right\}  
     382\end{equation} 
     383can be quite large. 
     384 
     385\item[$\bullet$ pure iso-neutral operator]  The iso-neutral flux of locally referenced  
     386potential density is zero, $i.e.$ 
     387\begin{align} \label{Gf_property2} 
     388\begin{matrix} 
     389&{_i^k {\mathbb{F}_u}_{i_p}^{k_p} (\rho)}  
     390   &=    &\alpha_i^k   &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)  
     391   &- \ \;  \beta _i^k    &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (S) & = \ 0   \\ 
     392&{_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)}  
     393   &=    &\alpha_i^k   &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T)  
     394   &- \  \; \beta _i^k    &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (S)  &= \ 0 
     395\end{matrix} 
     396\end{align} 
     397This result is trivially obtained using the \eqref{Gf_triads} applied to $T$ and $S$  
     398and the definition of the triads' slopes \eqref{Gf_slopes}. 
     399 
     400\item[$\bullet$ conservation of tracer] The iso-neutral diffusion term conserve the  
     401total tracer content, $i.e.$ 
     402\begin{equation} \label{Gf_property1} 
     403\sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0 
     404\end{equation} 
     405This property is trivially satisfied since the iso-neutral diffusive operator  
     406is written in flux form. 
     407 
     408\item[$\bullet$ decrease of tracer variance] The iso-neutral diffusion term does  
     409not increase the total tracer variance, $i.e.$ 
     410\begin{equation} \label{Gf_property1} 
     411\sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0 
     412\end{equation} 
     413The property is demonstrated in the Appendix~\ref{Apdx_Gf_operator}. It is a  
     414key property for a diffusion term. It means that the operator is also a dissipation  
     415term, $i.e.$ it is a sink term for the square of the quantity on which it is applied.  
     416It therfore ensure that, when the diffusivity coefficient is large enough, the field  
     417on which it is applied become free of grid-point noise. 
     418 
     419\item[$\bullet$ self-adjoint operator] The iso-neutral diffusion operator is self-adjoint,  
     420$i.e.$ 
     421\begin{equation} \label{Gf_property1} 
     422\sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}  
     423\end{equation} 
     424In other word, there is no needs to develop a specific routine from the adjoint of this  
     425operator. We just have to apply the same routine. This properties can be demonstrated  
     426quite easily in a similar way the "non increase of tracer variance" property has been  
     427proved (see Appendix~\ref{Apdx_Gf_operator}). 
     428\end{description} 
     429 
     430 
     431$\ $\newline      %force an empty line 
     432% ================================================================ 
     433% Skew flux formulation for Eddy Induced Velocity :  
     434% ================================================================ 
     435\subsection{ Eddy induced velocity and Skew flux formulation} 
     436 
     437When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined),  
     438an additional advection term is added. The associated velocity is the so called  
     439eddy induced velocity, the formulation of which depends on the slopes of iso- 
     440neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used  
     441here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo}  
     442is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo} 
     443+ \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates.  
     444 
     445The eddy induced velocity is given by:  
     446\begin{equation} \label{Eq_eiv_v} 
     447\begin{split} 
     448 u^* & = - \frac{1}{e_2\,e_{3}}          \;\partial_k \left( e_2 \, A_e \; r_i  \right)    
     449          = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_i  \right)            \\ 
     450 v^* & = - \frac{1}{e_1\,e_3}\;             \partial_k \left( e_1 \, A_e \; r_j  \right)    
     451          = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_j  \right)             \\ 
     452w^* & =    \frac{1}{e_1\,e_2}\; \left\{   \partial_i  \left( e_2 \, A_e \; r_i  \right)  
     453                             + \partial_j  \left( e_1 \, A_e \;r_j   \right) \right\}   \\ 
     454\end{split} 
     455\end{equation} 
     456where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the  
     457slopes between the iso-neutral and the geopotential surfaces. 
     458%%gm wrong: to be modified with 2 2D streamfunctions 
     459 In other words, 
     460the eddy induced velocity can be derived from a vector streamfuntion, $\phi$, which 
     461is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^*  = \textbf{k} \times \nabla \phi$ 
     462%%end gm 
     463 
     464A traditional way to implement this additional advection is to add it to the eulerian  
     465velocity prior to compute the tracer advection. This allows us to take advantage of  
     466all the advection schemes offered for the tracers (see \S\ref{TRA_adv}) and not just  
     467a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers  
     468where \emph{positivity} of the advection scheme is of paramount importance.  
     469% give here the expression using the triads. It is different from the one given in \eqref{Eq_ldfeiv} 
     470% see just below a copy of this equation: 
     471%\begin{equation} \label{Eq_ldfeiv} 
     472%\begin{split} 
     473% u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\ 
     474% v^* & = \frac{1}{e_{1u}e_{3v}}\; \delta_k \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right]\\ 
     475%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\} \\ 
     476%\end{split} 
     477%\end{equation} 
     478\begin{equation} \label{Eq_eiv_vd}   
     479\textbf{F}_{eiv}^T   \equiv   \left( \begin{aligned}                                 
     480 \sum_{\substack{i_p,\,k_p}} & 
     481 +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k}  
     482\ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}]      \\ 
     483    \\ 
     484 \sum_{\substack{i_p,\,k_p}} & 
     485 - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p}  
     486\ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]    \\    
     487\end{aligned}   \right) 
     488\end{equation} 
     489 
     490\ref{Griffies_JPO98} introduces another way to implement the eddy induced advection,  
     491the so-called skew form. It is based on a transformation of the advective fluxes  
     492using the non-divergent nature of the eddy induced velocity.  
     493For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be  
     494transformed as follows: 
     495\begin{flalign*} 
     496\begin{split} 
     497\textbf{F}_{eiv}^T =  
     498\begin{pmatrix}  
     499           {e_{2}\,e_{3}\;  u^*}       \\ 
     500      {e_{1}\,e_{2}\; w^*}  \\ 
     501\end{pmatrix}   \;   T 
     502&= 
     503\begin{pmatrix}  
     504           { - \partial_k \left( e_{2} \, A_{e} \; r_i \right) \; T \;}       \\ 
     505      {+ \partial_i  \left( e_{2} \, A_{e} \; r_i \right) \; T \;}    \\ 
     506\end{pmatrix}        \\ 
     507&=        
     508\begin{pmatrix}  
     509           { - \partial_k \left( e_{2} \, A_{e} \; r_i  \; T \right) \;}  \\ 
     510      {+ \partial_i  \left( e_{2} \, A_{e} \; r_i  \; T \right) \;}   \\ 
     511\end{pmatrix}         
     512 +  
     513\begin{pmatrix}  
     514           {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}  \\ 
     515      { - e_{2} \, A_{e} \; r_i  \; \partial_i  T}  \\ 
     516\end{pmatrix}    
     517\end{split} 
     518\end{flalign*} 
     519and since the eddy induces velocity field is no-divergent, we end up with the skew  
     520form of the eddy induced advective fluxes: 
     521\begin{equation} \label{Eq_eiv_skew_continuous} 
     522\textbf{F}_{eiv}^T = \begin{pmatrix}  
     523           {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}   \\ 
     524      { - e_{2} \, A_{e} \; r_i  \; \partial_i  T}  \\ 
     525                                 \end{pmatrix} 
     526\end{equation} 
     527The tendency associated with eddy induced velocity is then simply the divergence  
     528of the \eqref{Eq_eiv_skew_continuous} fluxes. It naturally conserves the tracer  
     529content, as it is expressed in flux form and, as the advective form, it preserve the  
     530tracer variance. Another interesting property of \eqref{Eq_eiv_skew_continuous}  
     531form is that when $A=A_e$, a simplification occurs in the sum of the iso-neutral  
     532diffusion and eddy induced velocity terms: 
     533\begin{flalign} \label{Eq_eiv_skew+eiv_continuous} 
     534\textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &=  
     535\begin{pmatrix}  
     536           + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T -  e_2 \, A \; r_i                              \;\partial_k T   \\ 
     537      -  e_2 \, A_{e} \; r_i           \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T   \\ 
     538\end{pmatrix} 
     539+ 
     540\begin{pmatrix}  
     541           {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}   \\ 
     542      { - e_{2} \, A_{e} \; r_i  \; \partial_i  T}  \\ 
     543\end{pmatrix}      \\ 
     544&= \begin{pmatrix}  
     545           + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T    \\ 
     546      -  2\; e_2 \, A_{e} \; r_i      \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T   \\ 
     547\end{pmatrix} 
     548\end{flalign} 
     549The horizontal component reduces to the one use for an horizontal laplacian  
     550operator and the vertical one keep the same complexity, but not more. This property 
     551has been used to reduce the computational time \citep{Griffies_JPO98}, but it is  
     552not of practical use as usually $A \neq A_e$. Nevertheless this property can be used to  
     553choose a discret form of  \eqref{Eq_eiv_skew_continuous} which is consistent with the  
     554iso-neutral operator \eqref{Gf_operator}. Using the slopes \eqref{Gf_slopes}  
     555and defining $A_e$ at $T$-point($i.e.$ as $A$, the eddy diffusivity coefficient), 
     556the resulting discret form is given by: 
     557\begin{equation} \label{Eq_eiv_skew}   
     558\textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left( \begin{aligned}                                 
     559 \sum_{\substack{i_p,\,k_p}} & 
     560 +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k}  
     561\ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}]      \\ 
     562    \\ 
     563 \sum_{\substack{i_p,\,k_p}} & 
     564 - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p}  
     565\ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]    \\    
     566\end{aligned}   \right) 
     567\end{equation} 
     568Note that \eqref{Eq_eiv_skew} is valid in $z$-coordinate with or without partial cells.  
     569In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces  
     570must be added to $\mathbb{R}$ for the discret form to be exact.  
     571 
     572Such a choice of discretisation is consistent with the iso-neutral operator as it uses the  
     573same definition for the slopes. It also ensures the conservation of the tracer variance  
     574(see Appendix \ref{Apdx_eiv_skew}), $i.e.$ it does not include a diffusive component  
     575but is a "pure" advection term. 
     576 
     577 
     578 
     579 
     580$\ $\newpage      %force an empty line 
     581% ================================================================ 
     582% Discrete Invariants of the iso-neutral diffrusion  
     583% ================================================================ 
     584\subsection{Discrete Invariants of the iso-neutral diffrusion} 
     585\label{Apdx_Gf_operator} 
     586 
     587Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane.  
     588 
     589This part will be moved in an Appendix. 
     590 
     591The continuous property to be demonstrated is : 
     592\begin{align*} 
     593\int_D  D_l^T \; T \;dv   \leq 0 
     594\end{align*} 
     595The discrete form of its left hand side is obtained using \eqref{Eq_iso_flux} 
     596 
     597\begin{align*} 
     598&\int_D  D_l^T \; T \;dv \equiv  \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\}    \\ 
     599&\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{   
     600      \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]  
     601        + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]  \ T \right\}    \\ 
     602&\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{   
     603                {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T] 
     604             + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\ 
     605&\equiv -\sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{   
     606   \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] 
     607 - { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \;  
     608   \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]    
     609                                                                       \right\}      \\ 
     610% 
     611\allowdisplaybreaks 
     612\intertext{ Expending the summation on $i_p$ and $k_p$, it becomes:} 
     613% 
     614&\equiv -\sum_{i,k} 
     615\begin{Bmatrix}   
     616&\ \ \Bigl(  { _{i+1}^{k} \mathbb{T}_{-1/2}^{-1/2} (T) }  
     617&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}  
     618& -\ \ {_{i}^{k+1} \mathbb{R}_{-1/2}^{-1/2}}  
     619&      {_{i}^{k+1} \mathbb{T}_{-1/2}^{-1/2} (T) }    
     620&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr) 
     621& \\ 
     622&+\Bigl( \ \;\; { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }   
     623&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}  
     624& -\ \ {_i^{k+1} \mathbb{R}_{+1/2}^{-1/2}} 
     625& { _i^{k+1} \mathbb{T}_{+1/2}^{-1/2} (T) }    
     626&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}      \Bigr) 
     627& \\ 
     628&+\Bigl(  { _{i+1}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }  
     629&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}  
     630& -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{-1/2}^{+1/2}}  
     631&      \ \;\;{_{i}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }    
     632&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr) 
     633& \\ 
     634&+\Bigl( \ \;\; { _{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }  
     635&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}  
     636& -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{+1/2}^{+1/2}}  
     637&      \ \;\;{_{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }    
     638&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\ 
     639\end{Bmatrix} 
     640% 
     641\allowdisplaybreaks 
     642\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: } 
     643% 
     644&\equiv -\sum_{i,k} 
     645\begin{Bmatrix}   
     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 \mathbb{R}_{-1/2}^{-1/2}}  
     649&      {_i^k \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^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& \\ 
     664&+\Bigl( { _i^k \mathbb{T}_{+1/2}^{+1/2} (T) }  
     665&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}  
     666& -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}}  
     667&      {_i^k \mathbb{T}_{+1/2}^{+1/2} (T) }    
     668&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\ 
     669\end{Bmatrix}   \\ 
     670% 
     671\allowdisplaybreaks 
     672\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: } 
     673% 
     674&\equiv -\sum_{i,k} 
     675\begin{Bmatrix}   
     676&\ \ \Bigl(  \frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}  
     677& -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}}  
     678&\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}     \Bigr)^2 
     679& \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k}  \  A_i^k 
     680& \\ 
     681&+\Bigl(  \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}  
     682& -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}} 
     683&\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}      \Bigr)^2 
     684& \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k}  \  A_i^k 
     685& \\ 
     686&+\Bigl(  \frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}  
     687& -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}}  
     688&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)^2 
     689& \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k}  \  A_i^k 
     690& \\ 
     691&+\Bigl( \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}  
     692& -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}}  
     693&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)^2 
     694& \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k}  \  A_i^k      \\ 
     695\end{Bmatrix}   \\ 
     696& \\ 
     697% 
     698&\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{   
     699\begin{matrix}   
     700&\Bigl( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}}  
     701& -\ \ {_i^k \mathbb{R}_{i_p}^{k_p}}  
     702&\frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \Bigr)^2 
     703& \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \ \  
     704\end{matrix} 
     705 \right\}    
     706\quad   \leq 0 
     707\end{align*}  
     708The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities. 
     709 
     710Note 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: 
     711\begin{align*} 
     712\int_D  S \; D_l^T  \;dv &\equiv  \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\}    \\ 
     713&\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{   
     714\left( \frac{ \delta_{i +i_p} [S] }{{e_{1u} }_{\,i+i_p}^{\,k}}  
     715 - {_i^k \mathbb{R}_{i_p}^{k_p}}  
     716\frac{ \delta_{k+k_p} [S] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right)  \right.     
     717\\   & \qquad \qquad \qquad \ \left. 
     718\left( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}}  
     719 - {_i^k \mathbb{R}_{i_p}^{k_p}}  
     720\frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right)  
     721 \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \ 
     722 \right\}    
     723% 
     724\allowdisplaybreaks 
     725\intertext{which, by applying the same operation as before but in reverse order, leads to: } 
     726% 
     727&\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\}    
     728\end{align*}  
     729This means that the iso-neutral operator is self-adjoint. There is no need to develop a specific to obtain it. 
     730 
     731 
     732 
     733$\ $\newpage      %force an empty line 
     734% ================================================================ 
     735% Discrete Invariants of the skew flux formulation 
     736% ================================================================ 
     737\subsection{Discrete Invariants of the skew flux formulation} 
     738\label{Apdx_eiv_skew} 
     739 
     740 
     741Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.  
     742 
     743This have to be moved in an Appendix. 
     744 
     745The continuous property to be demonstrated is : 
     746\begin{align*} 
     747\int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0 
     748\end{align*} 
     749The discrete form of its left hand side is obtained using \eqref{Eq_eiv_skew} 
     750\begin{align*} 
     751 \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\; 
     752 \delta_i  &\left[                                                     
     753{e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}  
     754\ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]          
     755   \right] \; T_i^k      \\ 
     756- \delta_k &\left[  
     757{e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}  
     758\ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]    
     759   \right] \; T_i^k      \         \Biggr\}    
     760\end{align*} 
     761apply the adjoint of delta operator, it becomes 
     762\begin{align*} 
     763 \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\; 
     764  &\left(                                                     
     765{e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}  
     766\ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]          
     767   \right) \; \delta_{i+1/2}[T^{k}]      \\ 
     768- &\left(  
     769{e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}  
     770\ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]    
     771     \right) \; \delta_{k+1/2}[T_{i}]       \         \Biggr\}        
     772\end{align*} 
     773Expending the summation on $i_p$ and $k_p$, it becomes: 
     774\begin{align*} 
     775 \begin{matrix}   
     776&\sum\limits_{i,k}   \Bigl\{  
     777  &+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}  
     778  &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}]    &\delta_{i+1/2}[T^{k}]   &\\ 
     779&&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}       
     780  &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}}  &\delta_{k-1/2}[T_{i\ \ \ \;}]  &\delta_{i+1/2}[T^{k}] &\\ 
     781&&+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}  
     782  &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}]     &\delta_{i+1/2}[T^{k}] &\\ 
     783&&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}        
     784    &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ 
     785% 
     786&&-{e_{2u}}_i^{k+1}                                &{A_{e}}_i^{k+1}  
     787  &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}}   &\delta_{i-1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\    
     788&&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}  
     789  &{\ \ \;_i^k  \mathbb{R}_{-1/2}^{+1/2}}   &\delta_{i-1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}] &\\ 
     790&&-{e_{2u}}_i^{k+1    }                             &{A_{e}}_i^{k+1}  
     791  &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}}   &\delta_{i+1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\    
     792&&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}  
     793  &{\ \ \;_i^k  \mathbb{R}_{+1/2}^{+1/2}}   &\delta_{i+1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}]  
     794&\Bigr\}  \\ 
     795\end{matrix}    
     796\end{align*} 
     797The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the  
     798same but of opposite signs, they cancel out.  
     799Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$.  
     800The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the  
     801same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$  
     802they cancel out with the neighbouring grid points.  
     803Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the  
     804$i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the  
     805tracer is preserved by the discretisation of the skew fluxes. 
     806 
Note: See TracChangeset for help on using the changeset viewer.