Index: NEMO/trunk/doc/latex/NEMO/subfiles/annex_E.tex
===================================================================
--- NEMO/trunk/doc/latex/NEMO/subfiles/annex_E.tex	(revision 10406)
+++ NEMO/trunk/doc/latex/NEMO/subfiles/annex_E.tex	(revision 10414)
@@ -1,3 +1,4 @@
-\documentclass[../tex_main/NEMO_manual]{subfiles}
+\documentclass[../main/NEMO_manual]{subfiles}
+
 \begin{document}
 % ================================================================
@@ -6,12 +7,10 @@
 \chapter{Note on some algorithms}
 \label{apdx:E}
+
 \minitoc
 
 \newpage
-$\ $\newline    % force a new ligne
 
 This appendix some on going consideration on algorithms used or planned to be used in \NEMO. 
-
-$\ $\newline    % force a new ligne
 
 % -------------------------------------------------------------------------------------------------------------
@@ -25,15 +24,19 @@
 It is also known as Cell Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective Kinematics).
 For example, in the $i$-direction:
-\begin{equation} \label{eq:tra_adv_ubs2}
-\tau_u^{ubs} = \left\{	 \begin{aligned}
-  & \tau_u^{cen4} + \frac{1}{12} \,\tau"_i	   & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
-  & \tau_u^{cen4} - \frac{1}{12} \,\tau"_{i+1} & \quad \text{if }\ u_{i+1/2}       <       0
-  						 \end{aligned}    \right.
+\begin{equation}
+  \label{eq:tra_adv_ubs2}
+  \tau_u^{ubs} = \left\{
+	 \begin{aligned}
+      & \tau_u^{cen4} + \frac{1}{12} \,\tau"_i	   & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
+      & \tau_u^{cen4} - \frac{1}{12} \,\tau"_{i+1} & \quad \text{if }\ u_{i+1/2}       <       0
+    \end{aligned}
+  \right.
 \end{equation}
 or equivalently, the advective flux is
-\begin{equation} \label{eq:tra_adv_ubs2}
-U_{i+1/2} \ \tau_u^{ubs} 
-=U_{i+1/2} \ \overline{ T_i - \frac{1}{6}\,\tau"_i }^{\,i+1/2}
-- \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
+\begin{equation}
+  \label{eq:tra_adv_ubs2}
+  U_{i+1/2} \ \tau_u^{ubs}
+  =U_{i+1/2} \ \overline{ T_i - \frac{1}{6}\,\tau"_i }^{\,i+1/2}
+  - \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
 \end{equation}
 where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and
@@ -44,5 +47,4 @@
 Alternative choice: introduce the scale factors:  
 $\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]$.
-
 
 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error
@@ -82,17 +84,21 @@
 
 NB 3: It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows:
-\begin{equation} \label{eq:tra_adv_ubs2}
-\tau_u^{ubs} = \left\{	 \begin{aligned}
-	& \tau_u^{cen4} + \frac{1}{12} \tau"_i		& \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
-	& \tau_u^{cen4} - \frac{1}{12} \tau"_{i+1}	& \quad \text{if }\ u_{i+1/2}       <       0
-  						 \end{aligned}    \right.
+\begin{equation}
+  \label{eq:tra_adv_ubs2}
+  \tau_u^{ubs} = \left\{
+	 \begin{aligned}
+      & \tau_u^{cen4} + \frac{1}{12} \tau"_i		& \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
+      & \tau_u^{cen4} - \frac{1}{12} \tau"_{i+1}	& \quad \text{if }\ u_{i+1/2}       <       0
+    \end{aligned}
+  \right.
 \end{equation}
 or equivalently 
-\begin{equation} \label{eq:tra_adv_ubs2}
-\begin{split}
-e_{2u} e_{3u}\,u_{i+1/2} \ \tau_u^{ubs} 
-&= e_{2u} e_{3u}\,u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\tau"_i }^{\,i+1/2} \\
-& - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
-\end{split}
+\begin{equation}
+  \label{eq:tra_adv_ubs2}
+  \begin{split}
+    e_{2u} e_{3u}\,u_{i+1/2} \ \tau_u^{ubs}
+    &= e_{2u} e_{3u}\,u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\tau"_i }^{\,i+1/2} \\
+    & - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
+  \end{split}
 \end{equation}
 \autoref{eq:tra_adv_ubs2} has several advantages.
@@ -105,55 +111,58 @@
 
 laplacian diffusion:
-\begin{equation} \label{eq:tra_ldf_lap}
-\begin{split}
-D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\;  e_{3T} } &\left[ {\quad \delta_i 
-\left[ {A_u^{lT} \frac{e_{2u} e_{3u} }{e_{1u} }\;\delta_{i+1/2} 
-\left[ T \right]} \right]} \right.
-\\
-&\ \left. {+\; \delta_j \left[ 
-{A_v^{lT} \left( {\frac{e_{1v} e_{3v} }{e_{2v} }\;\delta_{j+1/2} \left[ T 
-\right]} \right)} \right]\quad } \right]
-\end{split}
+\begin{equation}
+  \label{eq:tra_ldf_lap}
+  \begin{split}
+    D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\;  e_{3T} } &\left[ {\quad \delta_i
+        \left[ {A_u^{lT} \frac{e_{2u} e_{3u} }{e_{1u} }\;\delta_{i+1/2}
+            \left[ T \right]} \right]} \right. \\
+    &\ \left. {+\; \delta_j \left[
+          {A_v^{lT} \left( {\frac{e_{1v} e_{3v} }{e_{2v} }\;\delta_{j+1/2} \left[ T
+                \right]} \right)} \right]\quad } \right]
+  \end{split}
 \end{equation}
 
 bilaplacian:
-\begin{equation} \label{eq:tra_ldf_lap}
-\begin{split}
-D_T^{lT} =&-\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
-& \delta_i \left[  \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta_{i+1/2} 
-		  \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}}
-    \delta_i \left[ \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta_{i+1/2} 
-    		  [T] \right] \right] \right]
-\end{split}
+\begin{equation}
+  \label{eq:tra_ldf_lap}
+  \begin{split}
+    D_T^{lT} =&-\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
+    & \delta_i \left[  \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta_{i+1/2}
+      \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}}
+        \delta_i \left[ \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta_{i+1/2}
+          [T] \right] \right] \right]
+  \end{split}
 \end{equation}
 with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$, 
 $i.e.$ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$
 it comes:
