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