source: branches/dev_005_AWL/DOC/TexFiles/Chapters/Annex_C.tex @ 3050

% ================================================================
% Chapter Ñ Appendix C : Discrete Invariants of the Equations
% ================================================================
\chapter{Discrete Invariants of the Equations}
\label{Apdx_C}
\minitoc

%%%  Appendix put in gmcomment as it has not been updated for z* and s coordinate
I'm writting this appendix. It will be available in a forthcoming release of the documentation

%\gmcomment{

% ================================================================
% Conservation Properties on Ocean Dynamics
% ================================================================
\section{Conservation Properties on Ocean Dynamics}
\label{Apdx_C.1}

First, the boundary condition on the vertical velocity (no flux through the surface 
and the bottom) is established for the discrete set of momentum equations. 
Then, it is shown that the non-linear terms of the momentum equation are written 
such that the potential enstrophy of a horizontally non-divergent flow is preserved 
while all the other non-diffusive terms preserve the kinetic energy; in practice the 
energy is also preserved. In addition, an option is also offered for the vorticity term 
discretization which provides a total kinetic energy conserving discretization for 
that term. 

Nota Bene: these properties are established here in the rigid-lid case and for the 
2nd order centered scheme. A forthcoming update will be their generalisation to 
the free surface case and higher order scheme.

% -------------------------------------------------------------------------------------------------------------
%       Bottom Boundary Condition on Vertical Velocity Field
% -------------------------------------------------------------------------------------------------------------
\subsection{Bottom Boundary Condition on Vertical Velocity Field}
\label{Apdx_C.1.1}


The discrete set of momentum equations used in the rigid-lid approximation 
automatically satisfies the surface and bottom boundary conditions 
(no flux through the surface and the bottom: $w_{surface} =w_{bottom} =~0$). 
Indeed, taking the discrete horizontal divergence of the vertical sum of the 
horizontal momentum equations (!!!Eqs. (II.2.1) and (II.2.2)!!!) weighted by the 
vertical scale factors, it becomes:
\begin{flalign*}
\frac{\partial } {\partial t}  \left( \sum\limits_k    \chi    \right)
\equiv 
\frac{\partial } {\partial t}  \left(  w_{surface} -w_{bottom}     \right)&&&\\
\end{flalign*}
\begin{flalign*}
\equiv \frac{1} {e_{1T}\,e_{2T}\,e_{3T}} 
   \biggl\{ \quad
   \delta_i 
      &\left[ 
      e_{2u}\,H_u 
         \left( 
         M_u - M_u - \frac{1} {H_u\,e_{2u}} \delta_j 
            \left[ \partial_t\, \psi \right] 
         \right) 
      \right] &&
   \biggr. \\
   \biggl. 
   + \delta_j 
      &\left[ 
      e_{1v}\,H_v 
         \left( M_v - M_v - \frac{1} {H_v\,e_{1v}} \delta_i 
            \left[ \partial_i\, \psi \right] 
         \right) 
      \right] 
   \biggr\}&& \\
\end{flalign*}
\begin{flalign*}
\equiv \frac{1} {e_{1T} \,e_{2T} \,e_{3T}} \;
   \biggl\{ 
   - \delta_i 
      \Bigl[ 
      \delta_j 
         \left[ \partial_t \psi  \right] 
      \Bigr]
   + \delta_j 
      \Bigl[ 
      \delta_i 
         \left[ \partial_t \psi  \right] 
      \Bigr]\; 
   \biggr\}\;
   \equiv 0
   &&&\\
\end{flalign*}


The surface boundary condition associated with the rigid lid approximation 
($w_{surface} =0)$ is imposed in the computation of the vertical velocity (!!! II.2.5!!!!). 
Therefore, it turns out to be:
\begin{equation*}
\frac{\partial } {\partial t}w_{bottom} \equiv 0
\end{equation*}
As the bottom velocity is initially set to zero, it remains zero all the time. 
Symmetrically, if $w_{bottom} =0$ is used in the computation of the vertical 
velocity (upward integral of the horizontal divergence), the same computation 
leads to $w_{surface} =0$ as soon as the surface vertical velocity is initially 
set to zero.

% -------------------------------------------------------------------------------------------------------------
%       Coriolis and advection terms: vector invariant form
% -------------------------------------------------------------------------------------------------------------
\subsection{Coriolis and advection terms: vector invariant form}
\label{Apdx_C_vor_zad}

% -------------------------------------------------------------------------------------------------------------
%       Vorticity Term
% -------------------------------------------------------------------------------------------------------------
\subsubsection{Vorticity Term}
\label{Apdx_C_vor} 

Potential vorticity is located at $f$-points and defined as: $\zeta / e_{3f}$. 
The standard discrete formulation of the relative vorticity term obviously 
conserves potential vorticity (ENS scheme). It also conserves the potential 
enstrophy for a horizontally non-divergent flow (i.e. $\chi $=0) but not the 
total kinetic energy. Indeed, using the symmetry or skew symmetry properties 
of the operators (Eqs \eqref{DOM_mi_adj} and \eqref{DOM_di_adj}), it can 
be shown that:
\begin{equation} \label{Apdx_C_1.1}
\int_D {\zeta / e_3\,\;{\textbf{k}}\cdot \frac{1} {e_3} \nabla \times \left( {\zeta \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} \equiv 0
\end{equation}
where $dv=e_1\,e_2\,e_3 \; di\,dj\,dk$ is the volume element. Indeed, using 
\eqref{Eq_dynvor_ens}, the discrete form of the right hand side of \eqref{Apdx_C_1.1} 
can be transformed as follow:
\begin{flalign*} 
&\int_D \zeta / e_3\,\; \textbf{k} \cdot \frac{1} {e_3 } \nabla \times 
   \left( 
   \zeta \; \textbf{k} \times  \textbf{U}_h  
   \right)\;
   dv
   &&& \displaybreak[0] \\
%
\equiv& \sum\limits_{i,j,k} 
\frac{\zeta / e_{3f}} {e_{1f}\,e_{2f}\,e_{3f}} 
   \biggl\{ \quad
   \delta_{i+1/2} 
      \left[ 
         - \overline {\left( {\zeta / e_{3f}} \right)}^{\,i}\;
            \overline{\overline {\left( e_{2u}\,e_{3u}\,u \right)}}^{\,i,j+1/ 2} 
       \right]  
   &&  \\ & \qquad \qquad \qquad \;\;
   - \delta_{j+1/2} 
      \left[   \;\;\;
           \overline {\left( \zeta / e_{3f} \right)}^{\,j}\;
           \overline{\overline {\left( e_{1v}\,e_{3v}\,v \right)}}^{\,i+1/2,j} 
      \right] 
   \;\;\biggr\} \;  e_{1f}\,e_{2f}\,e_{3f}       && \displaybreak[0] \\ 
%
\equiv& \sum\limits_{i,j,k} 
   \biggl\{   \delta_i     \left[   \frac{\zeta} {e_{3f}}   \right] \;
           \overline{  \left(   \frac{\zeta} {e_{3f}}   \right)  }^{\,i}\; 
           \overline{  \overline{   \left( e_{1u}\,e_{3u}\,u \right)  }  }^{\,i,j+1/2} 
         + \delta_j   \left[   \frac{\zeta} {e_{3f}}   \right] \;
            \overline{   \left(  \frac{\zeta} {e_{3f}}    \right)  }^{\,j} \; 
      \overline{\overline {\left( e_{2v}\,e_{3v}\,v \right)}}^{\,i+1/2,j}         \biggr\} 
      &&&& \displaybreak[0] \\ 
%
\equiv& \frac{1} {2} \sum\limits_{i,j,k}  
   \biggl\{ \delta_i    \Bigl[    \left(  \frac{\zeta} {e_{3f}} \right)^2   \Bigr]\;
         \overline{\overline {\left( e_{2u}\,e_{3u}\,u \right)}}^{\,i,j+1/2} 
            + \delta_j  \Bigl[    \left( \zeta / e_{3f} \right)^2     \Bigr]\; 
         \overline{\overline {\left( e_{1v}\,e_{3v}\,v \right)}}^{\,i+1/2,j} 
   \biggr\} 
   && \displaybreak[0] \\ 
%
\equiv& - \frac{1} {2} \sum\limits_{i,j,k}   \left(  \frac{\zeta} {e_{3f}} \right)^2\;
   \biggl\{    \delta_{i+1/2} 
         \left[   \overline{\overline {\left( e_{2u}\,e_{3u}\,u \right)}}^{\,i,j+1/2}    \right]  
               + \delta_{j+1/2}
      \left[   \overline{\overline {\left( e_{1v}\,e_{3v}\,v \right)}}^{\,i+1/2,j}     \right]  
   \biggr\}    && \\ 
\end{flalign*}
Since $\overline {\;\cdot \;} $ and $\delta $ operators commute: $\delta_{i+1/2} 
\left[ {\overline a^{\,i}} \right] = \overline {\delta_i \left[ a \right]}^{\,i+1/2}$, 
and introducing the horizontal divergence $\chi $, it becomes:
\begin{align*}
\equiv& \sum\limits_{i,j,k} - \frac{1} {2} \left(  \frac{\zeta} {e_{3f}} \right)^2 \; \overline{\overline{ e_{1T}\,e_{2T}\,e_{3T}\, \chi}}^{\,i+1/2,j+1/2} \;\;\equiv 0 
\qquad \qquad \qquad \qquad \qquad \qquad \qquad &&&&\\
\end{align*}

Note that the derivation is demonstrated here for the relative potential 
vorticity but it applies also to the planetary ($f/e_3$) and the total 
potential vorticity $((\zeta +f) /e_3 )$. Another formulation of the two 
components of the vorticity term is optionally offered (ENE scheme) :
\begin{equation*}
   - \zeta \;{\textbf{k}}\times {\textbf {U}}_h 
\equiv 
   \left( {{
   \begin{array} {*{20}c}
      + \frac{1} {e_{1u}} \; 
      \overline {\left( \zeta / e_{3f}      \right)   
      \overline {\left( e_{1v}\,e_{3v}\,v \right)}^{\,i+1/2}} ^{\,j} 
      \hfill \\
      - \frac{1} {e_{2v}} \; 
      \overline {\left( \zeta / e_{3f}      \right)
      \overline {\left( e_{2u}\,e_{3u}\,u \right)}^{\,j+1/2}} ^{\,i} 
      \hfill \\
   \end{array}} } 
   \right)
\end{equation*}

This formulation does not conserve the enstrophy but it does conserve the 
total kinetic energy. It is also possible to mix the two formulations in order 
to conserve enstrophy on the relative vorticity term and energy on the 
Coriolis term.
\begin{flalign*}
&\int\limits_D - \textbf{U}_h \cdot   \left(  \zeta \;\textbf{k} \times \textbf{U}_h  \right)  \;  dv &&  \\
\equiv& \sum\limits_{i,j,k}   \biggl\{    
      \overline {\left(  \frac{\zeta} {e_{3f}}      \right) 
        \overline {\left( e_{1v}e_{3v}v \right)}^{\,i+1/2}} ^{\,j} \, e_{2u}e_{3u}u 
   - \overline {\left(  \frac{\zeta} {e_{3f}}       \right) 
       \overline {\left( e_{2u}e_{3u}u \right)}^{\,j+1/2}} ^{\,i} \, e_{1v}e_{3v}v \;
                                   \biggr\}     
\\
\equiv& \sum\limits_{i,j,k}  \frac{\zeta} {e_{3f}}
   \biggl\{  \overline {\left( e_{1v}e_{3v} v \right)}^{\,i+1/2}\;
             \overline {\left( e_{2u}e_{3u} u \right)}^{\,j+1/2}        
        - \overline {\left( e_{2u}e_{3u} u \right)}^{\,j+1/2}\;
               \overline {\left( e_{1v}e_{3v} v \right)}^{\,i+1/2}
   \biggr\}
   \equiv 0
\end{flalign*}


% -------------------------------------------------------------------------------------------------------------
%       Gradient of Kinetic Energy / Vertical Advection
% -------------------------------------------------------------------------------------------------------------
\subsubsection{Gradient of Kinetic Energy / Vertical Advection}
\label{Apdx_C_zad} 

The change of Kinetic Energy (KE) due to the vertical advection is exactly 
balanced by the change of KE due to the horizontal gradient of KE~:
\begin{equation*}
      \int_D \textbf{U}_h \cdot \nabla_h \left( \frac{1}{2}\;{\textbf{U}_h}^2 \right)\;dv
 = - \int_D \textbf{U}_h \cdot w \frac{\partial \textbf{U}_h} {\partial k}\;dv
\end{equation*}
Indeed, using successively \eqref{DOM_di_adj} ($i.e.$ the skew symmetry 
property of the $\delta$ operator) and the incompressibility, then 
\eqref{DOM_di_adj} again, then the commutativity of operators 
$\overline {\,\cdot \,}$ and $\delta$, and finally \eqref{DOM_mi_adj} 
($i.e.$ the symmetry property of the $\overline {\,\cdot \,}$ operator) 
applied in the horizontal and vertical directions, it becomes:
\begin{flalign*}
&\int_D \textbf{U}_h \cdot \nabla_h \left( \frac{1}{2}\;{\textbf{U}_h}^2 \right)\;dv   &&&\\
\equiv& \frac{1}{2} \sum\limits_{i,j,k} 
   \biggl\{ 
   \frac{1} {e_{1u}}  \delta_{i+1/2} 
   \left[   \overline {u^2}^{\,i} + \overline {v^2}^{\,j}    \right]  u\,e_{1u}e_{2u}e_{3u}  
     + \frac{1} {e_{2v}}  \delta_{j+1/2} 
   \left[   \overline {u^2}^{\,i} + \overline {v^2}^{\,j}   \right]  v\,e_{1v}e_{2v}e_{3v} 
   \biggr\} 
   &&& \displaybreak[0] \\ 
%
\equiv&  \frac{1}{2} \sum\limits_{i,j,k} 
   \left(   \overline {u^2}^{\,i} + \overline {v^2}^{\,j}   \right)\;
   \delta_k \left[  e_{1T}\,e_{2T} \,w   \right]
%
\;\; \equiv -\frac{1}{2} \sum\limits_{i,j,k}  \delta_{k+1/2} 
   \left[ 
      \overline{ u^2}^{\,i} 
   + \overline{ v^2}^{\,j} 
   \right] \;
   e_{1v}\,e_{2v}\,w 
   &&& \displaybreak[0]\\
%
\equiv &\frac{1} {2} \sum\limits_{i,j,k} 
   \left(    \overline {\delta_{k+1/2} \left[ u^2 \right]}^{\,i} 
      + \overline {\delta_{k+1/2} \left[ v^2 \right]}^{\,j}    \right) \; e_{1T}\,e_{2T} \,w 
   && \displaybreak[0] \\
%  
\equiv &\frac{1} {2} \sum\limits_{i,j,k} 
   \biggl\{  \overline {e_{1T}\,e_{2T} \,w}^{\,i+1/2}\;2    
         \overline {u}^{\,k+1/2}\; \delta_{k+1/2}         \left[ u \right]     %&&&  \\
    + \overline {e_{1T}\,e_{2T} \,w}^{\,j+1/2}\;2  \overline {v}^{\,k+1/2}\; \delta_{k+1/2}  \left[ v \right]  \;
   \biggr\} 
   &&\displaybreak[0] \\ 
%
\equiv& -\sum\limits_{i,j,k}  
   \biggl\{
   \quad \frac{1} {b_u } \;
   \overline {\Bigl\{ \overline {e_{1T}\,e_{2T} \,w}^{\,i+1/2}\,\delta_{k+1/2}
      \left[ u \right] 
             \Bigr\} }^{\,k} \;u\;e_{1u}\,e_{2u}\,e_{3u}  
   && \\
   &\qquad \quad\; + \frac{1} {b_v } \;
   \overline {\Bigl\{ \overline {e_{1T}\,e_{2T} \,w}^{\,j+1/2} \delta_{k+1/2}
       \left[ v \right] 
         \Bigr\} }^{\,k} \;v\;e_{1v}\,e_{2v}\,e_{3v}  \;
   \biggr\} 
   && \\ 
\equiv& -\int\limits_D \textbf{U}_h \cdot w \frac{\partial \textbf{U}_h} {\partial k}\;dv &&&\\
\end{flalign*}

The main point here is that the satisfaction of this property links the choice of 
the discrete formulation of the vertical advection and of the horizontal gradient 
of KE. Choosing one imposes the other. For example KE can also be discretized 
as $1/2\,({\overline u^{\,i}}^2 + {\overline v^{\,j}}^2)$. This leads to the following 
expression for the vertical advection:
\begin{equation*}
\frac{1} {e_3 }\; w\; \frac{\partial \textbf{U}_h } {\partial k}
\equiv \left( {{\begin{array} {*{20}c}
\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} \;
\overline{\overline {e_{1T}\,e_{2T} \,w\;\delta_{k+1/2} 
\left[ \overline u^{\,i+1/2} \right]}}^{\,i+1/2,k}  \hfill \\
\frac{1} {e_{1v}\,e_{2v}\,e_{3v}} \;
\overline{\overline {e_{1T}\,e_{2T} \,w\;\delta_{k+1/2}
\left[ \overline v^{\,j+1/2} \right]}}^{\,j+1/2,k} \hfill \\
\end{array}} } \right)
\end{equation*}
a formulation that requires an additional horizontal mean in contrast with 
the one used in NEMO. Nine velocity points have to be used instead of 3. 
This is the reason why it has not been chosen.

% -------------------------------------------------------------------------------------------------------------
%       Coriolis and advection terms: flux form
% -------------------------------------------------------------------------------------------------------------
\subsection{Coriolis and advection terms: flux form}
\label{Apdx_C.1.3}

% -------------------------------------------------------------------------------------------------------------
%       Coriolis plus ``metric'' Term
% -------------------------------------------------------------------------------------------------------------
\subsubsection{Coriolis plus ``metric'' Term}
\label{Apdx_C.1.3.1} 

In flux from the vorticity term reduces to a Coriolis term in which the Coriolis 
parameter has been modified to account for the ``metric'' term. This altered 
Coriolis parameter is discretised at an f-point. It is given by: 
\begin{equation*}
f+\frac{1} {e_1 e_2 } 
\left( v \frac{\partial e_2 } {\partial i} - u \frac{\partial e_1 } {\partial j}\right)\;
\equiv \;
f+\frac{1} {e_{1f}\,e_{2f}} 
\left( \overline v^{\,i+1/2} \delta_{i+1/2} \left[ e_{2u} \right] 
-\overline u^{\,j+1/2} \delta_{j+1/2} \left[ e_{1u}  \right] \right)
\end{equation*}

The ENE scheme is then applied to obtain the vorticity term in flux form. 
It therefore conserves the total KE. The derivation is the same as for the 
vorticity term in the vector invariant form (\S\ref{Apdx_C_vor}).

% -------------------------------------------------------------------------------------------------------------
%       Flux form advection
% -------------------------------------------------------------------------------------------------------------
\subsubsection{Flux form advection}
\label{Apdx_C.1.3.2} 

The flux form operator of the momentum advection is evaluated using a 
centered second order finite difference scheme. Because of the flux form, 
the discrete operator does not contribute to the global budget of linear 
momentum. Because of the centered second order scheme, it conserves 
the horizontal kinetic energy, that is :

\begin{equation} \label{Apdx_C_I.3.10}
\int_D \textbf{U}_h \cdot 
\left( {{\begin{array} {*{20}c}
\nabla \cdot \left( \textbf{U}\,u \right) \hfill \\
\nabla \cdot \left( \textbf{U}\,v \right) \hfill \\
\end{array}} } \right)\;dv =\;0
\end{equation}

Let us demonstrate this property for the first term of the scalar product 
($i.e.$ considering just the the terms associated with the i-component of 
the advection):
\begin{flalign*}
&\int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv    &&&\\
%
\equiv& \sum\limits_{i,j,k} 
\biggl\{    \frac{1} {e_{1u}\, e_{2u}\,e_{3u}}    \biggl(   
      \delta_{i+1/2}  \left[   \overline {e_{2u}\,e_{3u}\,u}^{\,i}      \;\overline u^{\,i}          \right]   
   + \delta_j           \left[   \overline {e_{1u}\,e_{3u}\,v}^{\,i+1/2}\;\overline u^{\,j+1/2}   \right] 
      &&& \\ & \qquad \qquad \qquad \qquad \qquad \qquad
   + \delta_k          \left[   \overline {e_{1w}\,e_{2w}\,w}^{\,i+1/2}\;\overline u^{\,k+1/2} \right]
         \biggr)   \;   \biggr\} \, e_{1u}\,e_{2u}\,e_{3u} \;u 
      &&& \displaybreak[0] \\ 
%
\equiv& \sum\limits_{i,j,k} 
   \biggl\{ 
      \delta_{i+1/2} \left[   \overline {e_{2u}\,e_{3u}\,u}^{\,i}\;  \overline u^{\,i}  \right]
   + \delta_j          \left[   \overline {e_{1u}\,e_{3u}\,v}^{\,i+1/2}\;\overline u^{\,j+1/2}   \right]
      &&& \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad
   + \delta_k         \left[   \overline {e_{1w}\,e_{2w}\,w}^{\,i+12}\;\overline u^{\,k+1/2}  \right]
      \; \biggr\} \; u    &&& \displaybreak[0] \\
%
\equiv& - \sum\limits_{i,j,k}
   \biggl\{ 
      \overline {e_{2u}\,e_{3u}\,u}^{\,i}\;        \overline u^{\,i}       \delta_i 
      \left[ u \right] 
        + \overline {e_{1u}\,e_{3u}\,v}^{\,i+1/2}\;   \overline u^{\,j+1/2}   \delta_{j+1/2} 
      \left[ u \right] 
      &&& \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad
       + \overline {e_{1w}\,e_{2w}\,w}^{\,i+1/2}\; \overline u^{\,k+1/2}   \delta_{k+1/2}    \left[ u \right]     \biggr\}     && \displaybreak[0] \\
%
\equiv& - \sum\limits_{i,j,k}
   \biggl\{ 
        \overline {e_{2u}\,e_{3u}\,u}^{\,i}        \delta_i       \left[ u^2 \right] 
    + \overline {e_{1u}\,e_{3u}\,v}^{\,i+1/2}      \delta_{j+/2}  \left[ u^2 \right]
    + \overline {e_{1w}\,e_{2w}\,w}^{\,i+1/2}   \delta_{k+1/2}    \left[ u^2 \right] 
   \biggr\} 
   && \displaybreak[0] \\
%
\equiv& \sum\limits_{i,j,k}
   \bigg\{ 
      e_{2u}\,e_{3u}\,u\;     \delta_{i+1/2}       \left[ \overline {u^2}^{\,i} \right]
        + e_{1u}\,e_{3u}\,v\; \delta_{j+1/2}    \; \left[ \overline {u^2}^{\,i} \right]
    + e_{1w}\,e_{2w}\,w\;  \delta_{k+1/2}       \left[ \overline {u^2}^{\,i} \right] 
   \biggr\} 
   && \displaybreak[0] \\
%
\equiv& \sum\limits_{i,j,k}
\overline {u^2}^{\,i} 
   \biggl\{ 
      \delta_{i+1/2}    \left[ e_{2u}\,e_{3u}\,u  \right]
   + \delta_{j+1/2}  \left[ e_{1u}\,e_{3u}\,v  \right]
   + \delta_{k+1/2}  \left[ e_{1w}\,e_{2w}\,w \right] 
   \biggr\}  \;\;  \equiv 0
   &&& \\
\end{flalign*}

When the UBS scheme is used to evaluate the flux form momentum advection, 
the discrete operator does not contribute to the global budget of linear momentum 
(flux form). The horizontal kinetic energy is not conserved, but forced to decay 
($i.e.$ the scheme is diffusive). 

% -------------------------------------------------------------------------------------------------------------
%       Hydrostatic Pressure Gradient Term
% -------------------------------------------------------------------------------------------------------------
\subsection{Hydrostatic Pressure Gradient Term}
\label{Apdx_C.1.4}


A pressure gradient has no contribution to the evolution of the vorticity as the 
curl of a gradient is zero. In the $z$-coordinate, this property is satisfied locally 
on a C-grid with 2nd order finite differences (property \eqref{Eq_DOM_curl_grad}). 
When the equation of state is linear ($i.e.$ when an advection-diffusion equation 
for density can be derived from those of temperature and salinity) the change of 
KE due to the work of pressure forces is balanced by the change of potential 
energy due to buoyancy forces: 
\begin{equation*}
\int_D - \frac{1} {\rho_o} \left. \nabla p^h \right|_z \cdot \textbf{U}_h \;dv 
= \int_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
\end{equation*}

This property can be satisfied in a discrete sense for both $z$- and $s$-coordinates. 
Indeed, defining the depth of a $T$-point, $z_T$, as the sum of the vertical scale 
factors at $w$-points starting from the surface, the work of pressure forces can be 
written as:
\begin{flalign*}
&\int_D - \frac{1} {\rho_o} \left. \nabla p^h \right|_z \cdot \textbf{U}_h \;dv   &&& \\
\equiv& \sum\limits_{i,j,k} \biggl\{ \;  - \frac{1} {\rho_o e_{1u}}   \Bigl( 
\delta_{i+1/2} \left[ p^h \right] - g\;\overline \rho^{\,i+1/2}\;\delta_{i+1/2} \left[ z_T \right] 
               \Bigr)  \; u\;e_{1u}\,e_{2u}\,e_{3u}  
   &&  \\ & \qquad \qquad
                                 - \frac{1} {\rho_o e_{2v}}    \Bigl( 
\delta_{j+1/2} \left[ p^h \right] - g\;\overline \rho^{\,j+1/2}\delta_{j+1/2}  \left[ z_T \right] 
               \Bigr)  \; v\;e_{1v}\,e_{2v}\,e_{3v} \;
   \biggr\}   && \\ 
\end{flalign*}

Using  \eqref{DOM_di_adj}, $i.e.$ the skew symmetry property of the $\delta$ 
operator, \eqref{Eq_wzv}, the continuity equation), and \eqref{Eq_dynhpg_sco}, 
the hydrostatic equation in the $s$-coordinate, it becomes: 
\begin{flalign*} 
\equiv& \frac{1} {\rho_o} \sum\limits_{i,j,k}    \biggl\{ 
      e_{2u}\,e_{3u} \;u\;g\; \overline \rho^{\,i+1/2}\;\delta_{i+1/2}[ z_T]    
   +  e_{1v}\,e_{3v} \;v\;g\; \overline \rho^{\,j+1/2}\;\delta_{j+1/2}[ z_T]     
&& \\  & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \;\;\,
   +\Bigl(  \delta_i[ e_{2u}\,e_{3u}\,u] + \delta_j [ e_{1v}\,e_{3v}\,v]  \Bigr)\;p^h \biggr\}  &&\\
%
\equiv& \frac{1} {\rho_o } \sum\limits_{i,j,k}
   \biggl\{ 
       e_{2u}\,e_{3u} \;u\;g\;   \overline \rho^{\,i+1/2} \delta_{i+1/2} \left[ z_T \right]
   +  e_{1v}\,e_{3v} \;v\;g\; \overline \rho^{\,j+1/2} \delta_{j+1/2} \left[ z_T \right] 
   &&&\\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \;\;\,
    - \delta_k \left[ e_{1w} e_{2w}\,w \right]\;p^h   \biggr\}   &&&\\ 
%
\equiv& \frac{1} {\rho_o } \sum\limits_{i,j,k}
   \biggl\{ 
      e_{2u}\,e_{3u} \;u\;g\; \overline \rho^{\,i+1/2}\;\delta_{i+1/2} \left[ z_T \right]
   +  e_{1v}\,e_{3v} \;v\;g\; \overline \rho^{\,j+1/2} \;\delta_{j+1/2} \left[ z_T \right] 
   &&& \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \;\;\,
   +   e_{1w} e_{2w} \;w\;\delta_{k+1/2} \left[ p_h \right] 
   \biggr\}  &&&\\ 
%
\equiv& \frac{g} {\rho_o}  \sum\limits_{i,j,k}
   \biggl\{ 
      e_{2u}\,e_{3u} \;u\; \overline \rho^{\,i+1/2}\;\delta_{i+1/2} \left[ z_T \right]
   +  e_{1v}\,e_{3v} \;v\; \overline \rho^{\,j+1/2}\;\delta_{j+1/2} \left[ z_T \right]    
   &&& \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \;\;\,
   -  e_{1w} e_{2w} \;w\;e_{3w} \overline \rho^{\,k+1/2} 
   \biggr\}   &&&\\ 
\end{flalign*}
noting that by definition of $z_T$, $\delta_{k+1/2} \left[ z_T \right] \equiv - e_{3w} $, thus:
\begin{multline*}
\equiv \frac{g} {\rho_o}  \sum\limits_{i,j,k}
   \biggl\{ 
      e_{2u}\,e_{3u} \;u\; \overline \rho^{\,i+1/2}\;\delta_{i+1/2} \left[ z_T \right]
   +  e_{1v}\,e_{3v} \;v\; \overline \rho^{\,j+1/2} \delta_{j+1/2} \left[ z_T \right]
   \biggr. \\ 
\shoveright{
   \biggl. 
   +  e_{1w} e_{2w} \;w\;  \overline \rho^{\,k+1/2}\;\delta_{k+1/2} \left[ z_T \right] 
   \biggr\} } \\ 
\end{multline*}
Using \eqref{DOM_di_adj}, it becomes:
\begin{flalign*}
\equiv& - \frac{g} {\rho_o} \sum\limits_{i,j,k} z_T 
   \biggl\{ 
      \delta_i    \left[ e_{2u}\,e_{3u}\,u\; \overline \rho^{\,i+1/2}   \right]
   +  \delta_j    \left[ e_{1v}\,e_{3v}\,v\; \overline \rho^{\,j+1/2}   \right]
   +  \delta_k    \left[ e_{1w} e_{2w}\,w\;  \overline \rho^{\,k+1/2}   \right] 
   \biggr\} 
   &&& \\
%
\equiv& -\int_D \nabla \cdot \left( \rho \, \textbf{U} \right)\;g\;z\;\;dv    &&& \\
\end{flalign*}

Note that this property strongly constrains the discrete expression of both 
the depth of $T-$points and of the term added to the pressure gradient in the
$s$-coordinate. Nevertheless, it is almost never satisfied since a linear equation 
of state is rarely used.

% -------------------------------------------------------------------------------------------------------------
%       Surface Pressure Gradient Term
% -------------------------------------------------------------------------------------------------------------
\subsection{Surface Pressure Gradient Term}
\label{Apdx_C.1.5}


The surface pressure gradient has no contribution to the evolution of the vorticity. 
This property is trivially satisfied locally since the equation verified by $\psi$ has 
been derived from the discrete formulation of the momentum equation and of the curl. 
But it has to be noted that since the elliptic equation satisfied by $\psi$ is solved 
numerically by an iterative solver (preconditioned conjugate gradient or successive 
over relaxation), the property is only satisfied at the precision requested for the 
solver used.

With the rigid-lid approximation, the change of KE due to the work of surface 
pressure forces is exactly zero. This is satisfied in discrete form, at the precision 
requested for the elliptic solver used to solve this equation. This can be 
demonstrated as follows:
\begin{flalign*}
\int\limits_D  - \frac{1} {\rho_o} \nabla_h \left( p_s \right) \cdot \textbf{U}_h \;dv%   &&& \\
%
&\equiv \sum\limits_{i,j,k}   \biggl\{    \;
    \left(  - M_u - \frac{1} {H_u \,e_{2u}}  \delta_j    \left[ \partial_t \psi  \right]   \right)\;
    u\;e_{1u}\,e_{2u}\,e_{3u}  
&&&\\&  \qquad \;\;\,
      + \left(  - M_v + \frac{1} {H_v \,e_{1v}}  \delta_i \left[ \partial_t \psi \right]     \right)\;
     v\;e_{1v}\,e_{2v}\,e_{3v}   \; \biggr\}      
&&&\\
\\
%
&\equiv \sum\limits_{i,j}  \Biggl\{   \;
   \biggl( - M_u - \frac{1} {H_u \,e_{2u}}  \delta_j \left[ \partial_t \psi  \right]   \biggr)
   \biggl( \sum\limits_k  u\;e_{3u}   \biggr)\;  e_{1u}\,e_{2u} 
&&&\\&  \qquad \;\;\,
   + \biggl( - M_v + \frac{1} {H_v \,e_{1v}}  \delta_i \left[ \partial_t \psi  \right]   \biggr)
      \biggl(   \sum\limits_k   v\;e_{3v}   \biggr)\;   e_{1v}\,e_{2v} \;   \Biggr\} 
   && \\ 
%
\intertext{using the relation between \textit{$\psi $} and the vertical sum of the velocity, it becomes:}
%
&\equiv \sum\limits_{i,j} 
   \biggl\{  \;   
      \left( \;\;\,
      M_u + \frac{1} {H_u \,e_{2u}}  \delta_j \left[ \partial_t \psi \right] 
      \right)\;
      e_{1u} \,\delta_j 
         \left[ \partial_t \psi  \right] 
   && \\ &  \qquad \;\;\,
      + \left( 
      - M_v + \frac{1} {H_v \,e_{1v}}  \delta_i \left[ \partial_t \psi \right] 
      \right)\;
      e_{2v} \,\delta_i \left[ \partial_t \psi \right]   \;
   \biggr\} 
   && \\ 
%
\intertext{applying the adjoint of the $\delta$ operator, it is now:}
%
&\equiv \sum\limits_{i,j}  - \partial_t \psi \;
   \biggl\{    \;
     \delta_{j+1/2} \left[ e_{1u} M_u \right] 
     - \delta_{i+1/2} \left[ e_{1v} M_v \right] 
   && \\ &  \qquad \;\;\,
   + \delta_{i+1/2} 
      \left[ \frac{e_{2v}} {H_v \,e_{2v}}  \delta_i \left[ \partial_t \psi \right] 
      \right] 
   + \delta_{j+1/2} 
       \left[ \frac{e_{1u}} {H_u \,e_{2u}}  \delta_j \left[ \partial_t \psi \right] 
       \right]  
   \biggr\}   &&&\\
   &\equiv 0                   && \\ 
\end{flalign*}

The last equality is obtained using \eqref{Eq_dynspg_rl}, the discrete barotropic 
streamfunction time evolution equation. By the way, this shows that 
\eqref{Eq_dynspg_rl} is the only way to compute the streamfunction, 
otherwise the surface pressure forces will do work. Nevertheless, since 
the elliptic equation satisfied by $\psi $ is solved numerically by an iterative 
solver, the property is only satisfied at the precision requested for the solver.

% ================================================================
% Conservation Properties on Tracers
% ================================================================
\section{Conservation Properties on Tracers}
\label{Apdx_C.2}


All the numerical schemes used in NEMO are written such that the tracer content 
is conserved by the internal dynamics and physics (equations in flux form). 
For advection, only the CEN2 scheme ($i.e.$ $2^{nd}$ order finite different scheme) 
conserves the global variance of tracer. Nevertheless the other schemes ensure 
that the global variance decreases ($i.e.$ they are at least slightly diffusive). 
For diffusion, all the schemes ensure the decrease of the total tracer variance, 
except the iso-neutral operator. There is generally no strict conservation of mass, 
as the equation of state is non linear with respect to $T$ and $S$. In practice, 
the mass is conserved to a very high accuracy. 
% -------------------------------------------------------------------------------------------------------------
%       Advection Term
% -------------------------------------------------------------------------------------------------------------
\subsection{Advection Term}
\label{Apdx_C.2.1}

Whatever the advection scheme considered it conserves of the tracer content as all 
the scheme are written in flux form. Let $\tau$ be the tracer interpolated at velocity point 
(whatever the interpolation is). The conservation of the tracer content is obtained as follows: 
\begin{flalign*}
&\int_D \nabla \cdot \left( T \textbf{U} \right)\;dv    &&&\\
&\equiv  \sum\limits_{i,j,k}    \biggl\{
    \frac{1} {e_{1T}\,e_{2T}\,e_{3T}} 
    \left(  \delta_i    \left[   e_{2u}\,e_{3u}\; u \;\tau_u   \right]
           + \delta_j    \left[   e_{1v}\,e_{3v}\; v  \;\tau_v   \right] \right) 
&&&\\& \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad
   + \frac{1} {e_{3T}} \delta_k \left[ w\;\tau \right]    \biggl\}  e_{1T}\,e_{2T}\,e_{3T} &&&\\
%
&\equiv  \sum\limits_{i,j,k}     \left\{
      \delta_i  \left[ e_{2u}\,e_{3u}  \,\overline T^{\,i+1/2}\,u \right]
         + \delta_j  \left[ e_{1v}\,e_{3v}  \,\overline T^{\,j+1/2}\,v \right] 
   + \delta_k \left[ e_{1T}\,e_{2T} \,\overline T^{\,k+1/2}\,w \right] \right\} 
    && \\
&\equiv 0 &&&
\end{flalign*}

The conservation of the variance of tracer can be achieved only with the CEN2 scheme. 
It can be demonstarted as follows:
\begin{flalign*}
&\int\limits_D T\;\nabla \cdot \left( T\; \textbf{U} \right)\;dv &&&\\
\equiv& \sum\limits_{i,j,k} T\;
   \left\{
      \delta_i  \left[ e_{2u}\,e_{3u} \overline T^{\,i+1/2}\,u \right]
   + \delta_j  \left[ e_{1v}\,e_{3v} \overline T^{\,j+1/2}\,v \right]
   + \delta_k \left[ e_{1T}\,e_{2T} \overline T^{\,k+1/2}w \right]
   \right\} 
   && \\
\equiv& \sum\limits_{i,j,k} 
   \left\{
   -           e_{2u}\,e_{3u} \overline T^{\,i+1/2}\,u\,\delta_{i+1/2}  \left[ T \right] \right.
   -           e_{1v}\,e_{3v}  \overline T^{\,j+1/2}\,v\;\delta_{j+1/2}  \left[ T \right]
&&&\\& \qquad \qquad \qquad \qquad \qquad \qquad \quad \;
   - \left. e_{1T}\,e_{2T} \overline T^{\,k+1/2}w\;\delta_{k+1/2} \left[ T \right]
   \right\} 
   &&\\
\equiv&  -\frac{1} {2}  \sum\limits_{i,j,k}
   \Bigl\{
      e_{2u}\,e_{3u} \;  u\;\delta_{i+1/2} \left[ T^2 \right]
   + e_{1v}\,e_{3v} \;  v\;\delta_{j+1/2}  \left[ T^2 \right]
   + e_{1T}\,e_{2T} \;w\;\delta_{k+1/2} \left[ T^2 \right]
   \Bigr\} 
   && \\ 
\equiv& \frac{1} {2}  \sum\limits_{i,j,k} T^2
   \Bigl\{
      \delta_i  \left[ e_{2u}\,e_{3u}\,u \right]
   + \delta_j  \left[ e_{1v}\,e_{3v}\,v \right]
   + \delta_k \left[ e_{1T}\,e_{2T}\,w \right]
   \Bigr\} 
\quad \equiv 0 &&&
\end{flalign*}


% ================================================================
% Conservation Properties on Lateral Momentum Physics
% ================================================================
\section{Conservation Properties on Lateral Momentum Physics}
\label{Apdx_dynldf_properties}


The discrete formulation of the horizontal diffusion of momentum ensures the 
conservation of potential vorticity and the horizontal divergence, and the 
dissipation of the square of these quantities (i.e. enstrophy and the 
variance of the horizontal divergence) as well as the dissipation of the 
horizontal kinetic energy. In particular, when the eddy coefficients are 
horizontally uniform, it ensures a complete separation of vorticity and 
horizontal divergence fields, so that diffusion (dissipation) of vorticity 
(enstrophy) does not generate horizontal divergence (variance of the 
horizontal divergence) and \textit{vice versa}. 

These properties of the horizontal diffusion operator are a direct consequence 
of properties \eqref{Eq_DOM_curl_grad} and \eqref{Eq_DOM_div_curl}. 
When the vertical curl of the horizontal diffusion of momentum (discrete sense) 
is taken, the term associated with the horizontal gradient of the divergence is 
locally zero. 

% -------------------------------------------------------------------------------------------------------------
%       Conservation of Potential Vorticity
% -------------------------------------------------------------------------------------------------------------
\subsection{Conservation of Potential Vorticity}
\label{Apdx_C.3.1}

The lateral momentum diffusion term conserves the potential vorticity :
\begin{flalign*}
&\int \limits_D \frac{1} {e_3 } \textbf{k} \cdot \nabla \times 
   \Bigl[ \nabla_h 
      \left( A^{\,lm}\;\chi  \right)
   - \nabla_h \times
      \left( A^{\,lm}\;\zeta \; \textbf{k} \right)
   \Bigr]\;dv  = 0
\end{flalign*}
%%%%%%%%%%  recheck here....  (gm)
\begin{flalign*}
= \int \limits_D  -\frac{1} {e_3 } \textbf{k} \cdot \nabla \times 
   \Bigl[ \nabla_h \times 
      \left( A^{\,lm}\;\zeta \; \textbf{k} \right) 
   \Bigr]\;dv &&& \\ 
\end{flalign*}
\begin{flalign*}
\equiv& \sum\limits_{i,j}
   \left\{
   \delta_{i+1/2} 
   \left[ 
   \frac {e_{2v}} {e_{1v}\,e_{3v}}  \delta_i
      \left[ A_f^{\,lm} e_{3f} \zeta  \right]
    \right]
   + \delta_{j+1/2} 
   \left[ 
   \frac {e_{1u}} {e_{2u}\,e_{3u}} \delta_j 
      \left[ A_f^{\,lm} e_{3f} \zeta  \right]
   \right]
   \right\} 
   && \\ 
%
\intertext{Using \eqref{DOM_di_adj}, it follows:}
%
\equiv& \sum\limits_{i,j,k} 
   -\,\left\{
      \frac{e_{2v}} {e_{1v}\,e_{3v}}  \delta_i
      \left[ A_f^{\,lm} e_{3f} \zeta  \right]\;\delta_i \left[ 1\right]
   + \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j 
      \left[ A_f^{\,lm} e_{3f} \zeta  \right]\;\delta_j \left[ 1\right]
   \right\} \quad \equiv 0 
   && \\ 
\end{flalign*}

% -------------------------------------------------------------------------------------------------------------
%       Dissipation of Horizontal Kinetic Energy
% -------------------------------------------------------------------------------------------------------------
\subsection{Dissipation of Horizontal Kinetic Energy}
\label{Apdx_C.3.2}


The lateral momentum diffusion term dissipates the horizontal kinetic energy:
%\begin{flalign*}
\begin{equation*}
\begin{split}
\int_D \textbf{U}_h \cdot 
   \left[ \nabla_h      \right.   &     \left.       \left( A^{\,lm}\;\chi \right)      
   - \nabla_h \times  \left( A^{\,lm}\;\zeta \;\textbf{k} \right)     \right] \; dv    \\
\\  %%%
\equiv& \sum\limits_{i,j,k} 
   \left\{
     \frac{1} {e_{1u}}               \delta_{i+1/2} \left[  A_T^{\,lm}          \chi     \right]
   - \frac{1} {e_{2u}\,e_{3u}}  \delta_j           \left[ A_f^{\,lm} e_{3f} \zeta   \right]
   \right\} \; e_{1u}\,e_{2u}\,e_{3u} \;u     \\
&\;\; +  \left\{
      \frac{1} {e_{2u}}             \delta_{j+1/2} \left[ A_T^{\,lm}          \chi    \right] 
   + \frac{1} {e_{1v}\,e_{3v}} \delta_i            \left[ A_f^{\,lm} e_{3f} \zeta  \right] 
   \right\} \; e_{1v}\,e_{2u}\,e_{3v} \;v     \qquad \\ 
\\  %%%
\equiv& \sum\limits_{i,j,k} 
   \Bigl\{
     e_{2u}\,e_{3u} \;u\;  \delta_{i+1/2} \left[ A_T^{\,lm}           \chi    \right]
   - e_{1u}             \;u\;  \delta_j           \left[ A_f^{\,lm} e_{3f} \zeta  \right]
    \Bigl\} 
    \\ 
&\;\; + \Bigl\{
      e_{1v}\,e_{3v} \;v\;  \delta_{j+1/2}  \left[ A_T^{\,lm}           \chi    \right]
   + e_{2v}             \;v\;  \delta_i            \left[ A_f^{\,lm} e_{3f} \zeta  \right]
   \Bigl\}      \\ 
\\  %%%
\equiv& \sum\limits_{i,j,k} 
   - \Bigl(
     \delta_i   \left[  e_{2u}\,e_{3u} \;u  \right]
   + \delta_j  \left[  e_{1v}\,e_{3v}  \;v  \right] 
        \Bigr) \;  A_T^{\,lm} \chi   \\ 
&\;\; - \Bigl(
     \delta_{i+1/2}  \left[  e_{2v}  \;v  \right]
   - \delta_{j+1/2}  \left[  e_{1u} \;u  \right] 
        \Bigr)\;  A_f^{\,lm} e_{3f} \zeta      \\ 
\\  %%%
\equiv& \sum\limits_{i,j,k} 
   - A_T^{\,lm}  \,\chi^2   \;e_{1T}\,e_{2T}\,e_{3T}
   - A_f ^{\,lm}  \,\zeta^2 \;e_{1f }\,e_{2f }\,e_{3f}  
\quad \leq 0       \\
\end{split}
\end{equation*}

% -------------------------------------------------------------------------------------------------------------
%       Dissipation of Enstrophy
% -------------------------------------------------------------------------------------------------------------
\subsection{Dissipation of Enstrophy}
\label{Apdx_C.3.3}


The lateral momentum diffusion term dissipates the enstrophy when the eddy 
coefficients are horizontally uniform:
\begin{flalign*}
&\int\limits_D  \zeta \; \textbf{k} \cdot \nabla \times 
   \left[ 
     \nabla_h 
      \left( A^{\,lm}\;\chi  \right)
   -\nabla_h \times 
      \left( A^{\,lm}\;\zeta \; \textbf{k} \right)
   \right]\;dv &&&\\
&= A^{\,lm} \int \limits_D \zeta \textbf{k} \cdot \nabla \times 
   \left[ 
   \nabla_h \times 
      \left( \zeta \; \textbf{k} \right)
   \right]\;dv &&&\displaybreak[0]\\
&\equiv A^{\,lm} \sum\limits_{i,j,k}  \zeta \;e_{3f} 
   \left\{
     \delta_{i+1/2} 
   \left[ 
   \frac{e_{2v}} {e_{1v}\,e_{3v}} \delta_i 
      \left[ e_{3f} \zeta  \right]
   \right]
   + \delta_{j+1/2} 
   \left[
   \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j 
      \left[ e_{3f} \zeta  \right]
   \right]
   \right\} 
   &&&\\ 
%
\intertext{Using \eqref{DOM_di_adj}, it follows:}
%
&\equiv  - A^{\,lm} \sum\limits_{i,j,k} 
   \left\{
     \left( 
     \frac{1} {e_{1v}\,e_{3v}}  \delta_i 
      \left[ e_{3f} \zeta  \right] 
     \right)^2   e_{1v}\,e_{2v}\,e_{3v}
   + \left(
     \frac{1} {e_{2u}\,e_{3u}}  \delta_j 
      \left[ e_{3f} \zeta  \right]
     \right)^2   e_{1u}\,e_{2u}\,e_{3u}
     \right\}      &&&\\
& \leq \;0       &&&\\ 
\end{flalign*}

% -------------------------------------------------------------------------------------------------------------
%       Conservation of Horizontal Divergence
% -------------------------------------------------------------------------------------------------------------
\subsection{Conservation of Horizontal Divergence}
\label{Apdx_C.3.4}

When the horizontal divergence of the horizontal diffusion of momentum 
(discrete sense) is taken, the term associated with the vertical curl of the 
vorticity is zero locally, due to (!!! II.1.8  !!!!!). The resulting term conserves the 
$\chi$ and dissipates $\chi^2$ when the eddy coefficients are 
horizontally uniform.
\begin{flalign*}
& \int\limits_D  \nabla_h \cdot 
   \Bigl[ 
     \nabla_h 
      \left( A^{\,lm}\;\chi \right)
   - \nabla_h \times 
      \left( A^{\,lm}\;\zeta \;\textbf{k} \right)
   \Bigr]
   dv
= \int\limits_D  \nabla_h \cdot \nabla_h 
   \left( A^{\,lm}\;\chi  \right)
   dv 
&&&\\
&\equiv \sum\limits_{i,j,k} 
   \left\{ 
     \delta_i
      \left[ 
      A_u^{\,lm} \frac{e_{2u}\,e_{3u}} {e_{1u}}  \delta_{i+1/2} 
         \left[ \chi \right] 
      \right]
   + \delta_j 
      \left[ 
      A_v^{\,lm} \frac{e_{1v}\,e_{3v}} {e_{2v}}  \delta_{j+1/2} 
         \left[ \chi \right] 
      \right] 
   \right\}
   &&&\\ 
%
\intertext{Using \eqref{DOM_di_adj}, it follows:}
%
&\equiv \sum\limits_{i,j,k} 
   - \left\{
   \frac{e_{2u}\,e_{3u}} {e_{1u}}  A_u^{\,lm} \delta_{i+1/2} 
      \left[ \chi \right]
   \delta_{i+1/2} 
      \left[ 1 \right] 
   + \frac{e_{1v}\,e_{3v}} {e_{2v}}  A_v^{\,lm} \delta_{j+1/2} 
      \left[ \chi \right]
   \delta_{j+1/2} 
      \left[ 1 \right]
   \right\} \;
   \equiv 0 
   &&& \\ 
\end{flalign*}

% -------------------------------------------------------------------------------------------------------------
%       Dissipation of Horizontal Divergence Variance
% -------------------------------------------------------------------------------------------------------------
\subsection{Dissipation of Horizontal Divergence Variance}
\label{Apdx_C.3.5}

\begin{flalign*}
&\int\limits_D \chi \;\nabla_h \cdot 
   \left[    \nabla_h              \left( A^{\,lm}\;\chi                    \right)
           - \nabla_h   \times  \left( A^{\,lm}\;\zeta \;\textbf{k} \right)    \right]\;  dv
 = A^{\,lm}   \int\limits_D \chi \;\nabla_h \cdot \nabla_h \left( \chi \right)\;  dv    &&&\\ 
%
&\equiv A^{\,lm}  \sum\limits_{i,j,k}  \frac{1} {e_{1T}\,e_{2T}\,e_{3T}}  \chi 
   \left\{
      \delta_i  \left[   \frac{e_{2u}\,e_{3u}} {e_{1u}}  \delta_{i+1/2} \left[ \chi \right]   \right]
   + \delta_j  \left[   \frac{e_{1v}\,e_{3v}} {e_{2v}}   \delta_{j+1/2} \left[ \chi \right]   \right]
   \right\} \;   e_{1T}\,e_{2T}\,e_{3T}    &&&\\ 
%
\intertext{Using \eqref{DOM_di_adj}, it turns out to be:}
%
&\equiv - A^{\,lm} \sum\limits_{i,j,k}
   \left\{ 
   \left(  \frac{1} {e_{1u}}  \delta_{i+1/2}  \left[ \chi \right]  \right)^2  e_{1u}\,e_{2u}\,e_{3u}
+ \left(  \frac{1} {e_{2v}}  \delta_{j+1/2}  \left[ \chi \right]  \right)^2  e_{1v}\,e_{2v}\,e_{3v}
   \right\} \;    &&&\\
&\leq 0              &&&\\
\end{flalign*}

% ================================================================
% Conservation Properties on Vertical Momentum Physics
% ================================================================
\section{Conservation Properties on Vertical Momentum Physics}
\label{Apdx_C_4}


As for the lateral momentum physics, the continuous form of the vertical diffusion 
of momentum satisfies several integral constraints. The first two are associated 
with the conservation of momentum and the dissipation of horizontal kinetic energy:
\begin{align*}
 \int\limits_D  
    \frac{1} {e_3 }\; \frac{\partial } {\partial k}
   \left( 
   \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}
   \right)\;
   dv \qquad \quad &= \vec{\textbf{0}}
   \\