-\begin{equation} \label{eq:tra_ldf_lap}
-\begin{split}
-D_T^{lT} =&-\frac{1}{12}\,\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
-& \delta_i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta_{i+1/2} 
-		 \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} 
-    \delta_i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta_{i+1/2} 
-    		[T] \right] \right] \right]
-\end{split}
+\begin{equation}
+  \label{eq:tra_ldf_lap}
+  \begin{split}
+    D_T^{lT} =&-\frac{1}{12}\,\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
+    & \delta_i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta_{i+1/2}
+      \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}}
+        \delta_i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta_{i+1/2}
+          [T] \right] \right] \right]
+  \end{split}
 \end{equation}
 if the velocity is uniform ($i.e.$ $|u|=cst$) then the diffusive flux is
-\begin{equation} \label{eq:tra_ldf_lap}
-\begin{split}
-F_u^{lT} = - \frac{1}{12}
- e_{2u}\,e_{3u}\,|u| \;\sqrt{ e_{1u}}\,\delta_{i+1/2} 
-		 \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} 
-    \delta_i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}}\:\delta_{i+1/2} 
-    		[T] \right] \right]
-\end{split}
+\begin{equation}
+  \label{eq:tra_ldf_lap}
+  \begin{split}
+    F_u^{lT} = - \frac{1}{12}
+    e_{2u}\,e_{3u}\,|u| \;\sqrt{ e_{1u}}\,\delta_{i+1/2}
+    \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}}
+      \delta_i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}}\:\delta_{i+1/2}
+        [T] \right] \right]
+  \end{split}
 \end{equation}
 beurk....  reverte the logic: starting from the diffusive part of the advective flux it comes:
 
-\begin{equation} \label{eq:tra_adv_ubs2}
-\begin{split}
-F_u^{lT}
-&= - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
-\end{split}
+\begin{equation}
+  \label{eq:tra_adv_ubs2}
+  \begin{split}
+    F_u^{lT} &= - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
+  \end{split}
 \end{equation}
 if the velocity is uniform ($i.e.$ $|u|=cst$) and
@@ -161,22 +170,21 @@
 
 sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$):
-\begin{equation} \label{eq:tra_adv_ubs2}
-\begin{split}
-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]
-\end{split}
+\begin{equation}
+  \label{eq:tra_adv_ubs2}
+  \begin{split}
+    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]
+  \end{split}
 \end{equation}
 which leads to ${A_T^{lT}}^2 = \frac{1}{12} {e_{1T}}^3\ \overline{|u|}^{\,i+1/2}$
 
 sol 2 coefficient at u-point: split $|u|$ into $\sqrt{|u|}$ and $e_{1T}$ into $\sqrt{e_{1u}}$
-\begin{equation} \label{eq:tra_adv_ubs2}
-\begin{split}
-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] \\
-&= - \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]
-\end{split}
+\begin{equation}
+  \label{eq:tra_adv_ubs2}
+  \begin{split}
+    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] \\
+    &= - \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]
+  \end{split}
 \end{equation}
 which leads to ${A_u^{lT}} = \frac{1}{12} {e_{1u}}^3\ |u|$
-
 
 % -------------------------------------------------------------------------------------------------------------
@@ -189,10 +197,11 @@
 Given the values of a variable $q$ at successive time step,
 the time derivation and averaging operators at the mid time step are:
-\begin{subequations} \label{eq:dt_mt}
-\begin{align}
- \delta_{t+\rdt/2} [q]     &=  \  \ \,   q^{t+\rdt}  - q^{t}		\\
- \overline q^{\,t+\rdt/2} &= \left\{ q^{t+\rdt} + q^{t} \right\} \; / \; 2
-\end{align}
-\end{subequations}
+\[
+  % \label{eq:dt_mt}
+  \begin{split}
+    \delta_{t+\rdt/2} [q]     &=  \  \ \,   q^{t+\rdt}  - q^{t}		\\
+    \overline q^{\,t+\rdt/2} &= \left\{ q^{t+\rdt} + q^{t} \right\} \; / \; 2
+  \end{split}
+\]
 As for space operator,
 the adjoint of the derivation and averaging time operators are $\delta_t^*=\delta_{t+\rdt/2}$ and
@@ -200,39 +209,38 @@
 
 The Leap-frog time stepping given by \autoref{eq:DOM_nxt} can be defined as:
-\begin{equation} \label{eq:LF}
-   \frac{\partial q}{\partial t} 
-   		\equiv \frac{1}{\rdt} \overline{ \delta_{t+\rdt/2}[q]}^{\,t} 
-		=         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt}
-\end{equation} 
+\[
+  % \label{eq:LF}
+  \frac{\partial q}{\partial t}
+  \equiv \frac{1}{\rdt} \overline{ \delta_{t+\rdt/2}[q]}^{\,t}
+  =         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt}
+\] 
 Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$,
 not $2\rdt$ as it can be found sometimes in literature.
 The leap-Frog time stepping is a second order centered scheme.
 As such it respects the quadratic invariant in integral forms, $i.e.$ the following continuous property,
-\begin{equation} \label{eq:Energy}
-\int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} 
-	=\int_{t_0}^{t_1} {\frac{1}{2}\, \frac{\partial q^2}{\partial t} \;dt} 
-	=  \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) ,
-\end{equation}
+\[
+  % \label{eq:Energy}
+  \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}
+  =\int_{t_0}^{t_1} {\frac{1}{2}\, \frac{\partial q^2}{\partial t} \;dt}
+  =  \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) ,
+\]
 is satisfied in discrete form.
 Indeed, 
