1 | \documentclass[../tex_main/NEMO_manual]{subfiles} |
---|
2 | \begin{document} |
---|
3 | % ================================================================ |
---|
4 | % Appendix E : Note on some algorithms |
---|
5 | % ================================================================ |
---|
6 | \chapter{Note on some algorithms} |
---|
7 | \label{apdx:E} |
---|
8 | \minitoc |
---|
9 | |
---|
10 | \newpage |
---|
11 | $\ $\newline % force a new ligne |
---|
12 | |
---|
13 | This appendix some on going consideration on algorithms used or planned to be used |
---|
14 | in \NEMO. |
---|
15 | |
---|
16 | $\ $\newline % force a new ligne |
---|
17 | |
---|
18 | % ------------------------------------------------------------------------------------------------------------- |
---|
19 | % UBS scheme |
---|
20 | % ------------------------------------------------------------------------------------------------------------- |
---|
21 | \section{Upstream Biased Scheme (UBS) (\protect\np{ln\_traadv\_ubs}\forcode{ = .true.})} |
---|
22 | \label{sec:TRA_adv_ubs} |
---|
23 | |
---|
24 | The UBS advection scheme is an upstream biased third order scheme based on |
---|
25 | an upstream-biased parabolic interpolation. It is also known as Cell Averaged |
---|
26 | QUICK scheme (Quadratic Upstream Interpolation for Convective |
---|
27 | Kinematics). For example, in the $i$-direction : |
---|
28 | \begin{equation} \label{eq:tra_adv_ubs2} |
---|
29 | \tau _u^{ubs} = \left\{ \begin{aligned} |
---|
30 | & \tau _u^{cen4} + \frac{1}{12} \,\tau"_i & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ |
---|
31 | & \tau _u^{cen4} - \frac{1}{12} \,\tau"_{i+1} & \quad \text{if }\ u_{i+1/2} < 0 |
---|
32 | \end{aligned} \right. |
---|
33 | \end{equation} |
---|
34 | or equivalently, the advective flux is |
---|
35 | \begin{equation} \label{eq:tra_adv_ubs2} |
---|
36 | U_{i+1/2} \ \tau _u^{ubs} |
---|
37 | =U_{i+1/2} \ \overline{ T_i - \frac{1}{6}\,\tau"_i }^{\,i+1/2} |
---|
38 | - \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] |
---|
39 | \end{equation} |
---|
40 | where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and |
---|
41 | $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. |
---|
42 | By choosing this expression for $\tau "$ we consider a fourth order approximation |
---|
43 | of $\partial_i^2$ with a constant i-grid spacing ($\Delta i=1$). |
---|
44 | |
---|
45 | Alternative choice: introduce the scale factors: |
---|
46 | $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} }\delta _{i+1/2}[\tau] \right]$. |
---|
47 | |
---|
48 | |
---|
49 | This results in a dissipatively dominant (i.e. hyper-diffusive) truncation |
---|
50 | error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of the |
---|
51 | advection scheme is similar to that reported in \cite{Farrow1995}. |
---|
52 | It is a relatively good compromise between accuracy and smoothness. It is |
---|
53 | not a \emph{positive} scheme meaning false extrema are permitted but the |
---|
54 | amplitude of such are significantly reduced over the centred second order |
---|
55 | method. Nevertheless it is not recommended to apply it to a passive tracer |
---|
56 | that requires positivity. |
---|
57 | |
---|
58 | The intrinsic diffusion of UBS makes its use risky in the vertical direction |
---|
59 | where the control of artificial diapycnal fluxes is of paramount importance. |
---|
60 | It has therefore been preferred to evaluate the vertical flux using the TVD |
---|
61 | scheme when \np{ln\_traadv\_ubs}\forcode{ = .true.}. |
---|
62 | |
---|
63 | For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds |
---|
64 | to a second order centred scheme is evaluated using the \textit{now} velocity |
---|
65 | (centred in time) while the second term which is the diffusive part of the scheme, |
---|
66 | is evaluated using the \textit{before} velocity (forward in time. This is discussed |
---|
67 | by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. UBS and QUICK |
---|
68 | schemes only differ by one coefficient. Substituting 1/6 with 1/8 in |
---|
69 | (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. |
---|
70 | This option is not available through a namelist parameter, since the 1/6 |
---|
71 | coefficient is hard coded. Nevertheless it is quite easy to make the |
---|
72 | substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme |
---|
73 | |
---|
74 | NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can |
---|
75 | be controlled by vertical advection (not vertical diffusion which is usually |
---|
76 | solved using an implicit scheme). Computer time can be saved by using a |
---|
77 | time-splitting technique on vertical advection. This possibility have been |
---|
78 | implemented and validated in ORCA05-L301. It is not currently offered in the |
---|
79 | current reference version. |
---|
80 | |
---|
81 | NB 2 : In a forthcoming release four options will be proposed for the |
---|
82 | vertical component used in the UBS scheme. $\tau _w^{ubs}$ will be |
---|
83 | evaluated using either \textit{(a)} a centred $2^{nd}$ order scheme , |
---|
84 | or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative |
---|
85 | parabolic splines following \citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, |
---|
86 | or \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an |
---|
87 | eight-order accurate conventional scheme. |
---|
88 | |
---|
89 | NB 3 : It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows: |
---|
90 | \begin{equation} \label{eq:tra_adv_ubs2} |
---|
91 | \tau _u^{ubs} = \left\{ \begin{aligned} |
---|
92 | & \tau _u^{cen4} + \frac{1}{12} \tau"_i & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ |
---|
93 | & \tau _u^{cen4} - \frac{1}{12} \tau"_{i+1} & \quad \text{if }\ u_{i+1/2} < 0 |
---|
94 | \end{aligned} \right. |
---|
95 | \end{equation} |
---|
96 | or equivalently |
---|
97 | \begin{equation} \label{eq:tra_adv_ubs2} |
---|
98 | \begin{split} |
---|
99 | e_{2u} e_{3u}\,u_{i+1/2} \ \tau _u^{ubs} |
---|
100 | &= e_{2u} e_{3u}\,u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\tau"_i }^{\,i+1/2} \\ |
---|
101 | & - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] |
---|
102 | \end{split} |
---|
103 | \end{equation} |
---|
104 | \autoref{eq:tra_adv_ubs2} has several advantages. First it clearly evidence that |
---|
105 | the UBS scheme is based on the fourth order scheme to which is added an |
---|
106 | upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order |
---|
107 | part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order |
---|
108 | part as stated above using \autoref{eq:tra_adv_ubs}. Third, the diffusive term is |
---|
109 | in fact a biharmonic operator with a eddy coefficient with is simply proportional |
---|
110 | to the velocity. |
---|
111 | |
---|
112 | laplacian diffusion: |
---|
113 | \begin{equation} \label{eq:tra_ldf_lap} |
---|
114 | \begin{split} |
---|
115 | D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\; e_{3T} } &\left[ {\quad \delta _i |
---|
116 | \left[ {A_u^{lT} \frac{e_{2u} e_{3u} }{e_{1u} }\;\delta _{i+1/2} |
---|
117 | \left[ T \right]} \right]} \right. |
---|
118 | \\ |
---|
119 | &\ \left. {+\; \delta _j \left[ |
---|
120 | {A_v^{lT} \left( {\frac{e_{1v} e_{3v} }{e_{2v} }\;\delta _{j+1/2} \left[ T |
---|
121 | \right]} \right)} \right]\quad } \right] |
---|
122 | \end{split} |
---|
123 | \end{equation} |
---|
124 | |
---|
125 | bilaplacian: |
---|
126 | \begin{equation} \label{eq:tra_ldf_lap} |
---|
127 | \begin{split} |
---|
128 | D_T^{lT} =&-\frac{1}{e_{1T} \; e_{2T}\; e_{3T}} \\ |
---|
129 | & \delta _i \left[ \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta _{i+1/2} |
---|
130 | \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} |
---|
131 | \delta _i \left[ \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta _{i+1/2} |
---|
132 | [T] \right] \right] \right] |
---|
133 | \end{split} |
---|
134 | \end{equation} |
---|
135 | with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$, |
---|
136 | $i.e.$ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$ |
---|
137 | it comes : |
---|
138 | \begin{equation} \label{eq:tra_ldf_lap} |
---|
139 | \begin{split} |
---|
140 | D_T^{lT} =&-\frac{1}{12}\,\frac{1}{e_{1T} \; e_{2T}\; e_{3T}} \\ |
---|
141 | & \delta _i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta _{i+1/2} |
---|
142 | \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} |
---|
143 | \delta _i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta _{i+1/2} |
---|
144 | [T] \right] \right] \right] |
---|
145 | \end{split} |
---|
146 | \end{equation} |
---|
147 | if the velocity is uniform ($i.e.$ $|u|=cst$) then the diffusive flux is |
---|
148 | \begin{equation} \label{eq:tra_ldf_lap} |
---|
149 | \begin{split} |
---|
150 | F_u^{lT} = - \frac{1}{12} |
---|
151 | e_{2u}\,e_{3u}\,|u| \;\sqrt{ e_{1u}}\,\delta _{i+1/2} |
---|
152 | \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} |
---|
153 | \delta _i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}}\:\delta _{i+1/2} |
---|
154 | [T] \right] \right] |
---|
155 | \end{split} |
---|
156 | \end{equation} |
---|
157 | beurk.... reverte the logic: starting from the diffusive part of the advective flux it comes: |
---|
158 | |
---|
159 | \begin{equation} \label{eq:tra_adv_ubs2} |
---|
160 | \begin{split} |
---|
161 | F_u^{lT} |
---|
162 | &= - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] |
---|
163 | \end{split} |
---|
164 | \end{equation} |
---|
165 | if the velocity is uniform ($i.e.$ $|u|=cst$) and choosing $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right]$ |
---|
166 | |
---|
167 | sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$): |
---|
168 | \begin{equation} \label{eq:tra_adv_ubs2} |
---|
169 | \begin{split} |
---|
170 | F_u^{lT} |
---|
171 | &= - \frac{1}{12} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{e_{1T}^3\,|u|}{e_{1T}e_{2T}\,e_{3T}}\,\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right] \right] |
---|
172 | \end{split} |
---|
173 | \end{equation} |
---|
174 | which leads to ${A_T^{lT}}^2 = \frac{1}{12} {e_{1T}}^3\ \overline{|u|}^{\,i+1/2}$ |
---|
175 | |
---|
176 | sol 2 coefficient at u-point: split $|u|$ into $\sqrt{|u|}$ and $e_{1T}$ into $\sqrt{e_{1u}}$ |
---|
177 | \begin{equation} \label{eq:tra_adv_ubs2} |
---|
178 | \begin{split} |
---|
179 | F_u^{lT} |
---|
180 | &= - \frac{1}{12} {e_{1u}}^1 \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{1}{e_{2T}\,e_{3T}}\,\delta _i \left[ \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right] \right] \\ |
---|
181 | &= - \frac{1}{12} e_{1u} \sqrt{e_{1u}|u|\,} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{1}{e_{1T}\,e_{2T}\,e_{3T}}\,\delta _i \left[ e_{1u} \sqrt{e_{1u}|u|\,} \frac{e_{2u} e_{3u} }{e_{1u}} \delta _{i+1/2}[\tau] \right] \right] |
---|
182 | \end{split} |
---|
183 | \end{equation} |
---|
184 | which leads to ${A_u^{lT}} = \frac{1}{12} {e_{1u}}^3\ |u|$ |
---|
185 | |
---|
186 | |
---|
187 | % ------------------------------------------------------------------------------------------------------------- |
---|
188 | % Leap-Frog energetic |
---|
189 | % ------------------------------------------------------------------------------------------------------------- |
---|
190 | \section{Leapfrog energetic} |
---|
191 | \label{sec:LF} |
---|
192 | |
---|
193 | We adopt the following semi-discrete notation for time derivative. Given the values of a variable $q$ at successive time step, the time derivation and averaging operators at the mid time step are: |
---|
194 | \begin{subequations} \label{eq:dt_mt} |
---|
195 | \begin{align} |
---|
196 | \delta _{t+\rdt/2} [q] &= \ \ \, q^{t+\rdt} - q^{t} \\ |
---|
197 | \overline q^{\,t+\rdt/2} &= \left\{ q^{t+\rdt} + q^{t} \right\} \; / \; 2 |
---|
198 | \end{align} |
---|
199 | \end{subequations} |
---|
200 | As for space operator, the adjoint of the derivation and averaging time operators are |
---|
201 | $\delta_t^*=\delta_{t+\rdt/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$ |
---|
202 | , respectively. |
---|
203 | |
---|
204 | The Leap-frog time stepping given by \autoref{eq:DOM_nxt} can be defined as: |
---|
205 | \begin{equation} \label{eq:LF} |
---|
206 | \frac{\partial q}{\partial t} |
---|
207 | \equiv \frac{1}{\rdt} \overline{ \delta _{t+\rdt/2}[q]}^{\,t} |
---|
208 | = \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt} |
---|
209 | \end{equation} |
---|
210 | Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$, not $2\rdt$ |
---|
211 | as it can be found sometime in literature. |
---|
212 | The leap-Frog time stepping is a second order centered scheme. As such it respects |
---|
213 | the quadratic invariant in integral forms, $i.e.$ the following continuous property, |
---|
214 | \begin{equation} \label{eq:Energy} |
---|
215 | \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} |
---|
216 | =\int_{t_0}^{t_1} {\frac{1}{2}\, \frac{\partial q^2}{\partial t} \;dt} |
---|
217 | = \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) , |
---|
218 | \end{equation} |
---|
219 | is satisfied in discrete form. Indeed, |
---|
220 | \begin{equation} \begin{split} |
---|
221 | \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} |
---|
222 | &\equiv \sum\limits_{0}^{N} |
---|
223 | {\frac{1}{\rdt} q^t \ \overline{ \delta _{t+\rdt/2}[q]}^{\,t} \ \rdt} |
---|
224 | \equiv \sum\limits_{0}^{N} { q^t \ \overline{ \delta _{t+\rdt/2}[q]}^{\,t} } \\ |
---|
225 | &\equiv \sum\limits_{0}^{N} { \overline{q}^{\,t+\Delta/2}{ \delta _{t+\rdt/2}[q]}} |
---|
226 | \equiv \sum\limits_{0}^{N} { \frac{1}{2} \delta _{t+\rdt/2}[q^2] }\\ |
---|
227 | &\equiv \sum\limits_{0}^{N} { \frac{1}{2} \delta _{t+\rdt/2}[q^2] } |
---|
228 | \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) |
---|
229 | \end{split} \end{equation} |
---|
230 | NB here pb of boundary condition when applying the adjoin! In space, setting to 0 |
---|
231 | the quantity in land area is sufficient to get rid of the boundary condition |
---|
232 | (equivalently of the boundary value of the integration by part). In time this boundary |
---|
233 | condition is not physical and \textbf{add something here!!!} |
---|
234 | |
---|
235 | |
---|
236 | |
---|
237 | |
---|
238 | |
---|
239 | |
---|
240 | % ================================================================ |
---|
241 | % Iso-neutral diffusion : |
---|
242 | % ================================================================ |
---|
243 | |
---|
244 | \section{Lateral diffusion operator} |
---|
245 | |
---|
246 | % ================================================================ |
---|
247 | % Griffies' iso-neutral diffusion operator : |
---|
248 | % ================================================================ |
---|
249 | \subsection{Griffies iso-neutral diffusion operator} |
---|
250 | |
---|
251 | Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98} |
---|
252 | scheme, but is formulated within the \NEMO framework ($i.e.$ using scale |
---|
253 | factors rather than grid-size and having a position of $T$-points that is not |
---|
254 | necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}). |
---|
255 | |
---|
256 | In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO, |
---|
257 | the off-diagonal terms of the small angle diffusion tensor contain several double |
---|
258 | spatial averages of a gradient, for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$. |
---|
259 | It is apparent that the combination of a $k$ average and a $k$ derivative of the tracer |
---|
260 | allows for the presence of grid point oscillation structures that will be invisible |
---|
261 | to the operator. These structures are \textit{computational modes}. They |
---|
262 | will not be damped by the iso-neutral operator, and even possibly amplified by it. |
---|
263 | In other word, the operator applied to a tracer does not warranties the decrease of |
---|
264 | its global average variance. To circumvent this, we have introduced a smoothing of |
---|
265 | the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). Nevertheless, this technique |
---|
266 | works fine for $T$ and $S$ as they are active tracers ($i.e.$ they enter the computation |
---|
267 | of density), but it does not work for a passive tracer. \citep{Griffies_al_JPO98} introduce |
---|
268 | a different way to discretise the off-diagonal terms that nicely solve the problem. |
---|
269 | The idea is to get rid of combinations of an averaged in one direction combined |
---|
270 | with a derivative in the same direction by considering triads. For example in the |
---|
271 | (\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows: |
---|
272 | \begin{equation} \label{eq:Gf_triads} |
---|
273 | _i^k \mathbb{T}_{i_p}^{k_p} (T) |
---|
274 | = \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k} \ A_i^k \left( |
---|
275 | \frac{ \delta_{i + i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } |
---|
276 | -\ {_i^k \mathbb{R}_{i_p}^{k_p}} \ \frac{ \delta_{k+k_p} [T^i] }{ {e_{3w}}_{\,i}^{\,k+k_p} } |
---|
277 | \right) |
---|
278 | \end{equation} |
---|
279 | where the indices $i_p$ and $k_p$ define the four triads and take the following value: |
---|
280 | $i_p = -1/2$ or $1/2$ and $k_p = -1/2$ or $1/2$, |
---|
281 | $b_u= e_{1u}\,e_{2u}\,e_{3u}$ is the volume of $u$-cells, |
---|
282 | $A_i^k$ is the lateral eddy diffusivity coefficient defined at $T$-point, |
---|
283 | and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad : |
---|
284 | \begin{equation} \label{eq:Gf_slopes} |
---|
285 | _i^k \mathbb{R}_{i_p}^{k_p} |
---|
286 | =\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} \ \frac |
---|
287 | {\left(\alpha / \beta \right)_i^k \ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] } |
---|
288 | {\left(\alpha / \beta \right)_i^k \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] } |
---|
289 | \end{equation} |
---|
290 | Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of |
---|
291 | multiplying the temperature derivative by $\alpha$ and the salinity derivative |
---|
292 | by $\beta$. This is more efficient as the ratio $\alpha / \beta$ can to be |
---|
293 | evaluated directly. |
---|
294 | |
---|
295 | Note that in \autoref{eq:Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of |
---|
296 | ${b_{uw}}_{\,i+i_p}^{\,k+k_p}$. This choice has been motivated by the decrease |
---|
297 | of tracer variance and the presence of partial cell at the ocean bottom |
---|
298 | (see \autoref{apdx:Gf_operator}). |
---|
299 | |
---|
300 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
301 | \begin{figure}[!ht] \begin{center} |
---|
302 | \includegraphics[width=0.70\textwidth]{Fig_ISO_triad} |
---|
303 | \caption{ \protect\label{fig:ISO_triad} |
---|
304 | Triads used in the Griffies's like iso-neutral diffision scheme for |
---|
305 | $u$-component (upper panel) and $w$-component (lower panel).} |
---|
306 | \end{center} |
---|
307 | \end{figure} |
---|
308 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
309 | |
---|
310 | The four iso-neutral fluxes associated with the triads are defined at $T$-point. |
---|
311 | They take the following expression : |
---|
312 | \begin{flalign} \label{eq:Gf_fluxes} |
---|
313 | \begin{split} |
---|
314 | {_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) |
---|
315 | &= \ \; \qquad \quad { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+i_p}^{\,k}} \\ |
---|
316 | {_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T) |
---|
317 | &= -\; { _i^k \mathbb{R}_{i_p}^{k_p} } |
---|
318 | \ \; { _i^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+k_p}} |
---|
319 | \end{split} |
---|
320 | \end{flalign} |
---|
321 | |
---|
322 | The resulting iso-neutral fluxes at $u$- and $w$-points are then given by the |
---|
323 | sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}): |
---|
324 | \begin{flalign} \label{eq:iso_flux} |
---|
325 | \textbf{F}_{iso}(T) |
---|
326 | &\equiv \sum_{\substack{i_p,\,k_p}} |
---|
327 | \begin{pmatrix} |
---|
328 | {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) \\ |
---|
329 | \\ |
---|
330 | {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T) \\ |
---|
331 | \end{pmatrix} \notag \\ |
---|
332 | & \notag \\ |
---|
333 | &\equiv \sum_{\substack{i_p,\,k_p}} |
---|
334 | \begin{pmatrix} |
---|
335 | && { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{1u}}_{\,i+1/2}^{\,k} } \\ |
---|
336 | \\ |
---|
337 | & -\; { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } |
---|
338 | & {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T) \;\ / \ { {e_{3w}}_{\,i}^{\,k+1/2} } \\ |
---|
339 | \end{pmatrix} % \\ |
---|
340 | % &\\ |
---|
341 | % &\equiv \sum_{\substack{i_p,\,k_p}} |
---|
342 | % \begin{pmatrix} |
---|
343 | % \qquad \qquad \qquad |
---|
344 | % \frac{1}{ {e_{1u}}_{\,i+1/2}^{\,k} } \ \; |
---|
345 | % { _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} }(T)\\ |
---|
346 | % \\ |
---|
347 | % -\frac{1}{ {e_{3w}}_{\,i}^{\,k+1/2} } \ \; |
---|
348 | % { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \; |
---|
349 | % {_i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} }(T)\\ |
---|
350 | % \end{pmatrix} |
---|
351 | \end{flalign} |
---|
352 | resulting in a iso-neutral diffusion tendency on temperature given by the divergence |
---|
353 | of the sum of all the four triad fluxes : |
---|
354 | \begin{equation} \label{eq:Gf_operator} |
---|
355 | D_l^T = \frac{1}{b_T} \sum_{\substack{i_p,\,k_p}} \left\{ |
---|
356 | \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] |
---|
357 | + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \right\} |
---|
358 | \end{equation} |
---|
359 | where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. |
---|
360 | |
---|
361 | This expression of the iso-neutral diffusion has been chosen in order to satisfy |
---|
362 | the following six properties: |
---|
363 | \begin{description} |
---|
364 | \item[$\bullet$ horizontal diffusion] The discretization of the diffusion operator |
---|
365 | recovers the traditional five-point Laplacian in the limit of flat iso-neutral direction : |
---|
366 | \begin{equation} \label{eq:Gf_property1a} |
---|
367 | D_l^T = \frac{1}{b_T} \ \delta_{i} |
---|
368 | \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right] |
---|
369 | \qquad \text{when} \quad |
---|
370 | { _i^k \mathbb{R}_{i_p}^{k_p} }=0 |
---|
371 | \end{equation} |
---|
372 | |
---|
373 | \item[$\bullet$ implicit treatment in the vertical] In the diagonal term associated |
---|
374 | with the vertical divergence of the iso-neutral fluxes (i.e. the term associated |
---|
375 | with a second order vertical derivative) appears only tracer values associated |
---|
376 | with a single water column. This is of paramount importance since it means |
---|
377 | that the implicit in time algorithm for solving the vertical diffusion equation can |
---|
378 | be used to evaluate this term. It is a necessity since the vertical eddy diffusivity |
---|
379 | associated with this term, |
---|
380 | \begin{equation} |
---|
381 | \sum_{\substack{i_p, \,k_p}} \left\{ |
---|
382 | A_i^k \; \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 |
---|
383 | \right\} |
---|
384 | \end{equation} |
---|
385 | can be quite large. |
---|
386 | |
---|
387 | \item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of locally referenced |
---|
388 | potential density is zero, $i.e.$ |
---|
389 | \begin{align} \label{eq:Gf_property2} |
---|
390 | \begin{matrix} |
---|
391 | &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} (\rho)} |
---|
392 | &= &\alpha_i^k &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) |
---|
393 | &- \ \; \beta _i^k &{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (S) & = \ 0 \\ |
---|
394 | &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)} |
---|
395 | &= &\alpha_i^k &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (T) |
---|
396 | &- \ \; \beta _i^k &{_i^k {\mathbb{F}_w}_{i_p}^{k_p} } (S) &= \ 0 |
---|
397 | \end{matrix} |
---|
398 | \end{align} |
---|
399 | This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$ |
---|
400 | and the definition of the triads' slopes \autoref{eq:Gf_slopes}. |
---|
401 | |
---|
402 | \item[$\bullet$ conservation of tracer] The iso-neutral diffusion term conserve the |
---|
403 | total tracer content, $i.e.$ |
---|
404 | \begin{equation} \label{eq:Gf_property1} |
---|
405 | \sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0 |
---|
406 | \end{equation} |
---|
407 | This property is trivially satisfied since the iso-neutral diffusive operator |
---|
408 | is written in flux form. |
---|
409 | |
---|
410 | \item[$\bullet$ decrease of tracer variance] The iso-neutral diffusion term does |
---|
411 | not increase the total tracer variance, $i.e.$ |
---|
412 | \begin{equation} \label{eq:Gf_property1} |
---|
413 | \sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0 |
---|
414 | \end{equation} |
---|
415 | The property is demonstrated in the \autoref{apdx:Gf_operator}. It is a |
---|
416 | key property for a diffusion term. It means that the operator is also a dissipation |
---|
417 | term, $i.e.$ it is a sink term for the square of the quantity on which it is applied. |
---|
418 | It therfore ensure that, when the diffusivity coefficient is large enough, the field |
---|
419 | on which it is applied become free of grid-point noise. |
---|
420 | |
---|
421 | \item[$\bullet$ self-adjoint operator] The iso-neutral diffusion operator is self-adjoint, |
---|
422 | $i.e.$ |
---|
423 | \begin{equation} \label{eq:Gf_property1} |
---|
424 | \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} |
---|
425 | \end{equation} |
---|
426 | In other word, there is no needs to develop a specific routine from the adjoint of this |
---|
427 | operator. We just have to apply the same routine. This properties can be demonstrated |
---|
428 | quite easily in a similar way the "non increase of tracer variance" property has been |
---|
429 | proved (see \autoref{apdx:Gf_operator}). |
---|
430 | \end{description} |
---|
431 | |
---|
432 | |
---|
433 | $\ $\newline %force an empty line |
---|
434 | % ================================================================ |
---|
435 | % Skew flux formulation for Eddy Induced Velocity : |
---|
436 | % ================================================================ |
---|
437 | \subsection{Eddy induced velocity and skew flux formulation} |
---|
438 | |
---|
439 | When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined), |
---|
440 | an additional advection term is added. The associated velocity is the so called |
---|
441 | eddy induced velocity, the formulation of which depends on the slopes of iso- |
---|
442 | neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used |
---|
443 | here are referenced to the geopotential surfaces, $i.e.$ \autoref{eq:ldfslp_geo} |
---|
444 | is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo} |
---|
445 | + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. |
---|
446 | |
---|
447 | The eddy induced velocity is given by: |
---|
448 | \begin{equation} \label{eq:eiv_v} |
---|
449 | \begin{split} |
---|
450 | u^* & = - \frac{1}{e_2\,e_{3}} \;\partial_k \left( e_2 \, A_e \; r_i \right) |
---|
451 | = - \frac{1}{e_3} \;\partial_k \left( A_e \; r_i \right) \\ |
---|
452 | v^* & = - \frac{1}{e_1\,e_3}\; \partial_k \left( e_1 \, A_e \; r_j \right) |
---|
453 | = - \frac{1}{e_3} \;\partial_k \left( A_e \; r_j \right) \\ |
---|
454 | w^* & = \frac{1}{e_1\,e_2}\; \left\{ \partial_i \left( e_2 \, A_e \; r_i \right) |
---|
455 | + \partial_j \left( e_1 \, A_e \;r_j \right) \right\} \\ |
---|
456 | \end{split} |
---|
457 | \end{equation} |
---|
458 | where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the |
---|
459 | slopes between the iso-neutral and the geopotential surfaces. |
---|
460 | %%gm wrong: to be modified with 2 2D streamfunctions |
---|
461 | In other words, |
---|
462 | the eddy induced velocity can be derived from a vector streamfuntion, $\phi$, which |
---|
463 | is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^* = \textbf{k} \times \nabla \phi$ |
---|
464 | %%end gm |
---|
465 | |
---|
466 | A traditional way to implement this additional advection is to add it to the eulerian |
---|
467 | velocity prior to compute the tracer advection. This allows us to take advantage of |
---|
468 | all the advection schemes offered for the tracers (see \autoref{sec:TRA_adv}) and not just |
---|
469 | a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers |
---|
470 | where \emph{positivity} of the advection scheme is of paramount importance. |
---|
471 | % give here the expression using the triads. It is different from the one given in \autoref{eq:ldfeiv} |
---|
472 | % see just below a copy of this equation: |
---|
473 | %\begin{equation} \label{eq:ldfeiv} |
---|
474 | %\begin{split} |
---|
475 | % u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\ |
---|
476 | % v^* & = \frac{1}{e_{1u}e_{3v}}\; \delta_k \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right]\\ |
---|
477 | %w^* & = \frac{1}{e_{1w}e_{2w}}\; \left\{ \delta_i \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right] + %\delta_j \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right] \right\} \\ |
---|
478 | %\end{split} |
---|
479 | %\end{equation} |
---|
480 | \begin{equation} \label{eq:eiv_vd} |
---|
481 | \textbf{F}_{eiv}^T \equiv \left( \begin{aligned} |
---|
482 | \sum_{\substack{i_p,\,k_p}} & |
---|
483 | +{e_{2u}}_{i+1/2-i_p}^{k} \ \ {A_{e}}_{i+1/2-i_p}^{k} |
---|
484 | \ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} } \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ |
---|
485 | \\ |
---|
486 | \sum_{\substack{i_p,\,k_p}} & |
---|
487 | - {e_{2u}}_i^{k+1/2-k_p} \ {A_{e}}_i^{k+1/2-k_p} |
---|
488 | \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \delta_{i+i_p}[T^{k+1/2-k_p}] \\ |
---|
489 | \end{aligned} \right) |
---|
490 | \end{equation} |
---|
491 | |
---|
492 | \citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, |
---|
493 | the so-called skew form. It is based on a transformation of the advective fluxes |
---|
494 | using the non-divergent nature of the eddy induced velocity. |
---|
495 | For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be |
---|
496 | transformed as follows: |
---|
497 | \begin{flalign*} |
---|
498 | \begin{split} |
---|
499 | \textbf{F}_{eiv}^T = |
---|
500 | \begin{pmatrix} |
---|
501 | {e_{2}\,e_{3}\; u^*} \\ |
---|
502 | {e_{1}\,e_{2}\; w^*} \\ |
---|
503 | \end{pmatrix} \; T |
---|
504 | &= |
---|
505 | \begin{pmatrix} |
---|
506 | { - \partial_k \left( e_{2} \, A_{e} \; r_i \right) \; T \;} \\ |
---|
507 | {+ \partial_i \left( e_{2} \, A_{e} \; r_i \right) \; T \;} \\ |
---|
508 | \end{pmatrix} \\ |
---|
509 | &= |
---|
510 | \begin{pmatrix} |
---|
511 | { - \partial_k \left( e_{2} \, A_{e} \; r_i \; T \right) \;} \\ |
---|
512 | {+ \partial_i \left( e_{2} \, A_{e} \; r_i \; T \right) \;} \\ |
---|
513 | \end{pmatrix} |
---|
514 | + |
---|
515 | \begin{pmatrix} |
---|
516 | {+ e_{2} \, A_{e} \; r_i \; \partial_k T} \\ |
---|
517 | { - e_{2} \, A_{e} \; r_i \; \partial_i T} \\ |
---|
518 | \end{pmatrix} |
---|
519 | \end{split} |
---|
520 | \end{flalign*} |
---|
521 | and since the eddy induces velocity field is no-divergent, we end up with the skew |
---|
522 | form of the eddy induced advective fluxes: |
---|
523 | \begin{equation} \label{eq:eiv_skew_continuous} |
---|
524 | \textbf{F}_{eiv}^T = \begin{pmatrix} |
---|
525 | {+ e_{2} \, A_{e} \; r_i \; \partial_k T} \\ |
---|
526 | { - e_{2} \, A_{e} \; r_i \; \partial_i T} \\ |
---|
527 | \end{pmatrix} |
---|
528 | \end{equation} |
---|
529 | The tendency associated with eddy induced velocity is then simply the divergence |
---|
530 | of the \autoref{eq:eiv_skew_continuous} fluxes. It naturally conserves the tracer |
---|
531 | content, as it is expressed in flux form and, as the advective form, it preserve the |
---|
532 | tracer variance. Another interesting property of \autoref{eq:eiv_skew_continuous} |
---|
533 | form is that when $A=A_e$, a simplification occurs in the sum of the iso-neutral |
---|
534 | diffusion and eddy induced velocity terms: |
---|
535 | \begin{flalign} \label{eq:eiv_skew+eiv_continuous} |
---|
536 | \textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &= |
---|
537 | \begin{pmatrix} |
---|
538 | + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T - e_2 \, A \; r_i \;\partial_k T \\ |
---|
539 | - e_2 \, A_{e} \; r_i \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T \\ |
---|
540 | \end{pmatrix} |
---|
541 | + |
---|
542 | \begin{pmatrix} |
---|
543 | {+ e_{2} \, A_{e} \; r_i \; \partial_k T} \\ |
---|
544 | { - e_{2} \, A_{e} \; r_i \; \partial_i T} \\ |
---|
545 | \end{pmatrix} \\ |
---|
546 | &= \begin{pmatrix} |
---|
547 | + \frac{e_2\,e_3\,}{e_1} A \;\partial_i T \\ |
---|
548 | - 2\; e_2 \, A_{e} \; r_i \;\partial_i T + \frac{e_1\,e_2}{e_3} \, A \; r_i^2 \;\partial_k T \\ |
---|
549 | \end{pmatrix} |
---|
550 | \end{flalign} |
---|
551 | The horizontal component reduces to the one use for an horizontal laplacian |
---|
552 | operator and the vertical one keep the same complexity, but not more. This property |
---|
553 | has been used to reduce the computational time \citep{Griffies_JPO98}, but it is |
---|
554 | not of practical use as usually $A \neq A_e$. Nevertheless this property can be used to |
---|
555 | choose a discret form of \autoref{eq:eiv_skew_continuous} which is consistent with the |
---|
556 | iso-neutral operator \autoref{eq:Gf_operator}. Using the slopes \autoref{eq:Gf_slopes} |
---|
557 | and defining $A_e$ at $T$-point($i.e.$ as $A$, the eddy diffusivity coefficient), |
---|
558 | the resulting discret form is given by: |
---|
559 | \begin{equation} \label{eq:eiv_skew} |
---|
560 | \textbf{F}_{eiv}^T \equiv \frac{1}{4} \left( \begin{aligned} |
---|
561 | \sum_{\substack{i_p,\,k_p}} & |
---|
562 | +{e_{2u}}_{i+1/2-i_p}^{k} \ \ {A_{e}}_{i+1/2-i_p}^{k} |
---|
563 | \ \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} } \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ |
---|
564 | \\ |
---|
565 | \sum_{\substack{i_p,\,k_p}} & |
---|
566 | - {e_{2u}}_i^{k+1/2-k_p} \ {A_{e}}_i^{k+1/2-k_p} |
---|
567 | \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \delta_{i+i_p}[T^{k+1/2-k_p}] \\ |
---|
568 | \end{aligned} \right) |
---|
569 | \end{equation} |
---|
570 | Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells. |
---|
571 | In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces |
---|
572 | must be added to $\mathbb{R}$ for the discret form to be exact. |
---|
573 | |
---|
574 | Such a choice of discretisation is consistent with the iso-neutral operator as it uses the |
---|
575 | same definition for the slopes. It also ensures the conservation of the tracer variance |
---|
576 | (see Appendix \autoref{apdx:eiv_skew}), $i.e.$ it does not include a diffusive component |
---|
577 | but is a "pure" advection term. |
---|
578 | |
---|
579 | |
---|
580 | |
---|
581 | |
---|
582 | $\ $\newpage %force an empty line |
---|
583 | % ================================================================ |
---|
584 | % Discrete Invariants of the iso-neutral diffrusion |
---|
585 | % ================================================================ |
---|
586 | \subsection{Discrete invariants of the iso-neutral diffrusion} |
---|
587 | \label{subsec:Gf_operator} |
---|
588 | |
---|
589 | Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane. |
---|
590 | |
---|
591 | This part will be moved in an Appendix. |
---|
592 | |
---|
593 | The continuous property to be demonstrated is : |
---|
594 | \begin{align*} |
---|
595 | \int_D D_l^T \; T \;dv \leq 0 |
---|
596 | \end{align*} |
---|
597 | The discrete form of its left hand side is obtained using \autoref{eq:iso_flux} |
---|
598 | |
---|
599 | \begin{align*} |
---|
600 | &\int_D D_l^T \; T \;dv \equiv \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\} \\ |
---|
601 | &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ |
---|
602 | \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] |
---|
603 | + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \ T \right\} \\ |
---|
604 | &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ |
---|
605 | {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T] |
---|
606 | + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \ \delta_{k+1/2} [T] \right\} \\ |
---|
607 | &\equiv -\sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ |
---|
608 | \frac{ _{i+1/2-i_p}^k \mathbb{T}_{i_p}^{k_p} (T) }{ {e_{1u}}_{\,i+1/2}^{\,k} } \ \delta_{i+1/2} [T] |
---|
609 | - { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \; |
---|
610 | \frac{ _i^{k+1/2-k_p} \mathbb{T}_{i_p}^{k_p} (T) }{ {e_{3w}}_{\,i}^{\,k+1/2} } \ \delta_{k+1/2} [T] |
---|
611 | \right\} \\ |
---|
612 | % |
---|
613 | \allowdisplaybreaks |
---|
614 | \intertext{ Expending the summation on $i_p$ and $k_p$, it becomes:} |
---|
615 | % |
---|
616 | &\equiv -\sum_{i,k} |
---|
617 | \begin{Bmatrix} |
---|
618 | &\ \ \Bigl( { _{i+1}^{k} \mathbb{T}_{-1/2}^{-1/2} (T) } |
---|
619 | &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} |
---|
620 | & -\ \ {_{i}^{k+1} \mathbb{R}_{-1/2}^{-1/2}} |
---|
621 | & {_{i}^{k+1} \mathbb{T}_{-1/2}^{-1/2} (T) } |
---|
622 | &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}} \Bigr) |
---|
623 | & \\ |
---|
624 | &+\Bigl( \ \;\; { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) } |
---|
625 | &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} |
---|
626 | & -\ \ {_i^{k+1} \mathbb{R}_{+1/2}^{-1/2}} |
---|
627 | & { _i^{k+1} \mathbb{T}_{+1/2}^{-1/2} (T) } |
---|
628 | &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}} \Bigr) |
---|
629 | & \\ |
---|
630 | &+\Bigl( { _{i+1}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) } |
---|
631 | &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} |
---|
632 | & -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{-1/2}^{+1/2}} |
---|
633 | & \ \;\;{_{i}^{k} \mathbb{T}_{-1/2}^{+1/2} (T) } |
---|
634 | &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}} \Bigr) |
---|
635 | & \\ |
---|
636 | &+\Bigl( \ \;\; { _{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) } |
---|
637 | &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} |
---|
638 | & -\ \ \ \;\;{_{i}^{k} \mathbb{R}_{+1/2}^{+1/2}} |
---|
639 | & \ \;\;{_{i}^{k} \mathbb{T}_{+1/2}^{+1/2} (T) } |
---|
640 | &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}} \Bigr) \\ |
---|
641 | \end{Bmatrix} |
---|
642 | % |
---|
643 | \allowdisplaybreaks |
---|
644 | \intertext{The summation is done over all $i$ and $k$ indices, it is therefore possible to introduce a shift of $-1$ either in $i$ or $k$ direction in order to regroup all the terms of the summation by triad at a ($i$,$k$) point. In other words, we regroup all the terms in the neighbourhood that contain a triad at the same ($i$,$k$) indices. It becomes: } |
---|
645 | % |
---|
646 | &\equiv -\sum_{i,k} |
---|
647 | \begin{Bmatrix} |
---|
648 | &\ \ \Bigl( {_i^k \mathbb{T}_{-1/2}^{-1/2} (T) } |
---|
649 | &\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} |
---|
650 | & -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}} |
---|
651 | & {_i^k \mathbb{T}_{-1/2}^{-1/2} (T) } |
---|
652 | &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}} \Bigr) |
---|
653 | & \\ |
---|
654 | &+\Bigl( { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) } |
---|
655 | &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} |
---|
656 | & -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}} |
---|
657 | & { _i^k \mathbb{T}_{+1/2}^{-1/2} (T) } |
---|
658 | &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}} \Bigr) |
---|
659 | & \\ |
---|
660 | &+\Bigl( {_i^k \mathbb{T}_{-1/2}^{+1/2} (T) } |
---|
661 | &\frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} |
---|
662 | & -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}} |
---|
663 | & {_i^k \mathbb{T}_{-1/2}^{+1/2} (T) } |
---|
664 | &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}} \Bigr) |
---|
665 | & \\ |
---|
666 | &+\Bigl( { _i^k \mathbb{T}_{+1/2}^{+1/2} (T) } |
---|
667 | &\frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} |
---|
668 | & -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}} |
---|
669 | & {_i^k \mathbb{T}_{+1/2}^{+1/2} (T) } |
---|
670 | &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}} \Bigr) \\ |
---|
671 | \end{Bmatrix} \\ |
---|
672 | % |
---|
673 | \allowdisplaybreaks |
---|
674 | \intertext{Then outing in factor the triad in each of the four terms of the summation and substituting the triads by their expression given in \autoref{eq:Gf_triads}. It becomes: } |
---|
675 | % |
---|
676 | &\equiv -\sum_{i,k} |
---|
677 | \begin{Bmatrix} |
---|
678 | &\ \ \Bigl( \frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} |
---|
679 | & -\ \ {_i^k \mathbb{R}_{-1/2}^{-1/2}} |
---|
680 | &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}} \Bigr)^2 |
---|
681 | & \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k} \ A_i^k |
---|
682 | & \\ |
---|
683 | &+\Bigl( \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} |
---|
684 | & -\ \ {_i^k \mathbb{R}_{+1/2}^{-1/2}} |
---|
685 | &\frac{ \delta_{k-1/2} [T] }{{e_{3w}}_{\,i}^{\,k-1/2}} \Bigr)^2 |
---|
686 | & \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k} \ A_i^k |
---|
687 | & \\ |
---|
688 | &+\Bigl( \frac{ \delta_{i -1/2} [T] }{{e_{1u} }_{\,i-1/2}^{\,k}} |
---|
689 | & -\ \ {_i^k \mathbb{R}_{-1/2}^{+1/2}} |
---|
690 | &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}} \Bigr)^2 |
---|
691 | & \frac{1}{4} \ {b_u}_{\,i-1/2}^{\,k} \ A_i^k |
---|
692 | & \\ |
---|
693 | &+\Bigl( \frac{ \delta_{i +1/2} [T] }{{e_{1u} }_{\,i+1/2}^{\,k}} |
---|
694 | & -\ \ {_i^k \mathbb{R}_{+1/2}^{+1/2}} |
---|
695 | &\frac{ \delta_{k+1/2} [T] }{{e_{3w}}_{\,i}^{\,k+1/2}} \Bigr)^2 |
---|
696 | & \frac{1}{4} \ {b_u}_{\,i+1/2}^{\,k} \ A_i^k \\ |
---|
697 | \end{Bmatrix} \\ |
---|
698 | & \\ |
---|
699 | % |
---|
700 | &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ |
---|
701 | \begin{matrix} |
---|
702 | &\Bigl( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}} |
---|
703 | & -\ \ {_i^k \mathbb{R}_{i_p}^{k_p}} |
---|
704 | &\frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}} \Bigr)^2 |
---|
705 | & \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k} \ A_i^k \ \ |
---|
706 | \end{matrix} |
---|
707 | \right\} |
---|
708 | \quad \leq 0 |
---|
709 | \end{align*} |
---|
710 | The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities. |
---|
711 | |
---|
712 | Note that, if instead of multiplying $D_l^T$ by $T$, we were using another tracer field, let say $S$, then the previous demonstration would have let to: |
---|
713 | \begin{align*} |
---|
714 | \int_D S \; D_l^T \;dv &\equiv \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\} \\ |
---|
715 | &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ |
---|
716 | \left( \frac{ \delta_{i +i_p} [S] }{{e_{1u} }_{\,i+i_p}^{\,k}} |
---|
717 | - {_i^k \mathbb{R}_{i_p}^{k_p}} |
---|
718 | \frac{ \delta_{k+k_p} [S] }{{e_{3w}}_{\,i}^{\,k+k_p}} \right) \right. |
---|
719 | \\ & \qquad \qquad \qquad \ \left. |
---|
720 | \left( \frac{ \delta_{i +i_p} [T] }{{e_{1u} }_{\,i+i_p}^{\,k}} |
---|
721 | - {_i^k \mathbb{R}_{i_p}^{k_p}} |
---|
722 | \frac{ \delta_{k+k_p} [T] }{{e_{3w}}_{\,i}^{\,k+k_p}} \right) |
---|
723 | \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k} \ A_i^k \ |
---|
724 | \right\} |
---|
725 | % |
---|
726 | \allowdisplaybreaks |
---|
727 | \intertext{which, by applying the same operation as before but in reverse order, leads to: } |
---|
728 | % |
---|
729 | &\equiv \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\} |
---|
730 | \end{align*} |
---|
731 | This means that the iso-neutral operator is self-adjoint. There is no need to develop a specific to obtain it. |
---|
732 | |
---|
733 | |
---|
734 | |
---|
735 | $\ $\newpage %force an empty line |
---|
736 | % ================================================================ |
---|
737 | % Discrete Invariants of the skew flux formulation |
---|
738 | % ================================================================ |
---|
739 | \subsection{Discrete invariants of the skew flux formulation} |
---|
740 | \label{subsec:eiv_skew} |
---|
741 | |
---|
742 | |
---|
743 | Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane. |
---|
744 | |
---|
745 | This have to be moved in an Appendix. |
---|
746 | |
---|
747 | The continuous property to be demonstrated is : |
---|
748 | \begin{align*} |
---|
749 | \int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv \equiv 0 |
---|
750 | \end{align*} |
---|
751 | The discrete form of its left hand side is obtained using \autoref{eq:eiv_skew} |
---|
752 | \begin{align*} |
---|
753 | \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}} \Biggl\{ \;\; |
---|
754 | \delta_i &\left[ |
---|
755 | {e_{2u}}_{i+i_p+1/2}^{k} \;\ \ {A_{e}}_{i+i_p+1/2}^{k} |
---|
756 | \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} } \quad \delta_{k+k_p}[T_{i+i_p+1/2}] |
---|
757 | \right] \; T_i^k \\ |
---|
758 | - \delta_k &\left[ |
---|
759 | {e_{2u}}_i^{k+k_p+1/2} \ \ {A_{e}}_i^{k+k_p+1/2} |
---|
760 | \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} } \ \ \delta_{i+i_p}[T^{k+k_p+1/2}] |
---|
761 | \right] \; T_i^k \ \Biggr\} |
---|
762 | \end{align*} |
---|
763 | apply the adjoint of delta operator, it becomes |
---|
764 | \begin{align*} |
---|
765 | \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}} \Biggl\{ \;\; |
---|
766 | &\left( |
---|
767 | {e_{2u}}_{i+i_p+1/2}^{k} \;\ \ {A_{e}}_{i+i_p+1/2}^{k} |
---|
768 | \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} } \quad \delta_{k+k_p}[T_{i+i_p+1/2}] |
---|
769 | \right) \; \delta_{i+1/2}[T^{k}] \\ |
---|
770 | - &\left( |
---|
771 | {e_{2u}}_i^{k+k_p+1/2} \ \ {A_{e}}_i^{k+k_p+1/2} |
---|
772 | \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} } \ \ \delta_{i+i_p}[T^{k+k_p+1/2}] |
---|
773 | \right) \; \delta_{k+1/2}[T_{i}] \ \Biggr\} |
---|
774 | \end{align*} |
---|
775 | Expending the summation on $i_p$ and $k_p$, it becomes: |
---|
776 | \begin{align*} |
---|
777 | \begin{matrix} |
---|
778 | &\sum\limits_{i,k} \Bigl\{ |
---|
779 | &+{e_{2u}}_{i+1}^{k} &{A_{e}}_{i+1 }^{k} |
---|
780 | &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}] &\delta_{i+1/2}[T^{k}] &\\ |
---|
781 | &&+{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_{i}^{k\ \ \ \:} |
---|
782 | &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}} &\delta_{k-1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ |
---|
783 | &&+{e_{2u}}_{i+1}^{k} &{A_{e}}_{i+1 }^{k} |
---|
784 | &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}] &\delta_{i+1/2}[T^{k}] &\\ |
---|
785 | &&+{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_{i}^{k\ \ \ \:} |
---|
786 | &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ |
---|
787 | % |
---|
788 | &&-{e_{2u}}_i^{k+1} &{A_{e}}_i^{k+1} |
---|
789 | &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}} &\delta_{i-1/2}[T^{k+1}] &\delta_{k+1/2}[T_{i}] &\\ |
---|
790 | &&-{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_i^{k\ \ \ \:} |
---|
791 | &{\ \ \;_i^k \mathbb{R}_{-1/2}^{+1/2}} &\delta_{i-1/2}[T^{k\ \ \ \:}] &\delta_{k+1/2}[T_{i}] &\\ |
---|
792 | &&-{e_{2u}}_i^{k+1 } &{A_{e}}_i^{k+1} |
---|
793 | &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}} &\delta_{i+1/2}[T^{k+1}] &\delta_{k+1/2}[T_{i}] &\\ |
---|
794 | &&-{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_i^{k\ \ \ \:} |
---|
795 | &{\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{i+1/2}[T^{k\ \ \ \:}] &\delta_{k+1/2}[T_{i}] |
---|
796 | &\Bigr\} \\ |
---|
797 | \end{matrix} |
---|
798 | \end{align*} |
---|
799 | The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the |
---|
800 | same but of opposite signs, they cancel out. |
---|
801 | Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$. |
---|
802 | The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the |
---|
803 | same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$ |
---|
804 | they cancel out with the neighbouring grid points. |
---|
805 | Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the |
---|
806 | $i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the |
---|
807 | tracer is preserved by the discretisation of the skew fluxes. |
---|
808 | |
---|
809 | \end{document} |
---|