# source:NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/annex_C.tex@10368

Last change on this file since 10368 was 10368, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10365, see #2133

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