-\begin{equation} \begin{split}
-\int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} 
-	&\equiv \sum\limits_{0}^{N} 
-			{\frac{1}{\rdt} q^t \ \overline{ \delta_{t+\rdt/2}[q]}^{\,t} \ \rdt} 
-	   \equiv \sum\limits_{0}^{N}  { q^t \ \overline{ \delta_{t+\rdt/2}[q]}^{\,t} } \\
-	&\equiv \sum\limits_{0}^{N}  { \overline{q}^{\,t+\Delta/2}{ \delta_{t+\rdt/2}[q]}}
-	   \equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta_{t+\rdt/2}[q^2] }\\
-	&\equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta_{t+\rdt/2}[q^2] }
-	   \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right)
-\end{split} \end{equation}
+\[
+  \begin{split}
+    \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}
+    &\equiv \sum\limits_{0}^{N}
+    {\frac{1}{\rdt} q^t \ \overline{ \delta_{t+\rdt/2}[q]}^{\,t} \ \rdt}
+    \equiv \sum\limits_{0}^{N}  { q^t \ \overline{ \delta_{t+\rdt/2}[q]}^{\,t} } \\
+    &\equiv \sum\limits_{0}^{N}  { \overline{q}^{\,t+\Delta/2}{ \delta_{t+\rdt/2}[q]}}
+    \equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta_{t+\rdt/2}[q^2] }\\
+    &\equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta_{t+\rdt/2}[q^2] }
+    \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right)
+  \end{split}
+\]
 NB here pb of boundary condition when applying the adjoint!
 In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition 
 (equivalently of the boundary value of the integration by part).
 In time this boundary condition is not physical and \textbf{add something here!!!}
-
-
-
-
-
 
 % ================================================================
@@ -269,10 +277,11 @@
 a derivative in the same direction by considering triads.
 For example in the (\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows:
-\begin{equation} \label{eq:Gf_triads}
-_i^k \mathbb{T}_{i_p}^{k_p} (T)
-= \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k  	\left(  
-                                                     \frac{ \delta_{i + i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
--\ {_i^k \mathbb{R}_{i_p}^{k_p}} \ \frac{ \delta_{k+k_p} [T^i] }{ {e_{3w}}_{\,i}^{\,k+k_p} } 
-									  \right)
+\begin{equation}
+  \label{eq:Gf_triads}
+  _i^k \mathbb{T}_{i_p}^{k_p} (T)
+  = \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k  	\left(
+    \frac{ \delta_{i + i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
+    -\ {_i^k \mathbb{R}_{i_p}^{k_p}} \ \frac{ \delta_{k+k_p} [T^i] }{ {e_{3w}}_{\,i}^{\,k+k_p} }
+  \right)
 \end{equation}
 where the indices $i_p$ and $k_p$ define the four triads and take the following value:
@@ -281,9 +290,10 @@
 $A_i^k$ is the lateral eddy diffusivity coefficient defined at $T$-point,
 and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad:
-\begin{equation} \label{eq:Gf_slopes}
-_i^k \mathbb{R}_{i_p}^{k_p} 
-=\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} \ \frac 
-{\left(\alpha / \beta \right)_i^k  \ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] }
-{\left(\alpha / \beta \right)_i^k  \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }
+\begin{equation}
+  \label{eq:Gf_slopes}
+  _i^k \mathbb{R}_{i_p}^{k_p}
+  =\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} \ \frac
+  {\left(\alpha / \beta \right)_i^k  \ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] }
+  {\left(\alpha / \beta \right)_i^k  \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }
 \end{equation}
 Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of
@@ -296,10 +306,13 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht] \begin{center}
-\includegraphics[width=0.70\textwidth]{Fig_ISO_triad}
-\caption{  \protect\label{fig:ISO_triad}
-  Triads used in the Griffies's like iso-neutral diffision scheme for
-  $u$-component (upper panel) and $w$-component (lower panel).}
-\end{center}
+\begin{figure}[!ht]
+  \begin{center}
+    \includegraphics[width=0.70\textwidth]{Fig_ISO_triad}
+    \caption{
+      \protect\label{fig:ISO_triad}
+      Triads used in the Griffies's like iso-neutral diffision scheme for
+      $u$-component (upper panel) and $w$-component (lower panel).
+    }
+  \end{center}
 \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -307,50 +320,52 @@
 The four iso-neutral fluxes associated with the triads are defined at $T$-point. 
 They take the following expression:
-\begin{flalign} \label{eq:Gf_fluxes}
-\begin{split}
-{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) 
-   &= \ \; \qquad  \quad    { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+i_p}^{\,k}}    \\
-{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T)
-   &=  -\; { _i^k \mathbb{R}_{i_p}^{k_p} }
-				 \ \; { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+k_p}}
-\end{split}
-\end{flalign}
+\begin{flalign*}
+  % \label{eq:Gf_fluxes}
+  \begin{split}
+    {_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)
+    &= \ \; \qquad  \quad    { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+i_p}^{\,k}}    \\
+    {_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T)
+    &=  -\; { _i^k \mathbb{R}_{i_p}^{k_p} }
+    \ \; { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+k_p}}
+  \end{split}
+\end{flalign*}
 
 The resulting iso-neutral fluxes at $u$- and $w$-points are then given by
 the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}):