%
\intertext{and}
%
\int\limits_D 
   \textbf{U}_h \cdot 
   \frac{1} {e_3 }\; \frac{\partial } {\partial k}
   \left(
   \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}
   \right)\;
   dv \quad &\leq 0
   \\
\end{align*}
The first property is obvious. The second results from:

\begin{flalign*}
\int\limits_D 
   \textbf{U}_h \cdot 
   \frac{1} {e_3 }\; \frac{\partial } {\partial k}
   \left( 
   \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} 
   \right)\;dv
   &&&\\
\end{flalign*}
\begin{flalign*}
&\equiv \sum\limits_{i,j,k} 
   \left( 
     u\; \delta_k 
      \left[ 
      \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2} 
         \left[ u \right] 
      \right]\;
      e_{1u}\,e_{2u} 
   + v\;\delta_k 
      \left[
      \frac{A_v^{\,vm}} {e_{3vw}} \delta_{k+1/2} 
         \left[ v \right] 
      \right]\;
      e_{1v}\,e_{2v} 
   \right)
   &&&\\ 
%
\intertext{since the horizontal scale factor does not depend on $k$, it follows:}
%
&\equiv - \sum\limits_{i,j,k} 
   \left( 
      \frac{A_u^{\,vm}} {e_{3uw}}
      \left( 
      \delta_{k+1/2} 
         \left[ u \right] 
      \right)^2\;
      e_{1u}\,e_{2u} 
   + \frac{A_v^{\,vm}} {e_{3vw}} 
      \left( \delta_{k+1/2} 
         \left[ v \right] 
      \right)^2\;
      e_{1v}\,e_{2v}
   \right)
    \quad \leq 0
    &&&\\
