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

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

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133

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