-\begin{flalign} \label{eq:iso_flux} 
-\textbf{F}_{iso}(T) 
-&\equiv  \sum_{\substack{i_p,\,k_p}} 
-   \begin{pmatrix} 
-      {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\
-      \\
-      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)      \\   
-   \end{pmatrix}    \notag \\
-&  \notag \\
-&\equiv  \sum_{\substack{i_p,\,k_p}} 
-   \begin{pmatrix} 
-      && { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+1/2}^{\,k} }    \\
-      \\
-      & -\; { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }
-        & {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+1/2} }   \\   
-   \end{pmatrix}      % \\
-% &\\
-% &\equiv  \sum_{\substack{i_p,\,k_p}} 
-%    \begin{pmatrix} 
-%       \qquad  \qquad  \qquad 
-%       \frac{1}{ {e_{1u}}_{\,i+1/2}^{\,k} }  \ \;  
-%        { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T)\\
-%       \\
-%       -\frac{1}{ {e_{3w}}_{\,i}^{\,k+1/2} }  \ \; 
-%        { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \;  
-%                  {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T)\\   
-%    \end{pmatrix}      
+\begin{flalign}
+  \label{eq:iso_flux}
+  \textbf{F}_{iso}(T)
+  &\equiv  \sum_{\substack{i_p,\,k_p}}
+  \begin{pmatrix}
+    {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) \\ \\
+    {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)
+  \end{pmatrix}
+  \notag \\
+  &  \notag \\
+  &\equiv  \sum_{\substack{i_p,\,k_p}}
+  \begin{pmatrix}
+    && { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+1/2}^{\,k} } \\ \\
+    & -\; { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }
+    & {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+1/2} }
+  \end{pmatrix}      % \\
+  % &\\
+  % &\equiv  \sum_{\substack{i_p,\,k_p}}
+  % \begin{pmatrix}
+  %   \qquad  \qquad  \qquad
+  %   \frac{1}{ {e_{1u}}_{\,i+1/2}^{\,k} }  \ \;
+  %   { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T)\\
+  %   \\
+  %   -\frac{1}{ {e_{3w}}_{\,i}^{\,k+1/2} }  \ \;
+  %   { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \;
+  %   {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T)\\
+  % \end{pmatrix}
 \end{flalign}
 resulting in a iso-neutral diffusion tendency on temperature given by
 the divergence of the sum of all the four triad fluxes:
-\begin{equation} \label{eq:Gf_operator}
-D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{  
-		 \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] 
-	     + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]   \right\}
+\begin{equation}
+  \label{eq:Gf_operator}
+  D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{
+    \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
+    + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]   \right\}
 \end{equation}
 where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. 
@@ -361,10 +376,11 @@
   The discretization of the diffusion operator recovers the traditional five-point Laplacian in
   the limit of flat iso-neutral direction:
-\begin{equation} \label{eq:Gf_property1a}
-D_l^T = \frac{1}{b_T}  \ \delta_{i} 
-	\left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right] 
-\qquad  \text{when} \quad 
-	{ _i^k \mathbb{R}_{i_p}^{k_p} }=0
-\end{equation}
+  \[
+    % \label{eq:Gf_property1a}
+    D_l^T = \frac{1}{b_T}  \ \delta_{i}
+    \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right]
+    \qquad  \text{when} \quad
+    { _i^k \mathbb{R}_{i_p}^{k_p} }=0
+  \]
 
 \item[$\bullet$ implicit treatment in the vertical]
@@ -374,39 +390,42 @@
   This is of paramount importance since it means that
   the implicit in time algorithm for solving the vertical diffusion equation can be used to evaluate this term.
-  It is a necessity since the vertical eddy diffusivity associated with this term,  
-\begin{equation}
-	 \sum_{\substack{i_p, \,k_p}} \left\{  
+  It is a necessity since the vertical eddy diffusivity associated with this term,
+  \[
+	 \sum_{\substack{i_p, \,k_p}} \left\{
 		A_i^k \; \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
-	\right\} 
-\end{equation}
-can be quite large.
+    \right\}
+  \]
+  can be quite large.
 
 \item[$\bullet$ pure iso-neutral operator]
   The iso-neutral flux of locally referenced potential density is zero, $i.e.$
-\begin{align} \label{eq:Gf_property2}
-\begin{matrix}
-&{_i^k {\mathbb{F}_u}_{i_p}^{k_p} (\rho)} 
-	&=    &\alpha_i^k   &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) 
-	&- \ \;  \beta _i^k    &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (S) & = \ 0   \\
-&{_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)} 
-	&=    &\alpha_i^k   &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T) 
-	&- \  \; \beta _i^k    &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (S)  &= \ 0
-\end{matrix}
-\end{align}
-This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$ and
-the definition of the triads' slopes \autoref{eq:Gf_slopes}.
+  \begin{align*}
+    % \label{eq:Gf_property2}
+    \begin{matrix}
+      &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} (\rho)}
+      &=    &\alpha_i^k   &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)
+      &- \ \;  \beta _i^k    &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (S) & = \ 0   \\
+      &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)}
+      &=    &\alpha_i^k   &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T)
+      &- \  \; \beta _i^k    &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (S)  &= \ 0
+    \end{matrix}
+  \end{align*}
+  This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$ and
+  the definition of the triads' slopes \autoref{eq:Gf_slopes}.
 
 \item[$\bullet$ conservation of tracer]
   The iso-neutral diffusion term conserve the total tracer content, $i.e.$
-\begin{equation} \label{eq:Gf_property1}
-\sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0
-\end{equation}
+  \[
+    % \label{eq:Gf_property1}
+    \sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0
+  \]
 This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form.
 
 \item[$\bullet$ decrease of tracer variance]
   The iso-neutral diffusion term does not increase the total tracer variance, $i.e.$
-\begin{equation} \label{eq:Gf_property1}
-\sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0
-\end{equation}
+  \[
+    % \label{eq:Gf_property1}
+    \sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0
+  \]
 The property is demonstrated in the \autoref{apdx:Gf_operator}.
 It is a key property for a diffusion term.
@@ -418,7 +437,8 @@
 \item[$\bullet$ self-adjoint operator]
   The iso-neutral diffusion operator is self-adjoint, $i.e.$
