# source:NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex@11693

Last change on this file since 11693 was 11693, checked in by nicolasmartin, 13 months ago

Macros renaming

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