\end{flalign*}

The vorticity is also conserved. Indeed:
\begin{flalign*}
\int \limits_D 
   \frac{1} {e_3 } \textbf{k} \cdot \nabla \times 
      \left( 
      \frac{1} {e_3 }\; \frac{\partial } {\partial k}
         \left( 
         \frac{A^{\,vm}} {e_3}\; \frac{\partial \textbf{U}_h } {\partial k} 
         \right)
      \right)\;
      dv
      &&&\\ 
\end{flalign*}
\begin{flalign*}
\equiv \sum\limits_{i,j,k}  \frac{1} {e_{3f}}\; \frac{1} {e_{1f}\,e_{2f}}
   \bigg\{    \biggr.   \quad
   \delta_{i+1/2} 
      &\left(
      \frac{e_{2v}} {e_{3v}} \delta_k 
         \left[
         \frac{1} {e_{3vw}} \delta_{k+1/2} 
            \left[ v \right] 
         \right]
      \right)
    &&\\
   \biggl. 
   - \delta_{j+1/2} 
      &\left(
      \frac{e_{1u}} {e_{3u}} \delta_k 
         \left[
         \frac{1} {e_{3uw}}\delta_{k+1/2} 
            \left[ u \right]
         \right]
      \right)
   \biggr\} \;
   e_{1f}\,e_{2f}\,e_{3f} \; \equiv 0
   && \\