-\begin{equation} \label{eq:Gf_property1}
-\sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 
-\end{equation}
+  \[
+    % \label{eq:Gf_property1}
+    \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}
+  \]
 In other word, there is no needs to develop a specific routine from the adjoint of this operator.
 We just have to apply the same routine.
@@ -427,6 +447,4 @@
 \end{description}
 
-
-$\ $\newline      %force an empty line
 % ================================================================
 % Skew flux formulation for Eddy Induced Velocity : 
@@ -443,13 +461,14 @@
 
 The eddy induced velocity is given by: 
-\begin{equation} \label{eq:eiv_v}
-\begin{split}
- u^* & = - \frac{1}{e_2\,e_{3}}          \;\partial_k \left( e_2 \, A_e \; r_i  \right)   
-          = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_i  \right)            \\
- v^* & = - \frac{1}{e_1\,e_3}\;             \partial_k \left( e_1 \, A_e \; r_j  \right)   
-          = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_j  \right)             \\
-w^* & =    \frac{1}{e_1\,e_2}\; \left\{   \partial_i  \left( e_2 \, A_e \; r_i  \right) 
-					              + \partial_j  \left( e_1 \, A_e \;r_j   \right) \right\}   \\
-\end{split}
+\begin{equation}
+  \label{eq:eiv_v}
+  \begin{split}
+    u^* & = - \frac{1}{e_2\,e_{3}}          \;\partial_k \left( e_2 \, A_e \; r_i  \right)
+    = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_i  \right)            \\
+    v^* & = - \frac{1}{e_1\,e_3}\;             \partial_k \left( e_1 \, A_e \; r_j  \right)
+    = - \frac{1}{e_3}                     \;\partial_k \left(           A_e \; r_j  \right)             \\
+    w^* & =    \frac{1}{e_1\,e_2}\; \left\{   \partial_i  \left( e_2 \, A_e \; r_i  \right)
+      + \partial_j  \left( e_1 \, A_e \;r_j   \right) \right\}
+  \end{split}
 \end{equation}
 where $A_{e}$ is the eddy induced velocity coefficient,
@@ -475,15 +494,17 @@
 %\end{split}
 %\end{equation}
-\begin{equation} \label{eq:eiv_vd}  
-\textbf{F}_{eiv}^T   \equiv   \left( \begin{aligned}                                
- \sum_{\substack{i_p,\,k_p}} &
- +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k} 
-\ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}]      \\
-    \\
- \sum_{\substack{i_p,\,k_p}} &
- - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p} 
-\ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]    \\   
-\end{aligned}   \right)
-\end{equation}
+\[
+  % \label{eq:eiv_vd}
+  \textbf{F}_{eiv}^T   \equiv   \left(
+    \begin{aligned}
+      \sum_{\substack{i_p,\,k_p}} &
+      +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k}
+      \ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ \\
+      \sum_{\substack{i_p,\,k_p}} &
+      - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p}
+      \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]
+    \end{aligned}
+  \right)
+\]
 
 \citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form.
@@ -491,34 +512,38 @@
 For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows:
 \begin{flalign*}
-\begin{split}
-\textbf{F}_{eiv}^T = 
-\begin{pmatrix} 
- 	        {e_{2}\,e_{3}\;  u^*} 	 	\\
- 		{e_{1}\,e_{2}\; w^*}	 \\
-\end{pmatrix}   \;   T
-&=
-\begin{pmatrix} 
- 	        { - \partial_k \left( e_{2} \, A_{e} \; r_i \right) \; T \;} 	 	\\
- 		{+ \partial_i  \left( e_{2} \, A_{e} \; r_i \right) \; T \;}	 \\
-\end{pmatrix} 			\\
-&=			
-\begin{pmatrix} 
- 	        { - \partial_k \left( e_{2} \, A_{e} \; r_i  \; T \right) \;}  \\
- 		{+ \partial_i  \left( e_{2} \, A_{e} \; r_i  \; T \right) \;}	 \\
-\end{pmatrix} 			
- + 
-\begin{pmatrix} 
- 	        {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}  \\
- 		{ - e_{2} \, A_{e} \; r_i  \; \partial_i  T}	 \\
-\end{pmatrix} 	 
-\end{split}
+  \begin{split}
+    \textbf{F}_{eiv}^T =
+    \begin{pmatrix}
+      {e_{2}\,e_{3}\;  u^*} 	 	\\
+ 		{e_{1}\,e_{2}\; w^*}
+    \end{pmatrix}
+    \;   T
+    &=
+    \begin{pmatrix}
+      { - \partial_k \left( e_{2} \, A_{e} \; r_i \right) \; T \;} 	 	\\
+ 		{+ \partial_i  \left( e_{2} \, A_{e} \; r_i \right) \; T \;}
+    \end{pmatrix}
+    \\
+    &=
+    \begin{pmatrix}
+      { - \partial_k \left( e_{2} \, A_{e} \; r_i  \; T \right) \;}  \\
+ 		{+ \partial_i  \left( e_{2} \, A_{e} \; r_i  \; T \right) \;}
+    \end{pmatrix}
+    +
+    \begin{pmatrix}
+      {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}  \\
+ 		{ - e_{2} \, A_{e} \; r_i  \; \partial_i  T}
+    \end{pmatrix}
+  \end{split}
 \end{flalign*}
 and since the eddy induces velocity field is no-divergent,
 we end up with the skew form of the eddy induced advective fluxes:
-\begin{equation} \label{eq:eiv_skew_continuous}
-\textbf{F}_{eiv}^T = \begin{pmatrix} 
- 	        {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}   \\
- 		{ - e_{2} \, A_{e} \; r_i  \; \partial_i  T}	 \\
-                                 \end{pmatrix}
+\begin{equation}
+  \label{eq:eiv_skew_continuous}
+  \textbf{F}_{eiv}^T =
+  \begin{pmatrix}
+    {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}   \\
+    { - e_{2} \, A_{e} \; r_i  \; \partial_i  T}
+  \end{pmatrix}
 \end{equation}
 The tendency associated with eddy induced velocity is then simply the divergence of
