source: branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Annex_C.tex @ 6040

Last change on this file since 6040 was 6040, checked in by gm, 5 years ago

#1613: vvl by default : start to update the DOC for change in vvl, LDF and solvers

File size: 76.3 KB
Line 
1% ================================================================
2% Chapter Ñ Appendix C : Discrete Invariants of the Equations
3% ================================================================
4\chapter{Discrete Invariants of the Equations}
5\label{Apdx_C}
6\minitoc
7
8%%%  Appendix put in gmcomment as it has not been updated for z* and s coordinate
9%I'm writting this appendix. It will be available in a forthcoming release of the documentation
10
11%\gmcomment{
12
13\newpage
14$\ $\newline    % force a new ligne
15
16% ================================================================
17% Introduction / Notations
18% ================================================================
19\section{Introduction / Notations}
20\label{Apdx_C.0}
21
22Notation used in this appendix in the demonstations :
23
24fluxes at the faces of a $T$-box:
25\begin{equation*}
26U = e_{2u}\,e_{3u}\; u  \qquad  V = e_{1v}\,e_{3v}\; v  \qquad W = e_{1w}\,e_{2w}\; \omega     \\
27\end{equation*}
28
29volume of cells at $u$-, $v$-, and $T$-points:
30\begin{equation*}
31b_u = e_{1u}\,e_{2u}\,e_{3u}  \qquad  b_v = e_{1v}\,e_{2v}\,e_{3v}  \qquad b_t = e_{1t}\,e_{2t}\,e_{3t}     \\
32\end{equation*}
33
34partial derivative notation: $\partial_\bullet = \frac{\partial}{\partial \bullet}$
35
36$dv=e_1\,e_2\,e_3 \,di\,dj\,dk$  is the volume element, with only $e_3$ that depends on time.
37$D$ and $S$ are the ocean domain volume and surface, respectively.
38No wetting/drying is allow ($i.e.$ $\frac{\partial S}{\partial t} = 0$)
39Let $k_s$ and $k_b$ be the ocean surface and bottom, resp.
40($i.e.$ $s(k_s) = \eta$ and $s(k_b)=-H$, where $H$ is the bottom depth).
41\begin{flalign*}
42 z(k) = \eta - \int\limits_{\tilde{k}=k}^{\tilde{k}=k_s}  e_3(\tilde{k}) \;d\tilde{k}
43        = \eta - \int\limits_k^{k_s}  e_3 \;d\tilde{k} 
44\end{flalign*}
45
46Continuity equation with the above notation:
47\begin{equation*}
48\frac{1}{e_{3t}} \partial_t (e_{3t})+ \frac{1}{b_t}  \biggl\{  \delta_i [U] + \delta_j [V] + \delta_k [W] \biggr\} = 0
49\end{equation*}
50
51A quantity, $Q$ is conserved when its domain averaged time change is zero, that is when:
52\begin{equation*}
53\partial_t \left( \int_D{ Q\;dv } \right) =0
54\end{equation*}
55Noting that the coordinate system used ....  blah blah
56\begin{equation*}
57\partial_t \left( \int_D {Q\;dv} \right) =  \int_D { \partial_t \left( e_3 \, Q \right) e_1e_2\;di\,dj\,dk }
58                                                       =  \int_D { \frac{1}{e_3} \partial_t \left( e_3 \, Q \right) dv } =0
59\end{equation*}
60equation of evolution of $Q$ written as the time evolution of the vertical content of $Q$
61like for tracers, or momentum in flux form, the quadratic quantity $\frac{1}{2}Q^2$ is conserved when :
62\begin{flalign*}
63\partial_t \left(   \int_D{ \frac{1}{2} \,Q^2\;dv }   \right)
64=&  \int_D{ \frac{1}{2} \partial_t \left( \frac{1}{e_3}\left( e_3 \, Q \right)^2 \right) e_1e_2\;di\,dj\,dk } \\
65=&  \int_D {         Q   \;\partial_t\left( e_3 \, Q \right) e_1e_2\;di\,dj\,dk } 
66\int_D { \frac{1}{2} Q^2 \,\partial_t  (e_3) \;e_1e_2\;di\,dj\,dk } \\
67\end{flalign*}
68that is in a more compact form :
69\begin{flalign} \label{Eq_Q2_flux}
70\partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right)
71=&                   \int_D { \frac{Q}{e_3}  \partial_t \left( e_3 \, Q \right) dv } 
72  -   \frac{1}{2} \int_D {  \frac{Q^2}{e_3} \partial_t (e_3) \;dv }
73\end{flalign}
74equation of evolution of $Q$ written as the time evolution of $Q$
75like for momentum in vector invariant form, the quadratic quantity $\frac{1}{2}Q^2$ is conserved when :
76\begin{flalign*}
77\partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right)
78=&  \int_D { \frac{1}{2} \partial_t \left( e_3 \, Q^2 \right) \;e_1e_2\;di\,dj\,dk } \\
79=& \int_D {         Q      \partial_t Q  \;e_1e_2e_3\;di\,dj\,dk } 
80\int_D { \frac{1}{2} Q^2 \, \partial_t e_\;e_1e_2\;di\,dj\,dk } \\
81\end{flalign*}
82that is in a more compact form :
83\begin{flalign} \label{Eq_Q2_vect}
84\partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right)
85=& \int_D {         Q   \,\partial_t Q  \;dv } 
86+   \frac{1}{2} \int_D { \frac{1}{e_3} Q^2 \partial_t e_3 \;dv }
87\end{flalign}
88
89
90% ================================================================
91% Continuous Total energy Conservation
92% ================================================================
93\section{Continuous conservation}
94\label{Apdx_C.1}
95
96
97The discretization of pimitive equation in $s$-coordinate ($i.e.$ time and space varying
98vertical coordinate) must be chosen so that the discrete equation of the model satisfy
99integral constrains on energy and enstrophy.
100
101
102Let us first establish those constraint in the continuous world.
103The total energy ($i.e.$ kinetic plus potential energies) is conserved :
104\begin{flalign} \label{Eq_Tot_Energy}
105  \partial_t \left( \int_D \left( \frac{1}{2} {\textbf{U}_h}^2 +  \rho \, g \, z \right) \;dv \right)  = & 0
106\end{flalign}
107under the following assumptions: no dissipation, no forcing
108(wind, buoyancy flux, atmospheric pressure variations), mass
109conservation, and closed domain.
110
111This equation can be transformed to obtain several sub-equalities.
112The transformation for the advection term depends on whether
113the vector invariant form or the flux form is used for the momentum equation.
114Using \eqref{Eq_Q2_vect} and introducing \eqref{Apdx_A_dyn_vect} in \eqref{Eq_Tot_Energy} 
115for the former form and
116Using \eqref{Eq_Q2_flux} and introducing \eqref{Apdx_A_dyn_flux} in \eqref{Eq_Tot_Energy} 
117for the latter form  leads to:
118
119\begin{subequations} \label{E_tot}
120
121advection term (vector invariant form):
122\begin{equation} \label{E_tot_vect_vor}
123\int\limits_\zeta \; \left( \textbf{k} \times \textbf{U}_\right) \cdot \textbf{U}_\;  dv   = 0   \\
124\end{equation}
125%
126\begin{equation} \label{E_tot_vect_adv}
127   \int\limits_\textbf{U}_h \cdot \nabla_h \left( \frac{{\textbf{U}_h}^2}{2} \right)     dv
128+ \int\limits_\textbf{U}_h \cdot \nabla_z \textbf{U}_\;dv   
129\int\limits_D { \frac{{\textbf{U}_h}^2}{2} \frac{1}{e_3} \partial_t e_3 \;dv }   = 0   \\
130\end{equation}
131
132advection term (flux form):
133\begin{equation} \label{E_tot_flux_metric}
134\int\limits_\frac{1} {e_1 e_2 } \left( v \,\partial_i e_2 - u \,\partial_j e_\right)\; 
135 \left\textbf{k} \times \textbf{U}_\right) \cdot \textbf{U}_\;  dv   = 0   \\
136\end{equation}
137
138\begin{equation} \label{E_tot_flux_adv}
139   \int\limits_D \textbf{U}_h \cdot     \left(                 {{\begin{array} {*{20}c}
140\nabla \cdot \left( \textbf{U}\,u \right) \hfill \\
141\nabla \cdot \left( \textbf{U}\,v \right) \hfill \\       \end{array}} }           \right)   \;dv
142+   \frac{1}{2} \int\limits_D {  {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t e_3 \;dv } =\;\\
143\end{equation}
144
145coriolis term
146\begin{equation} \label{E_tot_cor}
147\int\limits_D  f   \; \left( \textbf{k} \times \textbf{U}_\right) \cdot \textbf{U}_\;  dv   = 0   \\
148\end{equation}
149
150pressure gradient:
151\begin{equation} \label{E_tot_pg}
152   - \int\limits_\left. \nabla p \right|_z \cdot \textbf{U}_h \;dv
153= - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
154   + \int\limits_D g\, \rho \; \partial_t z  \;dv   \\
155\end{equation}
156\end{subequations}
157
158where $\nabla_h = \left. \nabla \right|_k$ is the gradient along the $s$-surfaces.
159
160blah blah....
161$\ $\newline    % force a new ligne
162The prognostic ocean dynamics equation can be summarized as follows:
163\begin{equation*}
164\text{NXT} = \dbinom {\text{VOR} + \text{KEG} + \text {ZAD} }
165                  {\text{COR} + \text{ADV}                       }
166         + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF}
167\end{equation*}
168$\ $\newline    % force a new ligne
169
170Vector invariant form:
171\begin{subequations} \label{E_tot_vect}
172\begin{equation} \label{E_tot_vect_vor}
173\int\limits_D   \textbf{U}_h \cdot \text{VOR} \;dv   = 0   \\
174\end{equation}
175\begin{equation} \label{E_tot_vect_adv}
176   \int\limits_\textbf{U}_h \cdot \text{KEG}  \;dv
177+ \int\limits_\textbf{U}_h \cdot \text{ZAD}  \;dv   
178\int\limits_D { \frac{{\textbf{U}_h}^2}{2} \frac{1}{e_3} \partial_t e_3 \;dv }   = 0   \\
179\end{equation}
180\begin{equation} \label{E_tot_pg}
181   - \int\limits_\textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv
182= - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
183   + \int\limits_D g\, \rho \; \partial_t z  \;dv   \\
184\end{equation}
185\end{subequations}
186
187Flux form:
188\begin{subequations} \label{E_tot_flux}
189\begin{equation} \label{E_tot_flux_metric}
190\int\limits_\textbf{U}_h \cdot \text {COR} \;  dv   = 0   \\
191\end{equation}
192\begin{equation} \label{E_tot_flux_adv}
193   \int\limits_D \textbf{U}_h \cdot \text{ADV}   \;dv
194+   \frac{1}{2} \int\limits_D {  {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t e_\;dv } =\;\\
195\end{equation}
196\begin{equation} \label{E_tot_pg}
197   - \int\limits_\textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv
198= - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
199   + \int\limits_D g\, \rho \; \partial_t  z  \;dv   \\
200\end{equation}
201\end{subequations}
202
203
204$\ $\newline    % force a new ligne
205
206
207\eqref{E_tot_pg} is the balance between the conversion KE to PE and PE to KE.
208Indeed the left hand side of \eqref{E_tot_pg} can be transformed as follows:
209\begin{flalign*}
210\partial_\left( \int\limits_D { \rho \, g \, z  \;dv} \right)
211&= + \int\limits_D \frac{1}{e_3} \partial_t (e_3\,\rho) \;g\;z\;\;dv
212     +  \int\limits_D g\, \rho \; \partial_t z  \;dv   &&&\\
213&= - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
214     + \int\limits_D g\, \rho \; \partial_t z \;dv   &&&\\
215&= + \int\limits_\rho \,g \left( \textbf {U}_h \cdot \nabla_h z + \omega \frac{1}{e_3} \partial_k z \right\;dv
216     + \int\limits_D g\, \rho \; \partial_t z \;dv   &&&\\
217&= + \int\limits_\rho \,g \left( \omega + \partial_t z + \textbf {U}_h \cdot \nabla_h z  \right\;dv  &&&\\
218&=+  \int\limits_D g\, \rho \; w \; dv   &&&\\
219\end{flalign*}
220where the last equality is obtained by noting that the brackets is exactly the expression of $w$,
221the vertical velocity referenced to the fixe $z$-coordinate system (see \eqref{Apdx_A_w_s}).
222 
223The left hand side of \eqref{E_tot_pg} can be transformed as follows:
224\begin{flalign*}
225- \int\limits_\left. \nabla p \right|_z & \cdot \textbf{U}_h \;dv 
226= - \int\limits_\left( \nabla_h p + \rho \, g \nabla_h z \right) \cdot \textbf{U}_h \;dv   &&&\\
227\allowdisplaybreaks
228&= - \int\limits_\nabla_h  p \cdot \textbf{U}_h \;dv   - \int\limits_\rho \, g \, \nabla_h z \cdot \textbf{U}_h \;dv   &&&\\
229\allowdisplaybreaks
230&= +\int\limits_D p \,\nabla_h \cdot \textbf{U}_h \;dv   + \int\limits_\rho \, g \left( \omega - w + \partial_t z \right) \;dv   &&&\\
231\allowdisplaybreaks
232&= -\int\limits_D p \left( \frac{1}{e_3} \partial_t e_3 + \frac{1}{e_3} \partial_k \omega  \right) \;dv
233    +\int\limits_\rho \, g \left( \omega - w + \partial_t z \right) \;dv   &&&\\
234\allowdisplaybreaks
235&= -\int\limits_D \frac{p}{e_3} \partial_t e_\;dv   
236     +\int\limits_D \frac{1}{e_3} \partial_k p\; \omega \;dv
237    +\int\limits_\rho \, g \left( \omega - w + \partial_t z \right) \;dv   &&&\\
238&= -\int\limits_D \frac{p}{e_3} \partial_t e_\;dv   
239     -\int\limits_D \rho \, g \, \omega \;dv
240    +\int\limits_\rho \, g \left( \omega - w + \partial_t z \right) \;dv   &&&\\
241&= - \int\limits_D \frac{p}{e_3} \partial_t e_3 \; \;dv   
242     - \int\limits_\rho \, g \, w \;dv
243     + \int\limits_D   \rho \, g \, \partial_t z \;dv   &&&\\
244\allowdisplaybreaks
245\intertext{introducing the hydrostatic balance $\partial_k p=-\rho \,g\,e_3$ in the last term,
246it becomes:}
247&= - \int\limits_D \frac{p}{e_3} \partial_t e_3 \;dv   
248     - \int\limits_\rho \, g \, w \;dv
249     - \int\limits_\frac{1}{e_3} \partial_k p\, \partial_t z \;dv   &&&\\
250&= - \int\limits_D \frac{p}{e_3} \partial_t e_3 \;dv   
251     - \int\limits_\rho \, g \, w \;dv
252     + \int\limits_D \,\frac{p}{e_3}\partial_t ( \partial_k z )  dv   &&&\\
253%
254&= - \int\limits_\rho \, g \, w \;dv   &&&\\
255\end{flalign*}
256
257
258%gm comment
259\gmcomment{
260%
261The last equality comes from the following equation,
262\begin{flalign*}
263\int\limits_D p \frac{1}{e_3} \frac{\partial e_3}{\partial t}\; \;dv   
264     = \int\limits_D   \rho \, g \, \frac{\partial z }{\partial t} \;dv \quad\\ 
265\end{flalign*}
266that can be demonstrated as follows:
267
268\begin{flalign*}
269\int\limits_D   \rho \, g \, \frac{\partial z }{\partial t} \;dv   
270&= \int\limits_D    \rho \, g \, \frac{\partial \eta}{\partial t} \;dv   
271  -  \int\limits_D    \rho \, g \, \frac{\partial}{\partial t} \left\int\limits_k^{k_s}  e_3 \;d\tilde{k} \right) \;dv   &&&\\
272&= \int\limits_D    \rho \, g \, \frac{\partial \eta}{\partial t} \;dv   
273  -  \int\limits_D    \rho \, g    \left\int\limits_k^{k_s}  \frac{\partial e_3}{\partial t} \;d\tilde{k} \right) \;dv   &&&\\
274%
275\allowdisplaybreaks
276\intertext{The second term of the right hand side can be transformed by applying the integration by part rule:
277$\left[ a\,b \right]_{k_b}^{k_s} = \int_{k_b}^{k_s}  a\,\frac{\partial b}{\partial k}       \;dk
278                                               + \int_{k_b}^{k_s}      \frac{\partial a}{\partial k} \,b \;dk $
279to the following function: $a=  \int_k^{k_s}  \frac{\partial e_3}{\partial t} \;d\tilde{k}$ 
280and $b=  \int_k^{k_s}  \rho \, e_3 \;d\tilde{k}$
281(note that $\frac{\partial}{\partial k} \left(  \int_k^{k_s}  a \;d\tilde{k}  \right) = - a$ as $k$ is the lower bound of the integral).
282This leads to:  }
283\end{flalign*}
284\begin{flalign*}
285&\left[ \int\limits_{k}^{k_s}  \frac{\partial e_3}{\partial t} \,dk \cdot \int\limits_{k}^{k_s}  \rho \, e_3 \,dk   \right]_{k_b}^{k_s}
286=-\int\limits_{k_b}^{k_s} \left\int\limits_k^{k_s}  \frac{\partial e_3}{\partial t} \;d\tilde{k} \right\rho \,e_3 \;dk
287  -\int\limits_{k_b}^{k_s}  \frac{\partial e_3}{\partial t}  \left\int\limits_k^{k_s}  \rho \, e_3 \;d\tilde{k} \right)   dk
288&&&\\
289\allowdisplaybreaks
290\intertext{Noting that $\frac{\partial \eta}{\partial t}
291      = \frac{\partial}{\partial t}  \left( \int_{k_b}^{k_s}   e_3  \;d\tilde{k}  \right)
292      = \int_{k_b}^{k_s}  \frac{\partial e_3}{\partial t} \;d\tilde{k}$
293and
294      $p(k) = \int_k^{k_s}  \rho \,g \, e_3 \;d\tilde{k} $,
295but also that $\frac{\partial \eta}{\partial t}$ does not depends on $k$, it comes:
296}
297& - \int\limits_{k_b}^{k_s}  \rho \, \frac{\partial \eta}{\partial t} \, e_3 \;dk
298= - \int\limits_{k_b}^{k_s} \left\int\limits_k^{k_s}  \frac{\partial e_3}{\partial t} \;d\tilde{k} \right)   \, \rho \, g   e_3\;dk
299   - \int\limits_{k_b}^{k_s}  \frac{\partial e_3}{\partial t} \frac{p}{g}         \;dk       &&&\\
300\end{flalign*}
301Mutliplying by $g$ and integrating over the $(i,j)$ domain it becomes:
302\begin{flalign*}
303\int\limits_\rho \, g \, \left\int\limits_k^{k_s}  \frac{\partial e_3}{\partial t} \;d\tilde{k} \right)    \;dv
304\int\limits_\rho \, g \, \frac{\partial \eta}{\partial t} dv
305  - \int\limits_\frac{p}{e_3}\frac{\partial e_3}{\partial t}         \;dv
306\end{flalign*}
307Using this property, we therefore have:
308\begin{flalign*}
309\int\limits_D   \rho \, g \, \frac{\partial z }{\partial t} \;dv   
310&= \int\limits_D    \rho \, g \, \frac{\partial \eta}{\partial t}   \;dv   
311  - \left\int\limits_\rho \, g \, \frac{\partial \eta}{\partial t} dv
312           - \int\limits_\frac{p}{e_3}\frac{\partial e_3}{\partial t}   \;dv  \right)    &&&\\
313%
314&=\int\limits_D \frac{p}{e_3} \frac{\partial (e_3\,\rho)}{\partial t}\; \;dv     
315\end{flalign*}
316% end gm comment
317}
318%
319
320
321% ================================================================
322% Discrete Total energy Conservation : vector invariant form
323% ================================================================
324\section{Discrete total energy conservation : vector invariant form}
325\label{Apdx_C.1}
326
327% -------------------------------------------------------------------------------------------------------------
328%       Total energy conservation
329% -------------------------------------------------------------------------------------------------------------
330\subsection{Total energy conservation}
331\label{Apdx_C_KE+PE}
332
333The discrete form of the total energy conservation, \eqref{Eq_Tot_Energy}, is given by:
334\begin{flalign*}
335\partial_\left\sum\limits_{i,j,k} \biggl\{ \frac{u^2}{2} \,b_u + \frac{v^2}{2}\, b_v +  \rho \, g \, z_t \,b_\biggr\} \right) &=0  \\
336\end{flalign*}
337which in vector invariant forms, it leads to:
338\begin{equation} \label{KE+PE_vect_discrete}   \begin{split}
339                        \sum\limits_{i,j,k} \biggl\{   u\,                        \partial_t u         \;b_u
340                                                              + v\,                        \partial_t v          \;b_\biggr\}
341  + \frac{1}{2} \sum\limits_{i,j,k} \biggl\{  \frac{u^2}{e_{3u}}\partial_t e_{3u} \;b_u
342                                                             +  \frac{v^2}{e_{3v}}\partial_t e_{3v} \;b_v   \biggr\}      \\
343= - \sum\limits_{i,j,k} \biggl\{ \frac{1}{e_{3t}}\partial_t (e_{3t} \rho) \, g \, z_t \;b_\biggr\} 
344     - \sum\limits_{i,j,k} \biggl\{ \rho \,g\,\partial_t (z_t) \,b_\biggr\}                               
345\end{split} \end{equation}
346
347Substituting the discrete expression of the time derivative of the velocity either in vector invariant,
348leads to the discrete equivalent of the four equations \eqref{E_tot_flux}.
349
350% -------------------------------------------------------------------------------------------------------------
351%       Vorticity term (coriolis + vorticity part of the advection)
352% -------------------------------------------------------------------------------------------------------------
353\subsection{Vorticity term (coriolis + vorticity part of the advection)}
354\label{Apdx_C_vor}
355
356Let $q$, located at $f$-points, be either the relative ($q=\zeta / e_{3f}$), or 
357the planetary ($q=f/e_{3f}$), or the total potential vorticity ($q=(\zeta +f) /e_{3f}$).
358Two discretisation of the vorticity term (ENE and EEN) allows the conservation of
359the kinetic energy.
360% -------------------------------------------------------------------------------------------------------------
361%       Vorticity Term with ENE scheme
362% -------------------------------------------------------------------------------------------------------------
363\subsubsection{Vorticity Term with ENE scheme (\np{ln\_dynvor\_ene}=.true.)}
364\label{Apdx_C_vorENE} 
365
366For the ENE scheme, the two components of the vorticity term are given by :
367\begin{equation*}
368- e_3 \, q \;{\textbf{k}}\times {\textbf {U}}_h    \equiv 
369   \left( {{  \begin{array} {*{20}c}
370      + \frac{1} {e_{1u}} \; 
371      \overline {\, q \ \overline {\left( e_{1v}\,e_{3v}\,v \right)}^{\,i+1/2}} ^{\,j}        \hfill \\
372      - \frac{1} {e_{2v}} \; 
373      \overline {\, q \ \overline {\left( e_{2u}\,e_{3u}\,u \right)}^{\,j+1/2}} ^{\,i}       \hfill \\
374   \end{array}} }    \right)
375\end{equation*}
376
377This formulation does not conserve the enstrophy but it does conserve the
378total kinetic energy. Indeed, the kinetic energy tendency associated to the
379vorticity term and averaged over the ocean domain can be transformed as
380follows:
381\begin{flalign*}
382&\int\limits_D -  \left(  e_3 \, q \;\textbf{k} \times \textbf{U}_\right) \cdot \textbf{U}_\;  dv &&  \\
383& \qquad \qquad {\begin{array}{*{20}l} 
384&\equiv \sum\limits_{i,j,k}   \biggl\{   
385     \frac{1} {e_{1u}} \overline { \,q\ \overline{ V }^{\,i+1/2}} ^{\,j} \, u \; b_u
386   - \frac{1} {e_{2v}}\overline { \, q\ \overline{ U }^{\,j+1/2}} ^{\,i} \, v \; b_v \; \biggr\}    \\ 
387&\equiv  \sum\limits_{i,j,k}  \biggl\{   
388     \overline { \,q\ \overline{ V }^{\,i+1/2}}^{\,j} \; U
389   - \overline { \,q\ \overline{ U }^{\,j+1/2}}^{\,i} \; V  \; \biggr\}     \\
390&\equiv \sum\limits_{i,j,k} q \  \biggl\{  \overline{ V }^{\,i+1/2}\; \overline{ U }^{\,j+1/2}       
391                                              - \overline{ U }^{\,j+1/2}\; \overline{ V }^{\,i+1/2}         \biggr\}  \quad  \equiv 0
392\end{array} }     
393\end{flalign*}
394In other words, the domain averaged kinetic energy does not change due to the vorticity term.
395
396
397% -------------------------------------------------------------------------------------------------------------
398%       Vorticity Term with EEN scheme
399% -------------------------------------------------------------------------------------------------------------
400\subsubsection{Vorticity Term with EEN scheme (\np{ln\_dynvor\_een}=.true.)}
401\label{Apdx_C_vorEEN} 
402
403With the EEN scheme, the vorticity terms are represented as:
404\begin{equation} \label{Eq_dynvor_een}
405\left\{ {    \begin{aligned}
406 +q\,e_3 \, v  &\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}} 
407                         {^{i+1/2-i_p}_j}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{1v} e_{3v} \ \right)^{i+i_p-1/2}_{j+j_p}   \\
408 - q\,e_3 \, u     &\equiv -\frac{1}{e_{2v} }    \sum_{\substack{i_p,\,k_p}} 
409                         {^i_{j+1/2-j_p}}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{2u} e_{3u} \ \right)^{i+i_p}_{j+j_p-1/2}   \\
410\end{aligned}   } \right.
411\end{equation} 
412where the indices $i_p$ and $j_p$ take the following value:
413$i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$,
414and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by:
415\begin{equation} \label{Q_triads}
416_i^j \mathbb{Q}^{i_p}_{j_p}
417= \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right)
418\end{equation}
419
420This formulation does conserve the total kinetic energy. Indeed,
421\begin{flalign*}
422&\int\limits_D - \textbf{U}_h \cdot   \left\zeta \;\textbf{k} \times \textbf{U}_\right\;  dv &&  \\
423\equiv \sum\limits_{i,j,k} &  \biggl\{ 
424       \left\sum_{\substack{i_p,\,k_p}} 
425                         {^{i+1/2-i_p}_j}\mathbb{Q}^{i_p}_{j_p} \; V^{i+1/2-i_p}_{j+j_p} \right] U^{i+1/2}_{j}    %   &&\\
426     - \left\sum_{\substack{i_p,\,k_p}} 
427                         {^i_{j+1/2-j_p}}\mathbb{Q}^{i_p}_{j_p} \; U^{i+i_p}_{j+1/2-j_p}  \right] V^{i}_{j+1/2}    \biggr\}     && \\
428\\
429\equiv \sum\limits_{i,j,k} &  \sum_{\substack{i_p,\,k_p}} \biggl\{  \ \  
430                         {^{i+1/2-i_p}_j}\mathbb{Q}^{i_p}_{j_p} \; V^{i+1/2-i_p}_{j+j_p}  \, U^{i+1/2}_{j}     %  &&\\
431                       - {^i_{j+1/2-j_p}}\mathbb{Q}^{i_p}_{j_p} \; U^{i+i_p}_{j+1/2-j_p} \, V^{i}_{j+1/2}     \ \;     \biggr\}     &&  \\
432%
433\allowdisplaybreaks
434\intertext{ Expending the summation on $i_p$ and $k_p$, it becomes:}
435%
436\equiv \sum\limits_{i,j,k} & \biggl\{  \ \  
437                    {^{i+1}_j     }\mathbb{Q}^{-1/2}_{+1/2} \;V^{i+1}_{j+1/2} \; U^{\,i+1/2}_{j}     
438                -  {^i_{j}\quad}\mathbb{Q}^{-1/2}_{+1/2} \; U^{i-1/2}_{j}    \; V^{\,i}_{j+1/2}         &&  \\
439        &       + {^{i+1}_j     }\mathbb{Q}^{-1/2}_{-1/2} \; V^{i+1}_{j-1/2} \; U^{\,i+1/2}_{j}       
440                  - {^i_{j+1}     }\mathbb{Q}^{-1/2}_{-1/2} \; U^{i-1/2}_{j+1} \; V^{\,i}_{j+1/2}        \biggr.     &&  \\
441        &       + {^{i}_j\quad}\mathbb{Q}^{+1/2}_{+1/2} \; V^{i}_{j+1/2}   \; U^{\,i+1/2}_{j}     
442                  - {^i_{j}\quad}\mathbb{Q}^{+1/2}_{+1/2} \; U^{i+1/2}_{j}   \; V^{\,i}_{j+1/2}          \biggr.        &&  \\
443        &       + {^{i}_j\quad}\mathbb{Q}^{+1/2}_{-1/2} \; V^{i}_{j-1/2}     \; U^{\,i+1/2}_{j}       
444                 -  {^i_{j+1}     }\mathbb{Q}^{+1/2}_{-1/2} \; U^{i+1/2}_{j+1}\; V^{\,i}_{j+1/2}  \ \;     \biggr\}     &&  \\
445%
446\allowdisplaybreaks
447\intertext{The summation is done over all $i$ and $j$ indices, it is therefore possible to introduce
448a shift of $-1$ either in $i$ or $j$ direction in some of the term of the summation (first term of the
449first and second lines, second term of the second and fourth lines). By doning so, we can regroup
450all the terms of the summation by triad at a ($i$,$j$) point. In other words, we regroup all the terms
451in the neighbourhood  that contain a triad at the same ($i$,$j$) indices. It becomes: }
452\allowdisplaybreaks
453%
454\equiv \sum\limits_{i,j,k} & \biggl\{  \ \  
455             {^{i}_j}\mathbb{Q}^{-1/2}_{+1/2}  \left[  V^{i}_{j+1/2}\, U^{\,i-1/2}_{j}   
456                                                                       -  U^{i-1/2}_{j} \, V^{\,i}_{j+1/2}      \right]    &&  \\
457 &       + {^{i}_j}\mathbb{Q}^{-1/2}_{-1/2}  \left[  V^{i}_{j-1/2} \, U^{\,i-1/2}_{j}       
458                                                                     -    U^{i-1/2}_{j} \, V^{\,i}_{j-1/2}      \right]    \biggr.   &&  \\
459 &      + {^{i}_j}\mathbb{Q}^{+1/2}_{+1/2}  \left[  V^{i}_{j+1/2} \, U^{\,i+1/2}_{j}     
460                                                                     -    U^{i+1/2}_{j} \, V^{\,i}_{j+1/2}     \right\biggr&&  \\
461 &     + {^{i}_j}\mathbb{Q}^{+1/2}_{-1/2}  \left[   V^{i}_{j-1/2} \, U^{\,i+1/2}_{j}                                                                     
462                                                                    -    U^{i+1/2}_{j-1} \, V^{\,i}_{j-1/2}  \right\ \;   \biggr\}   \qquad
463\equiv \ 0   &&
464\end{flalign*}
465
466
467% -------------------------------------------------------------------------------------------------------------
468%       Gradient of Kinetic Energy / Vertical Advection
469% -------------------------------------------------------------------------------------------------------------
470\subsubsection{Gradient of Kinetic Energy / Vertical Advection}
471\label{Apdx_C_zad} 
472
473The change of Kinetic Energy (KE) due to the vertical advection is exactly
474balanced by the change of KE due to the horizontal gradient of KE~:
475\begin{equation*}
476    \int_D \textbf{U}_h \cdot \frac{1}{e_3 } \omega \partial_k \textbf{U}_h \;dv
477=  -   \int_D \textbf{U}_h \cdot \nabla_h \left( \frac{1}{2}\;{\textbf{U}_h}^2 \right)\;dv
478   +   \frac{1}{2} \int_D {  \frac{{\textbf{U}_h}^2}{e_3} \partial_t ( e_3) \;dv }  \\
479\end{equation*}
480Indeed, using successively \eqref{DOM_di_adj} ($i.e.$ the skew symmetry
481property of the $\delta$ operator) and the continuity equation, then
482\eqref{DOM_di_adj} again, then the commutativity of operators
483$\overline {\,\cdot \,}$ and $\delta$, and finally \eqref{DOM_mi_adj} 
484($i.e.$ the symmetry property of the $\overline {\,\cdot \,}$ operator)
485applied in the horizontal and vertical directions, it becomes:
486\begin{flalign*}
487& - \int_D \textbf{U}_h \cdot \text{KEG}\;dv   
488= - \int_D \textbf{U}_h \cdot \nabla_h \left( \frac{1}{2}\;{\textbf{U}_h}^2 \right)\;dv    &&&\\
489%
490\equiv  & -  \sum\limits_{i,j,k} \frac{1}{2}  \biggl\{   
491   \frac{1} {e_{1u}}  \delta_{i+1/2}   \left[   \overline {u^2}^{\,i} + \overline {v^2}^{\,j}   \right]  u \ b_u
492     + \frac{1} {e_{2v}}  \delta_{j+1/2}   \left[   \overline {u^2}^{\,i} + \overline {v^2}^{\,j}   \right]  v \ b_v   \biggr\}     &&&  \\ 
493%
494\equiv & + \sum\limits_{i,j,k} \frac{1}{2}  \left(   \overline {u^2}^{\,i} + \overline {v^2}^{\,j}   \right)\;
495                                                              \biggl\{ \delta_{i} \left[  U   \right] +  \delta_{j} \left[  V  \right]    \biggr\}       &&&  \\ 
496\allowdisplaybreaks
497%
498\equiv   & - \sum\limits_{i,j,k}  \frac{1}{2}
499   \left(       \overline {u^2}^{\,i} + \overline {v^2}^{\,j}   \right\; 
500   \biggl\{   \frac{b_t}{e_{3t}} \partial_t (e_{3t})  +  \delta_k \left[  W   \right]    \biggr\}    &&&\\
501\allowdisplaybreaks
502%
503\equiv & +  \sum\limits_{i,j,k} \frac{1}{2} \delta_{k+1/2}   \left[ \overline{ u^2}^{\,i} + \overline{ v^2}^{\,j}   \right] \;  W
504          -  \sum\limits_{i,j,k} \frac{1}{2} \left(   \overline {u^2}^{\,i} + \overline {v^2}^{\,j}   \right) \;\partial_t b_t   &&& \\
505\allowdisplaybreaks
506%
507\equiv   & + \sum\limits_{i,j,k} \frac{1} {2} \left(    \overline{\delta_{k+1/2} \left[ u^2 \right]}^{\,i} 
508                                                                   + \overline{\delta_{k+1/2} \left[ v^2 \right]}^{\,j}    \right) \; W     
509          -  \sum\limits_{i,j,k}  \left\frac{u^2}{2}\,\partial_t \overline{b_t}^{\,{i+1/2}}
510                                                 + \frac{v^2}{2}\,\partial_t \overline{b_t}^{\,{j+1/2}}   \right)    &&& \\
511\allowdisplaybreaks
512\intertext{Assuming that $b_u= \overline{b_t}^{\,i+1/2}$ and $b_v= \overline{b_t}^{\,j+1/2}$, or at least that the time
513derivative of these two equations is satisfied, it becomes:}
514%
515\equiv &     \sum\limits_{i,j,k} \frac{1} {2}
516   \biggl\{ \; \overline{W}^{\,i+1/2}\;\delta_{k+1/2} \left[ u^2 \right]   
517               + \overline{W}^{\,j+1/2}\;\delta_{k+1/2} \left[ v^2 \right]  \;  \biggr\}   
518          -  \sum\limits_{i,j,k}  \left\frac{u^2}{2}\,\partial_t b_u
519                                                + \frac{v^2}{2}\,\partial_t b_v   \right)    &&& \\
520\allowdisplaybreaks
521
522\equiv &     \sum\limits_{i,j,k} 
523   \biggl\{ \; \overline{W}^{\,i+1/2}\; \overline {u}^{\,k+1/2}\; \delta_{k+1/2}[ u ]   
524               + \overline{W}^{\,j+1/2}\; \overline {v}^{\,k+1/2}\; \delta_{k+1/2}[ v ]  \;  \biggr\} 
525          -  \sum\limits_{i,j,k}  \left\frac{u^2}{2}\,\partial_t b_u
526                                                + \frac{v^2}{2}\,\partial_t b_v   \right)    &&& \\
527%
528\allowdisplaybreaks
529\equiv  &  \sum\limits_{i,j,k} 
530   \biggl\{ \; \frac{1} {b_u } \; \overline { \overline{W}^{\,i+1/2}\,\delta_{k+1/2}  \left[ u \right] }^{\,k} \;u\;b_
531                    + \frac{1} {b_v } \; \overline { \overline{W}^{\,j+1/2} \delta_{k+1/2}  \left[ v \right]  }^{\,k} \;v\;b_\; \biggr\} 
532          -  \sum\limits_{i,j,k}  \left\frac{u^2}{2}\,\partial_t b_u
533                                                + \frac{v^2}{2}\,\partial_t b_v   \right)    &&& \\
534%
535\intertext{The first term provides the discrete expression for the vertical advection of momentum (ZAD),
536while the second term corresponds exactly to \eqref{KE+PE_vect_discrete}, therefore:}
537\equiv&                   \int\limits_D \textbf{U}_h \cdot \text{ZAD} \;dv 
538           + \frac{1}{2} \int_D { {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t  (e_3)  \;dv }    &&&\\
539\equiv&                   \int\limits_D \textbf{U}_h \cdot w \partial_k \textbf{U}_h \;dv 
540           + \frac{1}{2} \int_D { {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t  (e_3)  \;dv }    &&&\\
541\end{flalign*}
542
543There is two main points here. First, the satisfaction of this property links the choice of
544the discrete formulation of the vertical advection and of the horizontal gradient
545of KE. Choosing one imposes the other. For example KE can also be discretized
546as $1/2\,({\overline u^{\,i}}^2 + {\overline v^{\,j}}^2)$. This leads to the following
547expression for the vertical advection:
548\begin{equation*}
549\frac{1} {e_3 }\; \omega\; \partial_k \textbf{U}_h
550\equiv \left( {{\begin{array} {*{20}c}
551\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} \;  \overline{\overline {e_{1t}\,e_{2t} \,\omega\;\delta_{k+1/2} 
552\left[ \overline u^{\,i+1/2} \right]}}^{\,i+1/2,k}  \hfill \\
553\frac{1} {e_{1v}\,e_{2v}\,e_{3v}} \;   \overline{\overline {e_{1t}\,e_{2t} \,\omega \;\delta_{k+1/2}
554\left[ \overline v^{\,j+1/2} \right]}}^{\,j+1/2,k} \hfill \\
555\end{array}} } \right)
556\end{equation*}
557a formulation that requires an additional horizontal mean in contrast with
558the one used in NEMO. Nine velocity points have to be used instead of 3.
559This is the reason why it has not been chosen.
560
561Second, as soon as the chosen $s$-coordinate depends on time, an extra constraint
562arises on the time derivative of the volume at $u$- and $v$-points:
563\begin{flalign*}
564e_{1u}\,e_{2u}\,\partial_t (e_{3u}) =\overline{ e_{1t}\,e_{2t}\;\partial_t (e_{3t}) }^{\,i+1/2}    \\
565e_{1v}\,e_{2v}\,\partial_t (e_{3v})  =\overline{ e_{1t}\,e_{2t}\;\partial_t (e_{3t}) }^{\,j+1/2}
566\end{flalign*}
567which is (over-)satified by defining the vertical scale factor as follows:
568\begin{flalign} \label{e3u-e3v}
569e_{3u} = \frac{1}{e_{1u}\,e_{2u}}\;\overline{ e_{1t}^{ }\,e_{2t}^{ }\,e_{3t}^{ } }^{\,i+1/2}    \\
570e_{3v} = \frac{1}{e_{1v}\,e_{2v}}\;\overline{ e_{1t}^{ }\,e_{2t}^{ }\,e_{3t}^{ } }^{\,j+1/2} 
571\end{flalign}
572
573Blah blah required on the the step representation of bottom topography.....
574
575
576% -------------------------------------------------------------------------------------------------------------
577%       Pressure Gradient Term
578% -------------------------------------------------------------------------------------------------------------
579\subsection{Pressure Gradient Term}
580\label{Apdx_C.1.4}
581
582\gmcomment{
583A pressure gradient has no contribution to the evolution of the vorticity as the
584curl of a gradient is zero. In the $z$-coordinate, this property is satisfied locally
585on a C-grid with 2nd order finite differences (property \eqref{Eq_DOM_curl_grad}).
586}
587
588When the equation of state is linear ($i.e.$ when an advection-diffusion equation
589for density can be derived from those of temperature and salinity) the change of
590KE due to the work of pressure forces is balanced by the change of potential
591energy due to buoyancy forces:
592\begin{equation*}
593- \int_\left. \nabla p \right|_z \cdot \textbf{U}_h \;dv
594= - \int_D \nabla \cdot \left( \rho \,\textbf {U} \right) \,g\,z \;dv
595  + \int_D g\, \rho \; \partial_t (z)  \;dv
596\end{equation*}
597
598This property can be satisfied in a discrete sense for both $z$- and $s$-coordinates.
599Indeed, defining the depth of a $T$-point, $z_t$, as the sum of the vertical scale
600factors at $w$-points starting from the surface, the work of pressure forces can be
601written as:
602\begin{flalign*}
603&- \int_\left. \nabla p \right|_z \cdot \textbf{U}_h \;dv   
604\equiv \sum\limits_{i,j,k} \biggl\{ \;  - \frac{1} {e_{1u}}   \Bigl(
605\delta_{i+1/2} [p_t] - g\;\overline \rho^{\,i+1/2}\;\delta_{i+1/2} [z_t]     \Bigr\; u\;b_
606   &&  \\ & \qquad \qquad  \qquad \qquad  \qquad \quad \ \,
607                        - \frac{1} {e_{2v}}    \Bigl(
608\delta_{j+1/2} [p_t] - g\;\overline \rho^{\,j+1/2}\delta_{j+1/2} [z_t]      \Bigr\; v\;b_v \;  \biggr\}   && \\ 
609%
610\allowdisplaybreaks
611\intertext{Using successively  \eqref{DOM_di_adj}, $i.e.$ the skew symmetry property of
612the $\delta$ operator, \eqref{Eq_wzv}, the continuity equation, \eqref{Eq_dynhpg_sco},
613the hydrostatic equation in the $s$-coordinate, and $\delta_{k+1/2} \left[ z_t \right] \equiv e_{3w} $,
614which comes from the definition of $z_t$, it becomes: }
615\allowdisplaybreaks
616%
617\equiv& +  \sum\limits_{i,j,k}   g  \biggl\{ 
618      \overline\rho^{\,i+1/2}\,U\,\delta_{i+1/2}[z_t]   
619   +     \overline\rho^{\,j+1/2}\,V\,\delta_{j+1/2}[z_t]     
620   +\Bigl\delta_i[U] + \delta_j [V]  \Bigr)\;\frac{p_t}{g} \biggr\}  &&\\
621%
622\equiv& +  \sum\limits_{i,j,k}   g   \biggl\{ 
623      \overline\rho^{\,i+1/2}\,U\,\delta_{i+1/2}[z_t]   
624   +     \overline\rho^{\,j+1/2}\,V\,\delta_{j+1/2}[z_t]     
625    -       \left(   \frac{b_t}{e_{3t}} \partial_t (e_{3t})  +  \delta_k \left[ W \right]    \right) \frac{p_t}{g}    \biggr\}   &&&\\ 
626%
627\equiv& +  \sum\limits_{i,j,k}  g   \biggl\{ 
628      \overline\rho^{\,i+1/2}\,U\,\delta_{i+1/2}[z_t]   
629   +     \overline\rho^{\,j+1/2}\,V\,\delta_{j+1/2}[z_t]     
630   +  \frac{W}{g}\;\delta_{k+1/2} [p_t]
631    -        \frac{p_t}{g}\,\partial_t b_t    \biggr\}    &&&\\
632%
633\equiv& +  \sum\limits_{i,j,k}  g   \biggl\{ 
634      \overline\rho^{\,i+1/2}\,U\,\delta_{i+1/2}[z_t]   
635   +     \overline\rho^{\,j+1/2}\,V\,\delta_{j+1/2}[z_t]     
636   -  W\;e_{3w} \overline \rho^{\,k+1/2}         
637    -        \frac{p_t}{g}\,\partial_t b_t    \biggr\}    &&&\\
638%
639\equiv& +  \sum\limits_{i,j,k}    g   \biggl\{ 
640      \overline\rho^{\,i+1/2}\,U\,\delta_{i+1/2}[z_t]   
641   +     \overline\rho^{\,j+1/2}\,V\,\delta_{j+1/2}[z_t]     
642   +  W\; \overline \rho^{\,k+1/2}\;\delta_{k+1/2} [z_t]   
643    -        \frac{p_t}{g}\,\partial_t b_t    \biggr\}    &&&\\
644%
645\allowdisplaybreaks
646%
647\equiv& - \sum\limits_{i,j,k}   g \; z_t      \biggl\{ 
648      \delta_i    \left[ U\;  \overline \rho^{\,i+1/2}   \right]
649   +  \delta_j    \left[ V\;  \overline \rho^{\,j+1/2}   \right]
650   +  \delta_k    \left[ W\;  \overline \rho^{\,k+1/2}   \right]       \biggr\}   
651             - \sum\limits_{i,j,k}       \biggl\{ p_t\;\partial_t b_t    \biggr\}   &&&\\ 
652%
653\equiv& + \sum\limits_{i,j,k}   g \; z_t    \biggl\{      \partial_t ( e_{3t} \,\rho)    \biggr\}  \; b_t   
654             -  \sum\limits_{i,j,k}                 \biggl\{  p_t\;\partial_t b_t                     \biggr\}              &&&\\   
655%
656\end{flalign*}
657The first term is exactly the first term of the right-hand-side of \eqref{KE+PE_vect_discrete}.
658It remains to demonstrate that the last term, which is obviously a discrete analogue of
659$\int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to the last term of \eqref{KE+PE_vect_discrete}.
660In other words, the following property must be satisfied:
661\begin{flalign*}
662           \sum\limits_{i,j,k}  \biggl\{  p_t\;\partial_t b_t                  \biggr\}       
663\equiv  \sum\limits_{i,j,k}  \biggl\{ \rho \,g\,\partial_t (z_t) \,b_\biggr\}                               
664\end{flalign*}
665
666Let introduce $p_w$ the pressure at $w$-point such that $\delta_k [p_w] = - \rho \,g\,e_{3t}$.
667The right-hand-side of the above equation can be transformed as follows:
668
669\begin{flalign*}
670               \sum\limits_{i,j,k}  \biggl\{ \rho \,g\,\partial_t (z_t) \,b_\biggr\} 
671&\equiv   - \sum\limits_{i,j,k}  \biggl\{ \delta_k [p_w]\,\partial_t (z_t) \,e_{1t}\,e_{2t}  \biggr\}        &&&\\
672%
673&\equiv  + \sum\limits_{i,j,k}  \biggl\{  p_w\, \delta_{k+1/2} [\partial_t (z_t)] \,e_{1t}\,e_{2t}  \biggr\} 
674  \equiv  + \sum\limits_{i,j,k}  \biggl\{  p_w\, \partial_t (e_{3w}) \,e_{1t}\,e_{2t}  \biggr\}        &&&\\
675&\equiv  + \sum\limits_{i,j,k}  \biggl\{  p_w\, \partial_t (b_w) \biggr\} 
676 %
677% & \equiv     \sum\limits_{i,j,k} \biggl\{  \frac{1}{e_{3t}} \delta_k [p_w]\;\partial_t (z_t) \,b_w   \right)   \biggr\}           &&&\\
678% & \equiv     \sum\limits_{i,j,k} \biggl\{   p_w\;\partial_t \left(    \delta_k [z_t]   \right)  e_{1w}\,e_{2w}   \biggr\}           &&&\\
679% & \equiv     \sum\limits_{i,j,k} \biggl\{   p_w\;\partial_t b_w   \biggr\}         
680\end{flalign*}
681therefore, the balance to be satisfied is:
682\begin{flalign*}
683           \sum\limits_{i,j,k}  \biggl\{  p_t\;\partial_t (b_t) \biggr\}  \equiv  \sum\limits_{i,j,k}  \biggl\{  p_w\, \partial_t (b_w) \biggr\} 
684\end{flalign*}
685which is a purely vertical balance:
686\begin{flalign*}
687           \sum\limits_{k}  \biggl\{  p_t\;\partial_t (e_{3t}) \biggr\}  \equiv  \sum\limits_{k}  \biggl\{  p_w\, \partial_t (e_{3w}) \biggr\} 
688\end{flalign*}
689Defining $p_w = \overline{p_t}^{\,k+1/2}$
690
691%gm comment
692\gmcomment{
693\begin{flalign*}
694 \sum\limits_{i,j,k} \biggl\{   p_t\;\partial_t b_t   \biggr\}                                &&&\\
695 %
696 & \equiv     \sum\limits_{i,j,k} \biggl\{  \frac{1}{e_{3t}} \delta_k [p_w]\;\partial_t (z_t) \,b_w    \biggr\}           &&&\\
697 & \equiv     \sum\limits_{i,j,k} \biggl\{   p_w\;\partial_t \left(    \delta_{k+1/2} [z_t]   \right)  e_{1w}\,e_{2w}   \biggr\}           &&&\\
698 & \equiv     \sum\limits_{i,j,k} \biggl\{   p_w\;\partial_t b_w   \biggr\}         
699\end{flalign*}
700
701
702\begin{flalign*}
703\int\limits_D   \rho \, g \, \frac{\partial z }{\partial t} \;dv
704\equiv&  \sum\limits_{i,j,k}   \biggl\{  \frac{1}{e_{3t}} \frac{\partial e_{3t}}{\partial t} p   \biggr\} \; b_t   &&&\\
705\equiv&  \sum\limits_{i,j,k}   \biggl\{      \biggr\} \; b_t   &&&\\
706\end{flalign*}
707
708%
709\begin{flalign*}
710\equiv& - \int_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv
711   + \int\limits_D g\, \rho \; \frac{\partial z}{\partial t}  \;dv     &&& \\
712\end{flalign*}
713%
714}
715%end gm comment
716
717
718Note that this property strongly constrains the discrete expression of both
719the depth of $T-$points and of the term added to the pressure gradient in the
720$s$-coordinate. Nevertheless, it is almost never satisfied since a linear equation
721of state is rarely used.
722
723
724
725
726
727
728
729% ================================================================
730% Discrete Total energy Conservation : flux form
731% ================================================================
732\section{Discrete total energy conservation : flux form}
733\label{Apdx_C.1}
734
735% -------------------------------------------------------------------------------------------------------------
736%       Total energy conservation
737% -------------------------------------------------------------------------------------------------------------
738\subsection{Total energy conservation}
739\label{Apdx_C_KE+PE}
740
741The discrete form of the total energy conservation, \eqref{Eq_Tot_Energy}, is given by:
742\begin{flalign*}
743\partial_t \left\sum\limits_{i,j,k} \biggl\{ \frac{u^2}{2} \,b_u + \frac{v^2}{2}\, b_v +  \rho \, g \, z_t \,b_\biggr\} \right) &=0  \\
744\end{flalign*}
745which in flux form, it leads to:
746\begin{flalign*}
747                        \sum\limits_{i,j,k} \biggl\{  \frac{u    }{e_{3u}}\,\frac{\partial (e_{3u}u)}{\partial t} \,b_u
748                                                             +  \frac{v    }{e_{3v}}\,\frac{\partial (e_{3v}v)}{\partial t} \,b_\biggr\}
749&  -  \frac{1}{2} \sum\limits_{i,j,k} \biggl\{  \frac{u^2}{e_{3u}}\frac{\partial    e_{3u}    }{\partial t} \,b_u
750                                                             +  \frac{v^2}{e_{3v}}\frac{\partial    e_{3v}    }{\partial t} \,b_v   \biggr\}      \\
751&= - \sum\limits_{i,j,k} \biggl\{ \frac{1}{e_3t}\frac{\partial e_{3t} \rho}{\partial t} \, g \, z_t \,b_\biggr\} 
752     - \sum\limits_{i,j,k} \biggl\{ \rho \,g\,\frac{\partial z_t}{\partial t} \,b_\biggr\}                                     \\
753\end{flalign*}
754
755Substituting the discrete expression of the time derivative of the velocity either in vector invariant or in flux form,
756leads to the discrete equivalent of the
757
758
759% -------------------------------------------------------------------------------------------------------------
760%       Coriolis and advection terms: flux form
761% -------------------------------------------------------------------------------------------------------------
762\subsection{Coriolis and advection terms: flux form}
763\label{Apdx_C.1.3}
764
765% -------------------------------------------------------------------------------------------------------------
766%       Coriolis plus ``metric'' Term
767% -------------------------------------------------------------------------------------------------------------
768\subsubsection{Coriolis plus ``metric'' Term}
769\label{Apdx_C.1.3.1} 
770
771In flux from the vorticity term reduces to a Coriolis term in which the Coriolis
772parameter has been modified to account for the ``metric'' term. This altered
773Coriolis parameter is discretised at an f-point. It is given by:
774\begin{equation*}
775f+\frac{1} {e_1 e_2 } \left( v \frac{\partial e_2 } {\partial i} - u \frac{\partial e_1 } {\partial j}\right)\;
776\equiv \;
777f+\frac{1} {e_{1f}\,e_{2f}} \left( \overline v^{\,i+1/2} \delta_{i+1/2} \left[ e_{2u} \right] 
778                                               -\overline u^{\,j+1/2} \delta_{j+1/2} \left[ e_{1u}  \right] \right)
779\end{equation*}
780
781Either the ENE or EEN scheme is then applied to obtain the vorticity term in flux form.
782It therefore conserves the total KE. The derivation is the same as for the
783vorticity term in the vector invariant form (\S\ref{Apdx_C_vor}).
784
785% -------------------------------------------------------------------------------------------------------------
786%       Flux form advection
787% -------------------------------------------------------------------------------------------------------------
788\subsubsection{Flux form advection}
789\label{Apdx_C.1.3.2} 
790
791The flux form operator of the momentum advection is evaluated using a
792centered second order finite difference scheme. Because of the flux form,
793the discrete operator does not contribute to the global budget of linear
794momentum. Because of the centered second order scheme, it conserves
795the horizontal kinetic energy, that is :
796
797\begin{equation} \label{Apdx_C_ADV_KE_flux}
798 -  \int_D \textbf{U}_h \cdot     \left(                 {{\begin{array} {*{20}c}
799\nabla \cdot \left( \textbf{U}\,u \right) \hfill \\
800\nabla \cdot \left( \textbf{U}\,v \right) \hfill \\       \end{array}} }           \right)   \;dv
801-   \frac{1}{2} \int_D {  {\textbf{U}_h}^2 \frac{1}{e_3} \frac{\partial  e_3 }{\partial t} \;dv } =\;0
802\end{equation}
803
804Let us first consider the first term of the scalar product ($i.e.$ just the the terms
805associated with the i-component of the advection) :
806\begin{flalign*}
807&  - \int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv   \\
808%
809\equiv& - \sum\limits_{i,j,k} \biggl\{    \frac{1} {b_u}    \biggl(   
810      \delta_{i+1/2}  \left[   \overline {U}^{\,i}      \;\overline u^{\,i}          \right]   
811   + \delta_j           \left[   \overline {V}^{\,i+1/2}\;\overline u^{\,j+1/2}   \right]   
812   + \delta_k          \left[   \overline {W}^{\,i+1/2}\;\overline u^{\,k+1/2} \right]  \biggr)   \;   \biggr\} \, b_u \;u &&&  \\ 
813%
814\equiv& - \sum\limits_{i,j,k} 
815   \biggl\{ 
816      \delta_{i+1/2} \left[   \overline {U}^{\,i}\;  \overline u^{\,i}  \right]
817   + \delta_j          \left[   \overline {V}^{\,i+1/2}\;\overline u^{\,j+1/2}   \right]
818   + \delta_k         \left[   \overline {W}^{\,i+12}\;\overline u^{\,k+1/2}  \right]
819      \; \biggr\} \; u     \\
820%
821\equiv& + \sum\limits_{i,j,k}
822   \biggl\{ 
823      \overline {U}^{\,i}\;   \overline u^{\,i}    \delta_i \left[ u \right] 
824        + \overline {V}^{\,i+1/2}\; \overline u^{\,j+1/2}   \delta_{j+1/2} \left[ u \right] 
825        + \overline {W}^{\,i+1/2}\; \overline u^{\,k+1/2}   \delta_{k+1/2}    \left[ u \right]     \biggr\}     && \\
826%
827\equiv& + \frac{1}{2} \sum\limits_{i,j,k}    \biggl\{ 
828       \overline{U}^{\,i}     \delta_i       \left[ u^2 \right] 
829    + \overline{V}^{\,i+1/2}  \delta_{j+/2}  \left[ u^2 \right]
830    + \overline{W}^{\,i+1/2}  \delta_{k+1/2}    \left[ u^2 \right]      \biggr\} && \\
831%
832\equiv& -  \sum\limits_{i,j,k}    \frac{1}{2}   \bigg\{ 
833       U  \; \delta_{i+1/2}    \left[ \overline {u^2}^{\,i} \right]
834         + V  \; \delta_{j+1/2}    \left[ \overline {u^2}^{\,i} \right]
835    + W \; \delta_{k+1/2}   \left[ \overline {u^2}^{\,i} \right]     \biggr\}    &&& \\
836%
837\equiv& - \sum\limits_{i,j,k}  \frac{1}{2}  \overline {u^2}^{\,i}     \biggl\{ 
838      \delta_{i+1/2}    \left[ U  \right]
839   + \delta_{j+1/2}  \left[ V  \right]
840   + \delta_{k+1/2}  \left[ W \right]     \biggr\}    &&& \\
841%
842\equiv& + \sum\limits_{i,j,k}  \frac{1}{2}  \overline {u^2}^{\,i} 
843   \biggl\{     \left(   \frac{1}{e_{3t}} \frac{\partial e_{3t}}{\partial t}   \right) \; b_t     \biggr\}    &&& \\
844\end{flalign*}
845Applying similar manipulation applied to the second term of the scalar product
846leads to :
847\begin{equation*}
848 -  \int_D \textbf{U}_h \cdot     \left(                 {{\begin{array} {*{20}c}
849\nabla \cdot \left( \textbf{U}\,u \right) \hfill \\
850\nabla \cdot \left( \textbf{U}\,v \right) \hfill \\       \end{array}} }           \right)   \;dv     
851\equiv + \sum\limits_{i,j,k}  \frac{1}{2}  \left( \overline {u^2}^{\,i} + \overline {v^2}^{\,j} \right)
852   \biggl\{     \left(   \frac{1}{e_{3t}} \frac{\partial e_{3t}}{\partial t}   \right) \; b_t     \biggr\}   
853\end{equation*}
854which is the discrete form of
855$ \frac{1}{2} \int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv $.
856\eqref{Apdx_C_ADV_KE_flux} is thus satisfied.
857
858
859When the UBS scheme is used to evaluate the flux form momentum advection,
860the discrete operator does not contribute to the global budget of linear momentum
861(flux form). The horizontal kinetic energy is not conserved, but forced to decay
862($i.e.$ the scheme is diffusive).
863
864
865
866
867
868
869
870
871
872
873% ================================================================
874% Discrete Enstrophy Conservation
875% ================================================================
876\section{Discrete enstrophy conservation}
877\label{Apdx_C.1}
878
879
880% -------------------------------------------------------------------------------------------------------------
881%       Vorticity Term with ENS scheme
882% -------------------------------------------------------------------------------------------------------------
883\subsubsection{Vorticity Term with ENS scheme  (\np{ln\_dynvor\_ens}=.true.)}
884\label{Apdx_C_vorENS} 
885
886In the ENS scheme, the vorticity term is descretized as follows:
887\begin{equation} \label{Eq_dynvor_ens}
888\left\{   \begin{aligned}
889+\frac{1}{e_{1u}} & \overline{q}^{\,i}  & {\overline{ \overline{\left( e_{1v}\,e_{3v}\;  v \right) } } }^{\,i, j+1/2}    \\
890- \frac{1}{e_{2v}} & \overline{q}^{\,j}  & {\overline{ \overline{\left( e_{2u}\,e_{3u}\; u \right) } } }^{\,i+1/2, j} 
891\end{aligned}  \right.
892\end{equation} 
893
894The scheme does not allow but the conservation of the total kinetic energy but the conservation
895of $q^2$, the potential enstrophy for a horizontally non-divergent flow ($i.e.$ when $\chi$=$0$).
896Indeed, using the symmetry or skew symmetry properties of the operators (Eqs \eqref{DOM_mi_adj} 
897and \eqref{DOM_di_adj}), it can be shown that:
898\begin{equation} \label{Apdx_C_1.1}
899\int_D {q\,\;{\textbf{k}}\cdot \frac{1} {e_3} \nabla \times \left( {e_3 \, q \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} \equiv 0
900\end{equation}
901where $dv=e_1\,e_2\,e_3 \; di\,dj\,dk$ is the volume element. Indeed, using
902\eqref{Eq_dynvor_ens}, the discrete form of the right hand side of \eqref{Apdx_C_1.1} 
903can be transformed as follow:
904\begin{flalign*} 
905&\int_D q \,\; \textbf{k} \cdot \frac{1} {e_3 } \nabla \times 
906   \left(  e_3 \, q \; \textbf{k} \times  \textbf{U}_h   \right)\;   dv       \\
907%
908& \qquad {\begin{array}{*{20}l} 
909&\equiv \sum\limits_{i,j,k} 
910q \ \left\{  \delta_{i+1/2}  \left[ - \,\overline {q}^{\,i}\;  \overline{\overline  U }^{\,i,j+1/ 2} \right] 
911             - \delta_{j+1/2} \left[    \overline {q}^{\,j}\;  \overline{\overline  V }^{\,i+1/2, j} \right]     \right\}    \\ 
912%
913&\equiv \sum\limits_{i,j,k} 
914   \left\{   \delta_i [q] \; \overline{q}^{\,i} \; \overline{  \overline U  }^{\,i,j+1/2} 
915      + \delta_j [q] \; \overline{q}^{\,j} \; \overline{\overline V }^{\,i+1/2,j}        \right\}       &&  \\ 
916%
917&\equiv \,\frac{1} {2} \sum\limits_{i,j,k} 
918   \left\{         \delta_\left[ q^2 \right] \; \overline{\overline U }^{\,i,j+1/2} 
919            + \delta_\left[ q^2 \right] \; \overline{\overline V }^{\,i+1/2,j}      \right\}    &&  \\ 
920%
921&\equiv - \frac{1} {2} \sum\limits_{i,j,k}   q^2 \;
922   \left\{    \delta_{i+1/2}   \left[   \overline{\overline{ U }}^{\,i,j+1/2}    \right] 
923            + \delta_{j+1/2}  \left[   \overline{\overline{ V }}^{\,i+1/2,j}     \right]    \right\}    && \\ 
924\end{array} }     
925%
926\allowdisplaybreaks
927\intertext{ Since $\overline {\;\cdot \;} $ and $\delta $ operators commute: $\delta_{i+1/2}
928\left[ {\overline a^{\,i}} \right] = \overline {\delta_i \left[ a \right]}^{\,i+1/2}$,
929and introducing the horizontal divergence $\chi $, it becomes: }
930\allowdisplaybreaks
931%
932& \qquad {\begin{array}{*{20}l} 
933&\equiv \sum\limits_{i,j,k} - \frac{1} {2} q^2 \; \overline{\overline{ e_{1t}\,e_{2t}\,e_{3t}^{}\, \chi}}^{\,i+1/2,j+1/2} 
934\quad \equiv 0 &&
935\end{array} }     
936\end{flalign*}
937The later equality is obtain only when the flow is horizontally non-divergent, $i.e.$ $\chi$=$0$.
938
939
940% -------------------------------------------------------------------------------------------------------------
941%       Vorticity Term with EEN scheme
942% -------------------------------------------------------------------------------------------------------------
943\subsubsection{Vorticity Term with EEN scheme (\np{ln\_dynvor\_een}=.true.)}
944\label{Apdx_C_vorEEN} 
945
946With the EEN scheme, the vorticity terms are represented as:
947\begin{equation} \label{Eq_dynvor_een}
948\left\{ {    \begin{aligned}
949 +q\,e_3 \, v  &\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}} 
950                         {^{i+1/2-i_p}_j}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{1v} e_{3v} \ \right)^{i+i_p-1/2}_{j+j_p}   \\
951 - q\,e_3 \, u     &\equiv -\frac{1}{e_{2v} }    \sum_{\substack{i_p,\,k_p}} 
952                         {^i_{j+1/2-j_p}}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{2u} e_{3u} \ \right)^{i+i_p}_{j+j_p-1/2}   \\
953\end{aligned}   } \right.
954\end{equation} 
955where the indices $i_p$ and $k_p$ take the following value:
956$i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$,
957and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by:
958\begin{equation} \label{Q_triads}
959_i^j \mathbb{Q}^{i_p}_{j_p}
960= \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right)
961\end{equation}
962
963
964This formulation does conserve the potential enstrophy for a horizontally non-divergent flow ($i.e.$ $\chi=0$).
965
966Let consider one of the vorticity triad, for example ${^{i}_j}\mathbb{Q}^{+1/2}_{+1/2} $,
967similar manipulation can be done for the 3 others. The discrete form of the right hand
968side of \eqref{Apdx_C_1.1} applied to this triad only can be transformed as follow:
969
970\begin{flalign*} 
971&\int_D {q\,\;{\textbf{k}}\cdot \frac{1} {e_3} \nabla \times \left( {e_3 \, q \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} \\
972%
973\equiv& \sum\limits_{i,j,k} 
974 {q} \    \biggl\{ \;\;
975   \delta_{i+1/2} \left[   -\, {{^i_j}\mathbb{Q}^{+1/2}_{+1/2}  \; U^{i+1/2}_{j}}    \right] 
976      - \delta_{j+1/2} \left[       {{^i_j}\mathbb{Q}^{+1/2}_{+1/2}  \; V^{i}_{j+1/2}}    \right] 
977   \;\;\biggr\}        &&  \\ 
978%
979\equiv& \sum\limits_{i,j,k} 
980   \biggl\{   \delta_i [q] \  {{^i_j}\mathbb{Q}^{+1/2}_{+1/2}  \; U^{i+1/2}_{j}}
981         + \delta_j [q] \  {{^i_j}\mathbb{Q}^{+1/2}_{+1/2}  \; V^{i}_{j+1/2}}   \biggr\}
982      && \\ 
983%
984... & &&\\
985&Demonstation \ to \ be \ done... &&\\
986... & &&\\
987%
988\equiv& \frac{1} {2} \sum\limits_{i,j,k} 
989   \biggl\{ \delta_i    \Bigl[    \left(  {{^i_j}\mathbb{Q}^{+1/2}_{+1/2}} \right)^2   \Bigr]\;
990         \overline{\overline {U}}^{\,i,j+1/2} 
991            + \delta_\Bigl[    \left( {{^i_j}\mathbb{Q}^{+1/2}_{+1/2}} \right)^2     \Bigr]\; 
992         \overline{\overline {V}}^{\,i+1/2,j} 
993   \biggr\} 
994   &&  \\ 
995%
996\equiv& - \frac{1} {2} \sum\limits_{i,j,k}   \left( {{^i_j}\mathbb{Q}^{+1/2}_{+1/2}} \right)^2\;
997   \biggl\{    \delta_{i+1/2} 
998         \left[   \overline{\overline {U}}^{\,i,j+1/2}    \right] 
999               + \delta_{j+1/2}
1000      \left[   \overline{\overline {V}}^{\,i+1/2,j}     \right] 
1001   \biggr\}    && \\ 
1002%
1003\equiv& \sum\limits_{i,j,k} - \frac{1} {2} \left( {{^i_j}\mathbb{Q}^{+1/2}_{+1/2}} \right)^2
1004                                                            \; \overline{\overline{ b_t^{}\, \chi}}^{\,i+1/2,\,j+1/2}  &&\\
1005%
1006\ \ \equiv& \ 0     &&\\
1007\end{flalign*}
1008
1009
1010
1011
1012
1013% ================================================================
1014% Conservation Properties on Tracers
1015% ================================================================
1016\section{Conservation Properties on Tracers}
1017\label{Apdx_C.2}
1018
1019
1020All the numerical schemes used in NEMO are written such that the tracer content
1021is conserved by the internal dynamics and physics (equations in flux form).
1022For advection, only the CEN2 scheme ($i.e.$ $2^{nd}$ order finite different scheme)
1023conserves the global variance of tracer. Nevertheless the other schemes ensure
1024that the global variance decreases ($i.e.$ they are at least slightly diffusive).
1025For diffusion, all the schemes ensure the decrease of the total tracer variance,
1026except the iso-neutral operator. There is generally no strict conservation of mass,
1027as the equation of state is non linear with respect to $T$ and $S$. In practice,
1028the mass is conserved to a very high accuracy.
1029% -------------------------------------------------------------------------------------------------------------
1030%       Advection Term
1031% -------------------------------------------------------------------------------------------------------------
1032\subsection{Advection Term}
1033\label{Apdx_C.2.1}
1034
1035conservation of a tracer, $T$:
1036\begin{equation*}
1037\frac{\partial }{\partial t} \left(   \int_D {T\;dv}   \right)
1038\int_D { \frac{1}{e_3}\frac{\partial \left( e_3 \, T \right)}{\partial t} \;dv }=0
1039\end{equation*}
1040
1041conservation of its variance:
1042\begin{flalign*} 
1043\frac{\partial }{\partial t} \left( \int_D {\frac{1}{2} T^2\;dv} \right)
1044=&  \int_D { \frac{1}{e_3} Q      \frac{\partial \left( e_3 \, T \right) }{\partial t} \;dv } 
1045-   \frac{1}{2} \int_D {  T^2 \frac{1}{e_3} \frac{\partial  e_3 }{\partial t} \;dv }
1046\end{flalign*}
1047
1048
1049Whatever the advection scheme considered it conserves of the tracer content as all
1050the scheme are written in flux form. Indeed,  let $T$ be the tracer and $\tau_u$, $\tau_v$,
1051and $\tau_w$ its interpolated values at velocity point (whatever the interpolation is),
1052the conservation of the tracer content due to the advection tendency is obtained as follows:
1053\begin{flalign*}
1054&\int_D { \frac{1}{e_3}\frac{\partial \left( e_3 \, T \right)}{\partial t} \;dv } = - \int_D \nabla \cdot \left( T \textbf{U} \right)\;dv    &&&\\
1055&\equiv - \sum\limits_{i,j,k}    \biggl\{
1056    \frac{1} {b_t}  \left\delta_i    \left[   U \;\tau_u   \right]
1057                                + \delta_j    \left[   V  \;\tau_v   \right] \right)
1058   + \frac{1} {e_{3t}} \delta_k \left[ w\;\tau_w \right]    \biggl\}  b_t   &&&\\
1059%
1060&\equiv - \sum\limits_{i,j,k}     \left\{
1061       \delta_\left[   U \;\tau_u   \right]
1062         + \delta_\left[   V  \;\tau_v   \right]
1063    + \delta_k \left[ W \;\tau_w \right] \right\}   && \\
1064&\equiv 0 &&&
1065\end{flalign*}
1066
1067The conservation of the variance of tracer due to the advection tendency
1068can be achieved only with the CEN2 scheme, $i.e.$ when
1069$\tau_u= \overline T^{\,i+1/2}$, $\tau_v= \overline T^{\,j+1/2}$, and $\tau_w= \overline T^{\,k+1/2}$.
1070It can be demonstarted as follows:
1071\begin{flalign*}
1072&\int_D { \frac{1}{e_3} Q      \frac{\partial \left( e_3 \, T \right) }{\partial t} \;dv }
1073= - \int\limits_D \tau\;\nabla \cdot \left( T\; \textbf{U} \right)\;dv &&&\\
1074\equiv& - \sum\limits_{i,j,k} T\;
1075   \left\{
1076      \delta_\left[ U  \overline T^{\,i+1/2}  \right]
1077   + \delta_\left[ V  \overline T^{\,j+1/2}  \right]
1078   + \delta_k \left[ W \overline T^{\,k+1/2} \right]          \right\} && \\
1079\equiv& + \sum\limits_{i,j,k} 
1080   \left\{     U  \overline T^{\,i+1/2} \,\delta_{i+1/2}  \left[ T \right] 
1081                 +  V  \overline T^{\,j+1/2} \;\delta_{j+1/2}  \left[ T \right]
1082                 +  W \overline T^{\,k+1/2}\;\delta_{k+1/2} \left[ T \right]     \right\}      &&\\
1083\equiv&  + \frac{1} {2}  \sum\limits_{i,j,k}
1084   \Bigl\{   U  \;\delta_{i+1/2} \left[ T^2 \right]
1085                 + V  \;\delta_{j+1/2}  \left[ T^2 \right]
1086                 + W \;\delta_{k+1/2} \left[ T^2 \right]   \Bigr\}     && \\ 
1087\equiv& - \frac{1} {2}  \sum\limits_{i,j,k} T^2
1088   \Bigl\{    \delta_\left[ U  \right]
1089                  + \delta_\left[ V  \right]
1090                  + \delta_k \left[ W \right]     \Bigr\}      &&&  \\
1091\equiv& + \frac{1} {2}  \sum\limits_{i,j,k} T^2
1092   \Bigl\{   \frac{1}{e_{3t}} \frac{\partial e_{3t}\,T }{\partial t}     \Bigr\}      &&& \\
1093\end{flalign*}
1094which is the discrete form of $ \frac{1}{2} \int_D {  T^2 \frac{1}{e_3} \frac{\partial  e_3 }{\partial t} \;dv }$.
1095
1096% ================================================================
1097% Conservation Properties on Lateral Momentum Physics
1098% ================================================================
1099\section{Conservation Properties on Lateral Momentum Physics}
1100\label{Apdx_dynldf_properties}
1101
1102
1103The discrete formulation of the horizontal diffusion of momentum ensures the
1104conservation of potential vorticity and the horizontal divergence, and the
1105dissipation of the square of these quantities (i.e. enstrophy and the
1106variance of the horizontal divergence) as well as the dissipation of the
1107horizontal kinetic energy. In particular, when the eddy coefficients are
1108horizontally uniform, it ensures a complete separation of vorticity and
1109horizontal divergence fields, so that diffusion (dissipation) of vorticity
1110(enstrophy) does not generate horizontal divergence (variance of the
1111horizontal divergence) and \textit{vice versa}.
1112
1113These properties of the horizontal diffusion operator are a direct consequence
1114of properties \eqref{Eq_DOM_curl_grad} and \eqref{Eq_DOM_div_curl}.
1115When the vertical curl of the horizontal diffusion of momentum (discrete sense)
1116is taken, the term associated with the horizontal gradient of the divergence is
1117locally zero.
1118
1119% -------------------------------------------------------------------------------------------------------------
1120%       Conservation of Potential Vorticity
1121% -------------------------------------------------------------------------------------------------------------
1122\subsection{Conservation of Potential Vorticity}
1123\label{Apdx_C.3.1}
1124
1125The lateral momentum diffusion term conserves the potential vorticity :
1126\begin{flalign*}
1127&\int \limits_D \frac{1} {e_3 } \textbf{k} \cdot \nabla \times 
1128   \Bigl[    \nabla_\left( A^{\,lm}\;\chi  \right)
1129             - \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right)    \Bigr]\;dv  = 0
1130\end{flalign*}
1131%%%%%%%%%%  recheck here....  (gm)
1132\begin{flalign*}
1133= \int \limits_D  -\frac{1} {e_3 } \textbf{k} \cdot \nabla \times 
1134   \Bigl[ \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right)  \Bigr]\;dv &&& \\ 
1135\end{flalign*}
1136\begin{flalign*}
1137\equiv& \sum\limits_{i,j}
1138   \left\{
1139   \delta_{i+1/2} 
1140   \left[
1141   \frac {e_{2v}} {e_{1v}\,e_{3v}}  \delta_i
1142      \left[ A_f^{\,lm} e_{3f} \zeta  \right]
1143    \right]
1144   + \delta_{j+1/2} 
1145   \left[
1146   \frac {e_{1u}} {e_{2u}\,e_{3u}} \delta_j
1147      \left[ A_f^{\,lm} e_{3f} \zeta  \right]
1148   \right]
1149   \right\} 
1150   && \\ 
1151%
1152\intertext{Using \eqref{DOM_di_adj}, it follows:}
1153%
1154\equiv& \sum\limits_{i,j,k} 
1155   -\,\left\{
1156      \frac{e_{2v}} {e_{1v}\,e_{3v}}  \delta_i
1157      \left[ A_f^{\,lm} e_{3f} \zeta  \right]\;\delta_i \left[ 1\right]
1158   + \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j
1159      \left[ A_f^{\,lm} e_{3f} \zeta  \right]\;\delta_j \left[ 1\right]
1160   \right\} \quad \equiv 0
1161   && \\ 
1162\end{flalign*}
1163
1164% -------------------------------------------------------------------------------------------------------------
1165%       Dissipation of Horizontal Kinetic Energy
1166% -------------------------------------------------------------------------------------------------------------
1167\subsection{Dissipation of Horizontal Kinetic Energy}
1168\label{Apdx_C.3.2}
1169
1170
1171The lateral momentum diffusion term dissipates the horizontal kinetic energy:
1172%\begin{flalign*}
1173\begin{equation*}
1174\begin{split}
1175\int_D \textbf{U}_h \cdot 
1176   \left[ \nabla_h      \right.   &     \left.       \left( A^{\,lm}\;\chi \right)     
1177   - \nabla_h \times  \left( A^{\,lm}\;\zeta \;\textbf{k} \right)     \right] \; dv    \\
1178\\  %%%
1179\equiv& \sum\limits_{i,j,k} 
1180   \left\{
1181     \frac{1} {e_{1u}}               \delta_{i+1/2} \left[  A_T^{\,lm}          \chi     \right]
1182   - \frac{1} {e_{2u}\,e_{3u}}  \delta_j           \left[ A_f^{\,lm} e_{3f} \zeta   \right]
1183   \right\} \; e_{1u}\,e_{2u}\,e_{3u} \;u     \\
1184&\;\; +  \left\{
1185      \frac{1} {e_{2u}}             \delta_{j+1/2} \left[ A_T^{\,lm}          \chi    \right] 
1186   + \frac{1} {e_{1v}\,e_{3v}} \delta_i            \left[ A_f^{\,lm} e_{3f} \zeta  \right] 
1187   \right\} \; e_{1v}\,e_{2u}\,e_{3v} \;v     \qquad \\ 
1188\\  %%%
1189\equiv& \sum\limits_{i,j,k} 
1190   \Bigl\{
1191     e_{2u}\,e_{3u} \;u\;  \delta_{i+1/2} \left[ A_T^{\,lm}           \chi    \right]
1192   - e_{1u}             \;u\;  \delta_j           \left[ A_f^{\,lm} e_{3f} \zeta  \right]
1193    \Bigl\} 
1194    \\ 
1195&\;\; + \Bigl\{
1196      e_{1v}\,e_{3v} \;v\;  \delta_{j+1/2}  \left[ A_T^{\,lm}           \chi    \right]
1197   + e_{2v}             \;v\;  \delta_i            \left[ A_f^{\,lm} e_{3f} \zeta  \right]
1198   \Bigl\}      \\ 
1199\\  %%%
1200\equiv& \sum\limits_{i,j,k} 
1201   - \Bigl(
1202     \delta_i   \left[  e_{2u}\,e_{3u} \;u  \right]
1203   + \delta_\left[  e_{1v}\,e_{3v}  \;v  \right] 
1204        \Bigr) \;  A_T^{\,lm} \chi   \\ 
1205&\;\; - \Bigl(
1206     \delta_{i+1/2}  \left[  e_{2v}  \;v  \right]
1207   - \delta_{j+1/2}  \left[  e_{1u} \;u  \right] 
1208        \Bigr)\;  A_f^{\,lm} e_{3f} \zeta      \\ 
1209\\  %%%
1210\equiv& \sum\limits_{i,j,k} 
1211   - A_T^{\,lm}  \,\chi^2   \;e_{1t}\,e_{2t}\,e_{3t}
1212   - A_f ^{\,lm}  \,\zeta^2 \;e_{1f }\,e_{2f }\,e_{3f} 
1213\quad \leq 0       \\
1214\end{split}
1215\end{equation*}
1216
1217% -------------------------------------------------------------------------------------------------------------
1218%       Dissipation of Enstrophy
1219% -------------------------------------------------------------------------------------------------------------
1220\subsection{Dissipation of Enstrophy}
1221\label{Apdx_C.3.3}
1222
1223
1224The lateral momentum diffusion term dissipates the enstrophy when the eddy
1225coefficients are horizontally uniform:
1226\begin{flalign*}
1227&\int\limits_\zeta \; \textbf{k} \cdot \nabla \times 
1228   \left[   \nabla_h \left( A^{\,lm}\;\chi  \right)
1229          - \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right)   \right]\;dv &&&\\
1230&= A^{\,lm} \int \limits_D \zeta \textbf{k} \cdot \nabla \times 
1231   \left[    \nabla_h \times \left( \zeta \; \textbf{k} \right)   \right]\;dv &&&\\
1232&\equiv A^{\,lm} \sum\limits_{i,j,k}  \zeta \;e_{3f} 
1233   \left\{     \delta_{i+1/2} \left[  \frac{e_{2v}} {e_{1v}\,e_{3v}} \delta_i \left[ e_{3f} \zeta  \right]   \right]
1234             + \delta_{j+1/2} \left[  \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j \left[ e_{3f} \zeta  \right]   \right]      \right\}   &&&\\ 
1235%
1236\intertext{Using \eqref{DOM_di_adj}, it follows:}
1237%
1238&\equiv  - A^{\,lm} \sum\limits_{i,j,k} 
1239   \left\{    \left\frac{1} {e_{1v}\,e_{3v}}  \delta_i \left[ e_{3f} \zeta  \right]  \right)^2   b_v
1240            + \left\frac{1} {e_{2u}\,e_{3u}}  \delta_j \left[ e_{3f} \zeta  \right] \right)^2   b_\right\}      &&&\\
1241& \leq \;0       &&&\\ 
1242\end{flalign*}
1243
1244% -------------------------------------------------------------------------------------------------------------
1245%       Conservation of Horizontal Divergence
1246% -------------------------------------------------------------------------------------------------------------
1247\subsection{Conservation of Horizontal Divergence}
1248\label{Apdx_C.3.4}
1249
1250When the horizontal divergence of the horizontal diffusion of momentum
1251(discrete sense) is taken, the term associated with the vertical curl of the
1252vorticity is zero locally, due to (!!! II.1.8  !!!!!). The resulting term conserves the
1253$\chi$ and dissipates $\chi^2$ when the eddy coefficients are
1254horizontally uniform.
1255\begin{flalign*}
1256& \int\limits_\nabla_h \cdot 
1257   \Bigl[     \nabla_h \left( A^{\,lm}\;\chi \right)
1258             - \nabla_h \times \left( A^{\,lm}\;\zeta \;\textbf{k} \right)    \Bigr]  dv
1259= \int\limits_\nabla_h \cdot \nabla_h \left( A^{\,lm}\;\chi  \right)   dv   &&&\\
1260%
1261&\equiv \sum\limits_{i,j,k} 
1262   \left\{   \delta_i \left[ A_u^{\,lm} \frac{e_{2u}\,e_{3u}} {e_{1u}}  \delta_{i+1/2} \left[ \chi \right]  \right]
1263           + \delta_j \left[ A_v^{\,lm} \frac{e_{1v}\,e_{3v}} {e_{2v}}  \delta_{j+1/2} \left[ \chi \right]  \right]    \right\}    &&&\\ 
1264%
1265\intertext{Using \eqref{DOM_di_adj}, it follows:}
1266%
1267&\equiv \sum\limits_{i,j,k} 
1268   - \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] 
1269             + \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\} 
1270   \qquad \equiv 0     &&& \\ 
1271\end{flalign*}
1272
1273% -------------------------------------------------------------------------------------------------------------
1274%       Dissipation of Horizontal Divergence Variance
1275% -------------------------------------------------------------------------------------------------------------
1276\subsection{Dissipation of Horizontal Divergence Variance}
1277\label{Apdx_C.3.5}
1278
1279\begin{flalign*}
1280&\int\limits_D \chi \;\nabla_h \cdot 
1281   \left[    \nabla_h              \left( A^{\,lm}\;\chi                    \right)
1282           - \nabla_h   \times  \left( A^{\,lm}\;\zeta \;\textbf{k} \right)    \right]\;  dv
1283 = A^{\,lm}   \int\limits_D \chi \;\nabla_h \cdot \nabla_h \left( \chi \right)\;  dv    &&&\\ 
1284%
1285&\equiv A^{\,lm}  \sum\limits_{i,j,k}  \frac{1} {e_{1t}\,e_{2t}\,e_{3t}}  \chi 
1286   \left\{
1287      \delta_\left[   \frac{e_{2u}\,e_{3u}} {e_{1u}}  \delta_{i+1/2} \left[ \chi \right]   \right]
1288   + \delta_\left[   \frac{e_{1v}\,e_{3v}} {e_{2v}}   \delta_{j+1/2} \left[ \chi \right]   \right]
1289   \right\} \;   e_{1t}\,e_{2t}\,e_{3t}    &&&\\ 
1290%
1291\intertext{Using \eqref{DOM_di_adj}, it turns out to be:}
1292%
1293&\equiv - A^{\,lm} \sum\limits_{i,j,k}
1294   \left\{    \left\frac{1} {e_{1u}}  \delta_{i+1/2}  \left[ \chi \right]  \right)^2  b_u
1295                 + \left\frac{1} {e_{2v}}  \delta_{j+1/2}  \left[ \chi \right]  \right)^2  b_v    \right\} \;    &&&\\
1296%
1297&\leq 0              &&&\\
1298\end{flalign*}
1299
1300% ================================================================
1301% Conservation Properties on Vertical Momentum Physics
1302% ================================================================
1303\section{Conservation Properties on Vertical Momentum Physics}
1304\label{Apdx_C_4}
1305
1306
1307As for the lateral momentum physics, the continuous form of the vertical diffusion
1308of momentum satisfies several integral constraints. The first two are associated
1309with the conservation of momentum and the dissipation of horizontal kinetic energy:
1310\begin{align*}
1311 \int\limits_D   \frac{1} {e_3 }\; \frac{\partial } {\partial k} 
1312       \left(   \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}   \right)\;  dv
1313   \qquad \quad &= \vec{\textbf{0}}    \\
1314%
1315\intertext{and}
1316%
1317\int\limits_D
1318   \textbf{U}_h \cdot   \frac{1} {e_3 }\; \frac{\partial } {\partial k}
1319   \left(   \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}   \right)\; dv    \quad &\leq 0     \\
1320\end{align*}
1321The first property is obvious. The second results from:
1322
1323\begin{flalign*}
1324\int\limits_D
1325   \textbf{U}_h \cdot   \frac{1} {e_3 }\; \frac{\partial } {\partial k}
1326   \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}   \right)\;dv    &&&\\
1327\end{flalign*}
1328\begin{flalign*}
1329&\equiv \sum\limits_{i,j,k} 
1330   \left(
1331      u\; \delta_k \left[ \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2}  \left[ u \right]  \right]\;  e_{1u}\,e_{2u} 
1332   + v\; \delta_k \left[ \frac{A_v^{\,vm}} {e_{3vw}} \delta_{k+1/2}   \left[ v \right]  \right]\;  e_{1v}\,e_{2v} \right)   &&&\\ 
1333%
1334\intertext{since the horizontal scale factor does not depend on $k$, it follows:}
1335%
1336&\equiv - \sum\limits_{i,j,k} 
1337   \left\frac{A_u^{\,vm}} {e_{3uw}} \left( \delta_{k+1/2} \left[ u \right] \right)^2\; e_{1u}\,e_{2u} 
1338         + \frac{A_v^{\,vm}} {e_{3vw}}  \left( \delta_{k+1/2} \left[ v \right] \right)^2\; e_{1v}\,e_{2v}  \right)
1339\quad \leq 0   &&&\\
1340\end{flalign*}
1341
1342The vorticity is also conserved. Indeed:
1343\begin{flalign*}
1344\int \limits_D
1345   \frac{1} {e_3 } \textbf{k} \cdot \nabla \times 
1346      \left( \frac{1} {e_3 }\; \frac{\partial } {\partial k}  \left(
1347              \frac{A^{\,vm}} {e_3}\; \frac{\partial \textbf{U}_h } {\partial k} 
1348        \right\right)\; dv   &&&\\ 
1349\end{flalign*}
1350\begin{flalign*}
1351\equiv \sum\limits_{i,j,k}  \frac{1} {e_{3f}}\; \frac{1} {e_{1f}\,e_{2f}}
1352   \bigg\{    \biggr.   \quad
1353   \delta_{i+1/2} 
1354      &\left(   \frac{e_{2v}} {e_{3v}} \delta_\left[  \frac{1} {e_{3vw}} \delta_{k+1/2} \left[ v \right]  \right\right)   &&\\
1355   \biggl.
1356   - \delta_{j+1/2} 
1357      &\left(   \frac{e_{1u}} {e_{3u}} \delta_k \left[  \frac{1} {e_{3uw}}\delta_{k+1/2} \left[ u \right]  \right]   \right)
1358   \biggr\} \;
1359   e_{1f}\,e_{2f}\,e_{3f} \; \equiv 0   && \\
1360\end{flalign*}
1361If the vertical diffusion coefficient is uniform over the whole domain, the
1362enstrophy is dissipated, $i.e.$
1363\begin{flalign*}
1364\int\limits_D \zeta \, \textbf{k} \cdot \nabla \times 
1365   \left(   \frac{1} {e_3}\; \frac{\partial } {\partial k}
1366      \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}   \right)   \right)\; dv = 0   &&&\\
1367\end{flalign*}
1368This property is only satisfied in $z$-coordinates:
1369
1370\begin{flalign*}
1371\int\limits_D \zeta \, \textbf{k} \cdot \nabla \times 
1372   \left\frac{1} {e_3}\; \frac{\partial } {\partial k}
1373      \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}  \right)   \right)\; dv   &&& \\ 
1374\end{flalign*}
1375\begin{flalign*}
1376\equiv \sum\limits_{i,j,k} \zeta \;e_{3f} \;
1377   \biggl\{    \biggr\quad
1378   \delta_{i+1/2} 
1379      &\left(   \frac{e_{2v}} {e_{3v}} \delta_k \left[ \frac{A_v^{\,vm}} {e_{3vw}} \delta_{k+1/2}[v]  \right]   \right)   &&\\
1380   - \delta_{j+1/2} 
1381      &\biggl.
1382      \left(   \frac{e_{1u}} {e_{3u}} \delta_k \left[ \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2} [u]  \right]    \right) \biggr\}   &&\\ 
1383\end{flalign*}
1384\begin{flalign*}
1385\equiv \sum\limits_{i,j,k} \zeta \;e_{3f} 
1386   \biggl\{    \biggr\quad
1387   \frac{1} {e_{3v}} \delta_k
1388      &\left[ \frac{A_v^{\,vm}} {e_{3vw}} \delta_{k+1/2} \left[ \delta_{i+1/2} \left[ e_{2v}\,v \right] \right]   \right]    &&\\ 
1389    \biggl.
1390   - \frac{1} {e_{3u}} \delta_k
1391      &\left[  \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2} \left[ \delta_{j+1/2} \left[ e_{1u}\,u \right] \right\right\biggr\}  &&\\ 
1392\end{flalign*}
1393Using the fact that the vertical diffusion coefficients are uniform, and that in
1394$z$-coordinate, the vertical scale factors do not depend on $i$ and $j$ so
1395that: $e_{3f} =e_{3u} =e_{3v} =e_{3t} $ and $e_{3w} =e_{3uw} =e_{3vw} $,
1396it follows:
1397\begin{flalign*}
1398\equiv A^{\,vm} \sum\limits_{i,j,k} \zeta \;\delta_k
1399   \left[   \frac{1} {e_{3w}} \delta_{k+1/2} \Bigl[   \delta_{i+1/2} \left[ e_{2v}\,v \right]
1400                                - \delta_{j+1/ 2} \left[ e_{1u}\,u \right]   \Bigr]    \right]    &&&\\
1401\end{flalign*}
1402\begin{flalign*}
1403\equiv - A^{\,vm} \sum\limits_{i,j,k} \frac{1} {e_{3w}}
1404   \left( \delta_{k+1/2} \left[ \zeta  \right] \right)^2 \; e_{1f}\,e_{2f}  \; \leq 0    &&&\\
1405\end{flalign*}
1406Similarly, the horizontal divergence is obviously conserved:
1407
1408\begin{flalign*}
1409\int\limits_D \nabla \cdot 
1410\left( \frac{1} {e_3 }\; \frac{\partial } {\partial k}
1411      \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} \right) \right)\; dv = 0    &&&\\
1412\end{flalign*}
1413and the square of the horizontal divergence decreases ($i.e.$ the horizontal
1414divergence is dissipated) if the vertical diffusion coefficient is uniform over the
1415whole domain:
1416
1417\begin{flalign*}
1418\int\limits_D \chi \;\nabla \cdot 
1419\left( \frac{1} {e_3 }\; \frac{\partial } {\partial k}
1420      \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k}  \right) \right)\;  dv = 0  &&&\\
1421\end{flalign*}
1422This property is only satisfied in the $z$-coordinate:
1423\begin{flalign*}
1424\int\limits_D \chi \;\nabla \cdot 
1425\left( \frac{1} {e_3 }\; \frac{\partial } {\partial k}
1426      \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} \right\right)\; dv    &&&\\
1427\end{flalign*}
1428\begin{flalign*}
1429\equiv \sum\limits_{i,j,k} \frac{\chi } {e_{1t}\,e_{2t}}
1430   \biggl\{    \Biggr\quad
1431   \delta_{i+1/2} 
1432      &\left(   \frac{e_{2u}} {e_{3u}} \delta_k
1433            \left[ \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2} [u] \right] \right)    &&\\ 
1434   \Biggl.
1435   + \delta_{j+1/2} 
1436      &\left( \frac{e_{1v}} {e_{3v}} \delta_k
1437         \left[ \frac{A_v^{\,vm}} {e_{3vw}} \delta_{k+1/2} [v] \right]   \right)
1438   \Biggr\} \;  e_{1t}\,e_{2t}\,e_{3t}   &&\\ 
1439\end{flalign*}
1440
1441\begin{flalign*}
1442\equiv A^{\,vm} \sum\limits_{i,j,k}  \chi \,
1443   \biggl\{ \biggr\quad
1444   \delta_{i+1/2}
1445      &\left(
1446         \delta_k \left[
1447         \frac{1} {e_{3uw}} \delta_{k+1/2} \left[ e_{2u}\,u \right] \right]   \right)    && \\ 
1448   \biggl.
1449   + \delta_{j+1/2} 
1450      &\left(    \delta_k \left[
1451         \frac{1} {e_{3vw}} \delta_{k+1/2} \left[ e_{1v}\,v \right] \right]   \right)   \biggr\}    && \\ 
1452\end{flalign*}
1453
1454\begin{flalign*}
1455\equiv -A^{\,vm} \sum\limits_{i,j,k} 
1456\frac{\delta_{k+1/2} \left[ \chi \right]} {e_{3w}}\; \biggl\{ 
1457   \delta_{k+1/2} \Bigl[
1458         \delta_{i+1/2} \left[ e_{2u}\,u \right]
1459      + \delta_{j+1/2} \left[ e_{1v}\,v \right]  \Bigr]    \biggr\}    &&&\\
1460\end{flalign*}
1461
1462\begin{flalign*}
1463\equiv -A^{\,vm} \sum\limits_{i,j,k}
1464 \frac{1} {e_{3w}} \delta_{k+1/2} \left[ \chi \right]\; \delta_{k+1/2} \left[ e_{1t}\,e_{2t} \;\chi \right]   &&&\\
1465\end{flalign*}
1466
1467\begin{flalign*}
1468\equiv -A^{\,vm} \sum\limits_{i,j,k} 
1469\frac{e_{1t}\,e_{2t}} {e_{3w}}\; \left( \delta_{k+1/2} \left[ \chi \right]  \right)^2     \quad  \equiv 0    &&&\\
1470\end{flalign*}
1471
1472% ================================================================
1473% Conservation Properties on Tracer Physics
1474% ================================================================
1475\section{Conservation Properties on Tracer Physics}
1476\label{Apdx_C.5}
1477
1478The numerical schemes used for tracer subgridscale physics are written such
1479that the heat and salt contents are conserved (equations in flux form, second
1480order centered finite differences). Since a flux form is used to compute the
1481temperature and salinity, the quadratic form of these quantities (i.e. their variance)
1482globally tends to diminish. As for the advection term, there is generally no strict
1483conservation of mass, even if in practice the mass is conserved to a very high
1484accuracy.
1485
1486% -------------------------------------------------------------------------------------------------------------
1487%       Conservation of Tracers
1488% -------------------------------------------------------------------------------------------------------------
1489\subsection{Conservation of Tracers}
1490\label{Apdx_C.5.1}
1491
1492constraint of conservation of tracers:
1493\begin{flalign*}
1494&\int\limits_\nabla  \cdot \left( A\;\nabla T \right)\;dv  &&&\\ 
1495\\
1496&\equiv \sum\limits_{i,j,k} 
1497   \biggl\{    \biggr.
1498   \delta_i
1499      \left[
1500      A_u^{\,lT} \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} 
1501         \left[ T \right]
1502      \right]
1503   + \delta_j
1504      \left[
1505      A_v^{\,lT} \frac{e_{1v}\,e_{3v}} {e_{2v}} \delta_{j+1/2} 
1506         \left[ T \right] 
1507      \right]
1508   &&\\  & \qquad \qquad \qquad \qquad \qquad \qquad \quad \;\;\;
1509   + \delta_k
1510      \left[
1511      A_w^{\,vT} \frac{e_{1t}\,e_{2t}} {e_{3t}} \delta_{k+1/2} 
1512         \left[ T \right] 
1513      \right]
1514   \biggr\}   \quad  \equiv 0
1515   &&\\ 
1516\end{flalign*}
1517
1518In fact, this property simply results from the flux form of the operator.
1519
1520% -------------------------------------------------------------------------------------------------------------
1521%       Dissipation of Tracer Variance
1522% -------------------------------------------------------------------------------------------------------------
1523\subsection{Dissipation of Tracer Variance}
1524\label{Apdx_C.5.2}
1525
1526constraint on the dissipation of tracer variance:
1527\begin{flalign*}
1528\int\limits_D T\;\nabla & \cdot \left( A\;\nabla T \right)\;dv &&&\\ 
1529&\equiv   \sum\limits_{i,j,k} \; T
1530\biggl\{  \biggr.
1531     \delta_i \left[ A_u^{\,lT} \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} \left[T\right] \right]
1532& + \delta_j \left[ A_v^{\,lT} \frac{e_{1v} \,e_{3v}} {e_{2v}} \delta_{j+1/2} \left[T\right] \right]
1533      \quad&& \\ 
1534 \biggl.
1535&&+ \delta_k \left[A_w^{\,vT}\frac{e_{1t}\,e_{2t}} {e_{3t}}\delta_{k+1/2}\left[T\right]\right]
1536\biggr\} && 
1537\end{flalign*}
1538\begin{flalign*}
1539\equiv - \sum\limits_{i,j,k} 
1540   \biggl\{    \biggr\quad
1541   &    A_u^{\,lT} \left( \frac{1} {e_{1u}} \delta_{i+1/2} \left[ T \right]  \right)^2   e_{1u}\,e_{2u}\,e_{3u}    && \\
1542   & + A_v^{\,lT} \left( \frac{1} {e_{2v}} \delta_{j+1/2} \left[ T \right]  \right)^2   e_{1v}\,e_{2v}\,e_{3v}     && \\ \biggl.
1543   & + A_w^{\,vT} \left( \frac{1} {e_{3w}} \delta_{k+1/2} \left[ T \right]   \right)^2    e_{1w}\,e_{2w}\,e_{3w}   \biggr\} 
1544   \quad      \leq 0      && \\ 
1545\end{flalign*}
1546
1547
1548%%%%  end of appendix in gm comment
1549%}
Note: See TracBrowser for help on using the repository browser.