\end{flalign*}
If the vertical diffusion coefficient is uniform over the whole domain, the 
enstrophy is dissipated, $i.e.$
\begin{flalign*}
\int\limits_D \zeta \, \textbf{k} \cdot \nabla \times 
   \left( 
   \frac{1} {e_3}\; \frac{\partial } {\partial k}
      \left( 
      \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} 
      \right) 
   \right)\;
   dv = 0
   &&&\\
\end{flalign*}
This property is only satisfied in $z$-coordinates:

\begin{flalign*}
\int\limits_D \zeta \, \textbf{k} \cdot \nabla \times 
   \left( 
   \frac{1} {e_3}\; \frac{\partial } {\partial k}
      \left( 
      \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} 
      \right)
   \right)\;
   dv
   &&& \\ 
\end{flalign*}
\begin{flalign*}
\equiv \sum\limits_{i,j,k} \zeta \;e_{3f} \;
   \biggl\{    \biggr.  \quad
   \delta_{i+1/2} 
      &\left( 
         \frac{e_{2v}} {e_{3v}} \delta_k 
         \left[ 
         \frac{A_v^{\,vm}} {e_{3vw}} \delta_{k+1/2} 
            \left[ v \right] 
         \right]
         \right)
         &&\\
   - \delta_{j+1/2} 
      &\biggl.
      \left(   
         \frac{e_{1u}} {e_{3u}} \delta_k 
         \left[ 
         \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2} 
            \left[ u \right]
          \right] 
         \right) 
   \biggr\} 
   &&\\ 