@@ -528,20 +553,23 @@
 Another interesting property of \autoref{eq:eiv_skew_continuous} form is that when $A=A_e$,
 a simplification occurs in the sum of the iso-neutral diffusion and eddy induced velocity terms:
-\begin{flalign} \label{eq:eiv_skew+eiv_continuous}
-\textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &= 
-\begin{pmatrix} 
- 	        + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T -  e_2 \, A \; r_i                              \;\partial_k T   \\
- 		-  e_2 \, A_{e} \; r_i           \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T   \\
-\end{pmatrix}
-+
-\begin{pmatrix} 
- 	        {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}   \\
- 		{ - e_{2} \, A_{e} \; r_i  \; \partial_i  T}	 \\
-\end{pmatrix}      \\
-&= \begin{pmatrix} 
- 	        + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T    \\
- 		-  2\; e_2 \, A_{e} \; r_i      \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T   \\
-\end{pmatrix}
-\end{flalign}
+\begin{flalign*}
+  % \label{eq:eiv_skew+eiv_continuous}
+  \textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &=
+  \begin{pmatrix}
+    + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T -  e_2 \, A \; r_i                              \;\partial_k T   \\
+    -  e_2 \, A_{e} \; r_i           \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T
+  \end{pmatrix}
+  +
+  \begin{pmatrix}
+    {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}   \\
+    { - e_{2} \, A_{e} \; r_i  \; \partial_i  T}
+  \end{pmatrix}
+  \\
+  &=
+  \begin{pmatrix}
+    + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T    \\
+    -  2\; e_2 \, A_{e} \; r_i      \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T
+  \end{pmatrix}
+\end{flalign*}
 The horizontal component reduces to the one use for an horizontal laplacian operator and
 the vertical one keeps the same complexity, but not more.
@@ -552,14 +580,16 @@
 Using the slopes \autoref{eq:Gf_slopes} and defining $A_e$ at $T$-point($i.e.$ as $A$,
 the eddy diffusivity coefficient), the resulting discret form is given by:
-\begin{equation} \label{eq:eiv_skew}  
-\textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left( \begin{aligned}                                
- \sum_{\substack{i_p,\,k_p}} &
- +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k} 
-\ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}]      \\
-    \\
- \sum_{\substack{i_p,\,k_p}} &
- - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p} 
-\ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]    \\   
-\end{aligned}   \right)
+\begin{equation}
+  \label{eq:eiv_skew}
+  \textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left(
+    \begin{aligned}
+      \sum_{\substack{i_p,\,k_p}} &
+      +{e_{2u}}_{i+1/2-i_p}^{k}                                  \ \ {A_{e}}_{i+1/2-i_p}^{k}
+      \ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ \\
+      \sum_{\substack{i_p,\,k_p}} &
+      - {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p}
+      \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]
+    \end{aligned}
+  \right)
 \end{equation}
 Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells. 
@@ -572,7 +602,4 @@
 $i.e.$ it does not include a diffusive component but is a "pure" advection term.
 
-
-
-
 $\ $\newpage      %force an empty line
 % ================================================================
@@ -587,125 +614,130 @@
 
 The continuous property to be demonstrated is:
+\[
+  \int_D  D_l^T \; T \;dv   \leq 0
+\]
+The discrete form of its left hand side is obtained using \autoref{eq:iso_flux}
+
 \begin{align*}
