New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
annex_C.tex in NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/doc/tex_sub – NEMO

source: NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/doc/tex_sub/annex_C.tex @ 9979

Last change on this file since 9979 was 9957, checked in by acc, 6 years ago

New development branch for the adaptive-implicit vertical advection (05_Gurvan-Vertical_advection)

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