\end{flalign*}
\begin{flalign*}
\equiv \sum\limits_{i,j,k} \zeta \;e_{3f} 
   \biggl\{    \biggr.  \quad
   \frac{1} {e_{3v}} \delta_k 
      &\left[ 
      \frac{A_v^{\,vm}} {e_{3vw}} \delta_{k+1/2} 
         \left[ \delta_{i+1/2} 
            \left[ e_{2v}\,v \right]
          \right]
      \right]
      &&\\ 
    \biggl. 
   - \frac{1} {e_{3u}} \delta_k 
      &\left[ 
      \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2} 
         \left[ \delta_{j+1/2} 
            \left[ e_{1u}\,u \right]
          \right]
      \right]
   \biggr\} 
   &&\\ 
\end{flalign*}
Using the fact that the vertical diffusion coefficients are uniform, and that in 
$z$-coordinate, the vertical scale factors do not depend on $i$ and $j$ so 
that: $e_{3f} =e_{3u} =e_{3v} =e_{3T} $ and $e_{3w} =e_{3uw} =e_{3vw} $, 
it follows:
\begin{flalign*}
\equiv A^{\,vm} \sum\limits_{i,j,k} \zeta \;\delta_k 
   \left[ 
   \frac{1} {e_{3w}} \delta_{k+1/2} 
      \Bigl[ 
      \delta_{i+1/2} 
         \left[ e_{2v}\,v \right]
      - \delta_{j+1/ 2} 
         \left[ e_{1u}\,u \right]
       \Bigr] 
   \right] 
   &&&\\