-\int_D  D_l^T \; T \;dv   \leq 0
-\end{align*}
-The discrete form of its left hand side is obtained using \autoref{eq:iso_flux}
-
-\begin{align*}
-&\int_D  D_l^T \; T \;dv \equiv  \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\}    \\
-&\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{  
-		\delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] 
-	     + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]  \ T \right\}    \\
-&\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{  
-                {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T]
-             + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\
-&\equiv -\sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{  
-	\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]
- - { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \; 
-	\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]   
-	                                                                    \right\}      \\
-%
-\allowdisplaybreaks
-\intertext{ Expending the summation on $i_p$ and $k_p$, it becomes:}
-%
-&\equiv -\sum_{i,k}
-\begin{Bmatrix}  
-&\ \ \Bigl(  { _{i+1}^{k} \mathbb{T}_{-1/2}^{-1/2} (T) } 
-&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
-& -\ \ {_{i}^{k+1} \mathbb{R}_{-1/2}^{-1/2}} 
-&      {_{i}^{k+1} \mathbb{T}_{-1/2}^{-1/2} (T) }   
-&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
-& \\
-&+\Bigl( \ \;\; { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }  
-&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
-& -\ \ {_i^{k+1} \mathbb{R}_{+1/2}^{-1/2}}
-& { _i^{k+1} \mathbb{T}_{+1/2}^{-1/2} (T) }   
-&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}      \Bigr)
-& \\
-&+\Bigl(  { _{i+1}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) } 
-&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
-& -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{-1/2}^{+1/2}} 
-&      \ \;\;{_{i}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }   
-&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
-& \\
-&+\Bigl( \ \;\; { _{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) } 
-&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
-& -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{+1/2}^{+1/2}} 
-&      \ \;\;{_{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }   
-&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
-\end{Bmatrix}
-%
-\allowdisplaybreaks
-  \intertext{The summation is done over all $i$ and $k$ indices,
+  &\int_D  D_l^T \; T \;dv \equiv  \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\}    \\
+  &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
+    \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
+    + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]  \ T \right\}    \\
+  &\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
+    {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T]
+    + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\
+  &\equiv -\sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
+    \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]
+    - { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \;
+    \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]
+    \right\}      \\
+    %
+  \allowdisplaybreaks
+  \intertext{ Expending the summation on $i_p$ and $k_p$, it becomes:}
+  %
+  &\equiv -\sum_{i,k}
+    \begin{Bmatrix}
+      &\ \ \Bigl(  { _{i+1}^{k} \mathbb{T}_{-1/2}^{-1/2} (T) }
+      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
+      & -\ \ {_{i}^{k+1} \mathbb{R}_{-1/2}^{-1/2}}
+      &      {_{i}^{k+1} \mathbb{T}_{-1/2}^{-1/2} (T) }
+      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
+      & \\
+      &+\Bigl( \ \;\; { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }
+      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
+      & -\ \ {_i^{k+1} \mathbb{R}_{+1/2}^{-1/2}}
+      & { _i^{k+1} \mathbb{T}_{+1/2}^{-1/2} (T) }
+      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}      \Bigr)
+      & \\
+      &+\Bigl(  { _{i+1}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }
+      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
+      & -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{-1/2}^{+1/2}}
+      &      \ \;\;{_{i}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) }
+      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
+      & \\
+      &+\Bigl( \ \;\; { _{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }
+      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
+      & -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{+1/2}^{+1/2}}
+      &      \ \;\;{_{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) }
+      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
+    \end{Bmatrix}
+    %
+  \allowdisplaybreaks
+  \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: }
-%
-&\equiv -\sum_{i,k}
-\begin{Bmatrix}  
-&\ \ \Bigl(  {_i^k \mathbb{T}_{-1/2}^{-1/2} (T) } 
-&\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} 
-& -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}} 
-&      {_i^k \mathbb{T}_{-1/2}^{-1/2} (T) }   
-&\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}     \Bigr)
-& \\
-&+\Bigl(  { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }  
-&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
-& -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}}
-&      { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }   
-&\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}      \Bigr)
-& \\
-&+\Bigl(  {_i^k \mathbb{T}_{-1/2}^{+1/2} (T) } 
-&\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} 
-& -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}} 
-&      {_i^k \mathbb{T}_{-1/2}^{+1/2} (T) }   
-&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
-& \\
-&+\Bigl( { _i^k \mathbb{T}_{+1/2}^{+1/2} (T) } 
-&\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
-& -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}} 
-&      {_i^k \mathbb{T}_{+1/2}^{+1/2} (T) }   
-&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
-\end{Bmatrix}   \\
-%
-\allowdisplaybreaks
-  \intertext{Then outing in factor the triad in each of the four terms of the summation and
+  It becomes:
+  }
+  %
+  &\equiv -\sum_{i,k}
+    \begin{Bmatrix}
+      &\ \ \Bigl(  {_i^k \mathbb{T}_{-1/2}^{-1/2} (T) }
+      &\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
+      & -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}}
+      &      {_i^k \mathbb{T}_{-1/2}^{-1/2} (T) }
+      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}     \Bigr)
+      & \\
+      &+\Bigl(  { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }
+      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
+      & -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}}
+      &      { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) }
+      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}      \Bigr)
+      & \\
+      &+\Bigl(  {_i^k \mathbb{T}_{-1/2}^{+1/2} (T) }
+      &\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
+      & -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}}
+      &      {_i^k \mathbb{T}_{-1/2}^{+1/2} (T) }
+      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)
+      & \\
+      &+\Bigl( { _i^k \mathbb{T}_{+1/2}^{+1/2} (T) }
+      &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
+      & -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}}
+      &      {_i^k \mathbb{T}_{+1/2}^{+1/2} (T) }
+      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)   \\
+    \end{Bmatrix}   \\
+    %
+  \allowdisplaybreaks
+  \intertext{
+  Then outing in factor the triad in each of the four terms of the summation and
   substituting the triads by their expression given in \autoref{eq:Gf_triads}.
-  It becomes: }
-%
-&\equiv -\sum_{i,k}
-\begin{Bmatrix}  
-&\ \ \Bigl(  \frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} 
-& -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}} 
-&\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}     \Bigr)^2
-& \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k}  \  A_i^k
-& \\
-&+\Bigl(  \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
-& -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}}
-&\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}      \Bigr)^2
-& \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k}  \  A_i^k
-& \\
-&+\Bigl(  \frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} 
-& -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}} 
-&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)^2
-& \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k}  \  A_i^k
-& \\
-&+\Bigl( \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} 
-& -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}} 
-&\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)^2
-& \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k}  \  A_i^k      \\
-\end{Bmatrix}   \\
-& \\
-%
-&\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{  
-\begin{matrix}  
-&\Bigl( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}} 
-& -\ \ {_i^k \mathbb{R}_{i_p}^{k_p}} 
-&\frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \Bigr)^2
-& \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \ \ 
-\end{matrix}
- \right\}   
-\quad   \leq 0
+  It becomes:
+  }
+  %
+  &\equiv -\sum_{i,k}
+    \begin{Bmatrix}
+      &\ \ \Bigl(  \frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
+      & -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}}
+      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}     \Bigr)^2
+      & \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k}  \  A_i^k
+      & \\
+      &+\Bigl(  \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
+      & -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}}
+      &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}}      \Bigr)^2
+      & \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k}  \  A_i^k
+      & \\
+      &+\Bigl(  \frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}}
+      & -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}}
+      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)^2
+      & \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k}  \  A_i^k
+      & \\
+      &+\Bigl( \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}}
+      & -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}}
+      &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}}     \Bigr)^2
+      & \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k}  \  A_i^k      \\
+    \end{Bmatrix}
+  \\
+  & \\
+  %
+  &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
+    \begin{matrix}
+      &\Bigl( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}}
+      & -\ \ {_i^k \mathbb{R}_{i_p}^{k_p}}
+      &\frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \Bigr)^2
+      & \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \ \
+    \end{matrix}
+        \right\}
+        \quad   \leq 0
 \end{align*} 
 The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities.
@@ -714,27 +746,28 @@
 then the previous demonstration would have let to:
 \begin{align*}
-\int_D  S \; D_l^T  \;dv &\equiv  \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\}    \\
-&\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{  
-\left( \frac{ \delta_{i +i_p} [S] }{{e_{1u} }_{\,i+i_p}^{\,k}} 
- - {_i^k \mathbb{R}_{i_p}^{k_p}} 
-\frac{ \delta_{k+k_p} [S] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right)  \right.    
-\\   & \qquad \qquad \qquad \ \left.
-\left( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}} 
- - {_i^k \mathbb{R}_{i_p}^{k_p}} 
-\frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right) 
- \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \
- \right\}   
-%
-\allowdisplaybreaks
-\intertext{which, by applying the same operation as before but in reverse order, leads to: }
-%
-&\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\}   
+  \int_D  S \; D_l^T  \;dv &\equiv  \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\}    \\
+                           &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
+                             \left( \frac{ \delta_{i +i_p} [S] }{{e_{1u} }_{\,i+i_p}^{\,k}}
+                             - {_i^k \mathbb{R}_{i_p}^{k_p}}
+                             \frac{ \delta_{k+k_p} [S] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right)  \right. \\
+                           & \qquad \qquad \qquad \ \left.
+                             \left( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}}
+                             - {_i^k \mathbb{R}_{i_p}^{k_p}}
+                             \frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}}     \right)
+                             \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k   \
+                             \right\}
+                             %
+                             \allowdisplaybreaks
+                             \intertext{
+                             which, by applying the same operation as before but in reverse order, leads to:
+                             }
+                             %
+                           &\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\}
 \end{align*} 
 This means that the iso-neutral operator is self-adjoint.
 There is no need to develop a specific to obtain it.
 
-
-
-$\ $\newpage      %force an empty line
+\newpage
+
 % ================================================================
 % Discrete Invariants of the skew flux formulation
@@ -743,5 +776,4 @@
 \label{subsec:eiv_skew}
 
-
 Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane. 
 
@@ -750,53 +782,53 @@
 The continuous property to be demonstrated is:
 \begin{align*}
-\int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0
+  \int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0
 \end{align*}
 The discrete form of its left hand side is obtained using \autoref{eq:eiv_skew}
 \begin{align*}
- \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
- \delta_i  &\left[                                                    
-{e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k} 
-\ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]         
-   \right] \; T_i^k      \\
-- \delta_k &\left[ 
-{e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2} 
-\ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]   
-   \right] \; T_i^k      \         \Biggr\}   
+  \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
+  \delta_i  &\left[
+              {e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}
+              \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]
+              \right] \; T_i^k      \\
+  - \delta_k &\left[
+               {e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}
+               \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]
+               \right] \; T_i^k      \         \Biggr\}
 \end{align*}
 apply the adjoint of delta operator, it becomes
 \begin{align*}
- \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
-  &\left(                                                    
-{e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k} 
-\ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]         
-   \right) \; \delta_{i+1/2}[T^{k}]      \\
-- &\left( 
-{e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2} 
-\ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]   
-     \right) \; \delta_{k+1/2}[T_{i}]       \         \Biggr\}       
+  \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
+  &\left(
+    {e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}
+    \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]
+    \right) \; \delta_{i+1/2}[T^{k}]      \\
+  - &\left(
+      {e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}
+      \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]
+      \right) \; \delta_{k+1/2}[T_{i}]       \         \Biggr\}
 \end{align*}
 Expending the summation on $i_p$ and $k_p$, it becomes:
 \begin{align*}
- \begin{matrix}  
-&\sum\limits_{i,k}   \Bigl\{ 
-  &+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k} 
-  &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}]    &\delta_{i+1/2}[T^{k}]   &\\
-&&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}      
-  &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}}  &\delta_{k-1/2}[T_{i\ \ \ \;}]  &\delta_{i+1/2}[T^{k}] &\\
-&&+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k} 
-  &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}]     &\delta_{i+1/2}[T^{k}] &\\
-&&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}       
+  \begin{matrix}
+    &\sum\limits_{i,k}   \Bigl\{
+    &+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}
+    &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}]    &\delta_{i+1/2}[T^{k}]   &\\
+    &&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}
+    &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}}  &\delta_{k-1/2}[T_{i\ \ \ \;}]  &\delta_{i+1/2}[T^{k}] &\\
+    &&+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}
+    &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}]     &\delta_{i+1/2}[T^{k}] &\\
+    &&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}
     &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\
-%
-&&-{e_{2u}}_i^{k+1}                                &{A_{e}}_i^{k+1} 
-  &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}}   &\delta_{i-1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\   
-&&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:} 
-  &{\ \ \;_i^k  \mathbb{R}_{-1/2}^{+1/2}}   &\delta_{i-1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}] &\\
-&&-{e_{2u}}_i^{k+1    }                             &{A_{e}}_i^{k+1} 
-  &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}}   &\delta_{i+1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\   
-&&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:} 
-  &{\ \ \;_i^k  \mathbb{R}_{+1/2}^{+1/2}}   &\delta_{i+1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}] 
-&\Bigr\}  \\
-\end{matrix}   
+    %
+    &&-{e_{2u}}_i^{k+1}                                &{A_{e}}_i^{k+1}
+    &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}}   &\delta_{i-1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\
+    &&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}
+    &{\ \ \;_i^k  \mathbb{R}_{-1/2}^{+1/2}}   &\delta_{i-1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}] &\\
+    &&-{e_{2u}}_i^{k+1    }                             &{A_{e}}_i^{k+1}
+    &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}}   &\delta_{i+1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\
+    &&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}
+    &{\ \ \;_i^k  \mathbb{R}_{+1/2}^{+1/2}}   &\delta_{i+1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}]
+    &\Bigr\}  \\
+  \end{matrix}   
 \end{align*}
 The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs,
@@ -810,3 +842,5 @@
 $i.e.$ the variance of the tracer is preserved by the discretisation of the skew fluxes.
 
+\biblio
+
 \end{document}