\end{flalign*}
\begin{flalign*}
\equiv - A^{\,vm} \sum\limits_{i,j,k} \frac{1} {e_{3w}}
   \left( 
   \delta_{k+1/2} 
      \left[ \zeta  \right]
    \right)^2 \;
    e_{1f}\,e_{2f} 
    \; \leq 0
    &&&\\
\end{flalign*}
Similarly, the horizontal divergence is obviously conserved:

\begin{flalign*}
\int\limits_D \nabla \cdot 
   \left( 
   \frac{1} {e_3 }\; \frac{\partial } {\partial k}
      \left( 
      \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} 
      \right) 
   \right)\;
   dv = 0
   &&&\\
\end{flalign*}
and the square of the horizontal divergence decreases ($i.e.$ the horizontal 
divergence is dissipated) if the vertical diffusion coefficient is uniform over the 
whole domain:

\begin{flalign*}
\int\limits_D \chi \;\nabla \cdot 
   \left( 
   \frac{1} {e_3 }\; \frac{\partial } {\partial k}
      \left( 
      \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} 
      \right) 
   \right)\;
   dv = 0
   &&&\\
\end{flalign*}
This property is only satisfied in the $z$-coordinate:
\begin{flalign*}
\int\limits_D \chi \;\nabla \cdot 
   \left( 
   \frac{1} {e_3 }\; \frac{\partial } {\partial k}
      \left( 
      \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} 
      \right) 
   \right)\;
   dv
   &&&\\
\end{flalign*}
\begin{flalign*}
\equiv \sum\limits_{i,j,k} \frac{\chi } {e_{1T}\,e_{2T}}
   \biggl\{    \Biggr.  \quad
   \delta_{i+1/2} 
      &\left( 
         \frac{e_{2u}} {e_{3u}} \delta_k 
            \left[ 
         \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2} 
            \left[ u \right] 
         \right]
       \right)
       &&\\ 
   \Biggl. 
   + \delta_{j+1/2} 
      &\left( 
      \frac{e_{1v}} {e_{3v}} \delta_k 
         \left[ 
         \frac{A_v^{\,vm}} {e_{3vw}} \delta_{k+1/2} 
            \left[ v \right]
          \right]
      \right)
   \Biggr\} \;
   e_{1T}\,e_{2T}\,e_{3T} 
   &&\\ 
\end{flalign*}

\begin{flalign*}
\equiv A^{\,vm} \sum\limits_{i,j,k}  \chi \,
   \biggl\{ \biggr.  \quad
   \delta_{i+1/2}
      &\left( 
         \delta_k 
         \left[ 
         \frac{1} {e_{3uw}} \delta_{k+1/2} 
            \left[ e_{2u}\,u \right]
          \right]
         \right) 
         && \\ 
   \biggl. 
   + \delta_{j+1/2} 
      &\left( 
         \delta_k 
         \left[ 
         \frac{1} {e_{3vw}} \delta_{k+1/2} 
            \left[ e_{1v}\,v \right]
          \right] 
         \right)
   \biggr\} 
   && \\ 
\end{flalign*}

\begin{flalign*}
\equiv -A^{\,vm} \sum\limits_{i,j,k} 
\frac{\delta_{k+1/2} \left[ \chi \right]} {e_{3w}}\;
   \biggl\{ 
   \delta_{k+1/2} 
      \Bigl[
         \delta_{i+1/2} 
         \left[ e_{2u}\,u \right]
      + \delta_{j+1/2} 
         \left[ e_{1v}\,v \right]
      \Bigr]
   \biggr\} 
   &&&\\
\end{flalign*}

\begin{flalign*}
\equiv -A^{\,vm} \sum\limits_{i,j,k}
 \frac{1} {e_{3w}} 
\delta_{k+1/2} 
   \left[ \chi \right]\;
\delta_{k+1/2} 
   \left[ e_{1T}\,e_{2T} \;\chi \right]
&&&\\
\end{flalign*}

\begin{flalign*}
\equiv -A^{\,vm} \sum\limits_{i,j,k} 
\frac{e_{1T}\,e_{2T}} {e_{3w}}\;
   \left( 
   \delta_{k+1/2} 
      \left[ \chi \right]
   \right)^2
   \quad  \equiv 0
&&&\\
\end{flalign*}

% ================================================================
% Conservation Properties on Tracer Physics
% ================================================================
\section{Conservation Properties on Tracer Physics}
\label{Apdx_C.5}

The numerical schemes used for tracer subgridscale physics are written such 
that the heat and salt contents are conserved (equations in flux form, second 
order centered finite differences). Since a flux form is used to compute the 
temperature and salinity, the quadratic form of these quantities (i.e. their variance) 
globally tends to diminish. As for the advection term, there is generally no strict 
conservation of mass, even if in practice the mass is conserved to a very high 
accuracy. 

% -------------------------------------------------------------------------------------------------------------
%       Conservation of Tracers
% -------------------------------------------------------------------------------------------------------------
\subsection{Conservation of Tracers}
\label{Apdx_C.5.1}

constraint of conservation of tracers:
\begin{flalign*}
&\int\limits_D  \nabla  \cdot \left( A\;\nabla T \right)\;dv  &&&\\ 
\\
&\equiv \sum\limits_{i,j,k} 
   \biggl\{    \biggr.
   \delta_i 
      \left[ 
      A_u^{\,lT} \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} 
         \left[ T \right]
      \right]
   + \delta_j 
      \left[ 
      A_v^{\,lT} \frac{e_{1v}\,e_{3v}} {e_{2v}} \delta_{j+1/2} 
         \left[ T \right] 
      \right]
   &&\\  & \qquad \qquad \qquad \qquad \qquad \qquad \quad \;\;\;
   + \delta_k 
      \left[ 
      A_w^{\,vT} \frac{e_{1T}\,e_{2T}} {e_{3T}} \delta_{k+1/2} 
         \left[ T \right] 
      \right]
   \biggr\}   \quad  \equiv 0
   &&\\ 
\end{flalign*}

In fact, this property simply results from the flux form of the operator. 

% -------------------------------------------------------------------------------------------------------------
%       Dissipation of Tracer Variance
% -------------------------------------------------------------------------------------------------------------
\subsection{Dissipation of Tracer Variance}
\label{Apdx_C.5.2}

constraint on the dissipation of tracer variance:
\begin{flalign*}
\int\limits_D T\;\nabla & \cdot \left( A\;\nabla T \right)\;dv &&&\\ 
&\equiv   \sum\limits_{i,j,k} \; T
\biggl\{  \biggr.
     \delta_i \left[ A_u^{\,lT} \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} \left[T\right] \right]
& + \delta_j \left[ A_v^{\,lT} \frac{e_{1v} \,e_{3v}} {e_{2v}} \delta_{j+1/2} \left[T\right] \right] 
      \quad&& \\ 
 \biggl. 
&&+ \delta_k \left[A_w^{\,vT}\frac{e_{1T}\,e_{2T}} {e_{3T}}\delta_{k+1/2}\left[T\right]\right]
\biggr\} && 
\end{flalign*}
\begin{flalign*}
\equiv - \sum\limits_{i,j,k} 
   \biggl\{    \biggr.  \quad
   & A_u^{\,lT} 
      \left( 
      \frac{1} {e_{1u}} \delta_{i+1/2} 
         \left[ T \right]
      \right)^2
      e_{1u}\,e_{2u}\,e_{3u}
   && \\
   & + A_v^{\,lT} 
      \left( 
      \frac{1} {e_{2v}} \delta_{j+1/2} 
         \left[ T \right] 
      \right)^2
      e_{1v}\,e_{2v}\,e_{3v}
   && \\ 
   \biggl. 
   & + A_w^{\,vT} 
      \left( 
      \frac{1} {e_{3w}} \delta_{k+1/2} 
         \left[ T \right] 
      \right)^2
      e_{1w}\,e_{2w}\,e_{3w} 
   \biggr\} 
   \quad \leq 0 
   && \\ 
\end{flalign*}


%%%%  end of appendix in gm comment
%}