The aim of this chapter is to derive how differential and algebraic operators and boundary conditions are to be discretized. To do this, consider the atmosphere $A \subseteq \mathbb {R}^3$ as an open set and use the reversible equations of dry air as the system of equations, with the boundary conditions
\[ \begin{align} \mathbf{v}\cdot\mathbf{n} = 0 \end{align} \]
at the top and
\[ \begin{align} \mathbf{v} = \mathbf{0} \end{align} \]
at the bottom.
The following are chosen as prognostic variables:
the mass density $\rho$ in the centers of the boxes,
the potential temperature density $\newtilde{\theta}$ in the centers of the boxes,
the horizontal velocity $u$ at the edges of the cells
the vertical velocity $w$ in the levels as well
the temperature $T$ in the centers of the boxes as a semi-prognostic variable
The Hamilton function $H$ is then written as a function of $\rho, \newtilde{s}$ and $\mathbf{v}$,
\[ \begin{align} H = H\left(\rho, \newtilde{\theta}, \mathbf{v}\right). \end{align} \]
The prognostic equations for this case were derived in Sect. 10.1.2; they read
\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} &= -c^{(p)}\theta\nabla\Pi - \nabla k + \mathbf{v}\times\etabi - \nabla\phi,\tag{34.4}\label{eq:prog_rev_0}\\ \frac{\partial\rho}{\partial t} &= -\nabla\cdot\left(\rho\mathbf{v}\right),\tag{34.5}\label{eq:prog_rev_1}\\ \frac{\partial\newtilde{\theta}}{\partial t} &= -\nabla\cdot\left(\newtilde{\theta}\mathbf{v}\right),\tag{34.6}\label{eq:prog_rev_2}\\ T &= f\left(\rho, \newtilde{\theta}\right).\tag{34.7}\label{eq:prog_rev_3} \end{align} \]
The function $f$ in Eq. (34.7) is the thermal equation of state of ideal gases, $\etabi\coloneqq\zetabi + \mathbf{f}$ is the absolute vorticity.
The Poisson brackets are global integrals. One therefore defines a discrete integral operator according to the substitution
\[ \begin{align} \int_A\psi\left(\mathbf{r}\right)d^3r \to \sum_{i = 1}^N\psi_iV_i, \end{align} \]
here $\psi\left(\mathbf{r}\right)$ is an integrable scalar or vector field, $\psi_i$ is the value of the discretized version of $\psi$ in grid box $i$, $V_i$ is the volume of grid box $i$ and $N$ is the number of grid boxes. The Poisson brackets are discretized according to the basic ansatz
\[ \begin{align} \frac{dF\left[\mathbf{u}\right]}{dt} = \frac{d}{dt}\int_Af\left(\mathbf{u}\left(t\right)\right)d^3r = \int_A\frac{\delta F}{\delta\mathbf{u}}\cdot\frac{\partial\mathbf{u}}{\partial t}d^3r \to \sum_{c, k}\left(\frac{\delta F}{\delta\mathbf{u}}\cdot\frac{\partial\mathbf{u}}{\partial t}\right)_{c, k}V_{c, k} \end{align} \]
First, the Hamilton functional is needed. One discretizes Eq. (10.81) in the form
\[ \begin{align} H\left(\rho, \newtilde{s}, \mathbf{v}\right) = \sum_{c, k}\left[\frac{1}{2}\rho_{c, k}\left(\mathbf{v}\cdot\mathbf{v}\right)_{c, k} + \rho_{c, k}\phi_{c, k} + \newtilde{I}_{c, k}\right]V_{c, k}. \end{align} \]
This is divided into two parts
\[ \begin{align} H_\text{kin}\left(\mathbf{v}\right) &= \sum_{c, k}\frac{1}{2}\rho_{c, k}\left(\mathbf{v}\cdot\mathbf{v}\right)_{c, k}V_{c, k},\\ H_{I'}\left(\rho, \newtilde{s}\right) &= \sum_{c, k}\left(\rho_{c, k}\phi_{c, k} + \newtilde{I}_{c, k}\right)V_{c, k}. \end{align} \]
The kinetic-energy part can be split even further:
\[ \begin{align} H_\text{kin}^{(u)}\left(u\right) &= \sum_{c, k}\frac{1}{2}\rho_{c, k}\sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Delta z_{e, k}u_{e, k}^2V_{c, k} = \sum_{c, k}\sum_{e\in c}\frac{\rho_{c, k}l_ed_e\Delta z_{e, k}}{4}u_{e, k}^2,\tag{34.13}\label{eq:h_kin_u}\\ H_\text{kin}^{(w)}\left(w\right) &= \sum_{c, k}\frac{1}{2}\rho_{c, k}\left(\frac{\Delta z_{c, k - 1/2}}{2\Delta z_{c, k}}w_u^2 + \frac{\Delta z_{c, k + 1/2}}{2\Delta z_{c, k}}w_l^2\right)V_{c, k}\nonumber\\ &= \sum_{c, k}\left(\frac{\rho_{c, k}\Delta z_{c, k - 1/2}}{4}w_{c, k + 1/2}^2 + \frac{\rho_{c, k}\Delta z_{c, k + 1/2}}{4}w_{c, k + 1/2}^2\right)A_{c}\tag{34.14}\label{eq:h_kin_w} \end{align} \]
An alternative formulation of kinetic energy is
\[ \begin{align} H_\text{kin}^\star\left(\mathbf{v}\right) &\coloneqq \sum_{c, k}\frac{1}{2}\left(\rho\mathbf{v}\cdot\mathbf{v}\right)_{c, k}V_{c, k},\tag{34.15}\label{eq:h_kin_mod}\\ H_\text{kin}^{(u, \star)}\left(u\right) &\coloneqq \sum_{c, k}\sum_{e\in c}\frac{l_ed_e\Delta z_{e, k}}{4}\newtilde{\rho_{c, k}}^{(e)}u_{e, k}^2,\tag{34.16}\label{eq:h_kin_u_mod}\\ H_\text{kin}^{(w, \star)}\left(w\right) &\coloneqq \sum_{c, k}\left(\frac{\Delta z_{c, k - 1/2}}{4}\newtilde{\rho_{c, k}}^{(k - 1/2)}w_{c, k - 1/2}^2 + \frac{\Delta z_{c, k + 1/2}}{4}\newtilde{\rho_{c, k}}^{(k + 1/2)}w_{c, k + 1/2}^2\right)A_{c}.\tag{34.17}\label{eq:h_kin_w_mod} \end{align} \]
Here $\newtilde{\rho _{c, k}}^{(e)}$ is an averaging operator on the edges and $\newtilde{\rho _{c, k}}^{(k + 1/2)}$ is an averaging operator on the points where $w$ is placed. These operators must satisfy
\[ \begin{align} \sum_{e, k}\newtilde{\rho_{c, k}}^{(e)}V_{e, k} &= \sum_{c, k}\rho_{c, k}V_{c, k},\\ \sum_{c, k}\newtilde{\rho_{c, k}}^{(k + 1/2)}V_{c, k + 1/2} &= \sum_{c, k}\rho_{c, k}V_{c, k} \end{align} \]
in order not to jeopardize mass conservation.
Eq. (A.124) is discretized in the form
\[ \begin{align} \lim_{\epsilon\to 0}\frac{1}{\epsilon}\sum_{c, k}\left[f\left(\mathbf{u}_{c, k} + \epsilon\mathbf{h}_{c, k}\right) - f\left(\mathbf{u}_{c, k}\right)\right]V_{c, k} \hastobe \sum_{c, k}\left(\frac{\delta F}{\delta\mathbf{u}}\right)_{c, k}\cdot\mathbf{h}_{c, k}V_{c, k}. \end{align} \]
The functional derivative is the discrete function that fulfills this for all auxiliary functions $\mathbf{h}_{c, k}$. $\mathbf{h}$ is a tuple of state variables. It can also be summed over volumes other than $V_{c, k}$.
For the functional derivatives with respect to the thermodynamic quantities one obtains
\[ \begin{align} \left(\frac{\delta H_\text{kin}}{\delta\rho}\right)_{c, k} = \frac{1}{2}\left(\mathbf{v}\cdot\mathbf{v}\right)_{c, k}, & {} & \left(\frac{\delta H_{I'}}{\delta\rho}\right)_{c, k} = \phi_{c, k} + \frac{\partial\newtilde{I}_{c, k}}{\partial\rho}, & {} & \left(\frac{\delta H_{I'}}{\delta q_1}\right)_{c, k} = \frac{\partial\newtilde{I}_{c, k}}{\partial q_1}. \end{align} \]
The equations (34.16) - (34.17) can be rewritten in the form
\[ \begin{align} H_\text{kin}^{(u, \star)}\left(u\right) &\coloneqq \sum_{c, k}\sum_{e\in c}V_{e, k}\newtilde{\rho_{c, k}}^{(e)}\frac{u_{e, k}^2}{2},\\ H_\text{kin}^{(w, \star)}\left(w\right) &\coloneqq \sum_{c, k}\left(\frac{\Delta z_{c, k - 1/2}}{4}\newtilde{\rho_{c, k}}^{(k - 1/2)}w_{c, k - 1/2}^2 + \frac{\Delta z_{c, k + 1/2}}{4}\newtilde{\rho_{c, k}}^{(k + 1/2)}w_{c, k + 1/2}^2\right)A_{c}\nonumber\\ &= \sum_{c, k}\left(\frac{\Delta z_{c, k - 1/2}}{4}\newtilde{\rho_{c, k}}^{(k - 1/2)}w_{c, k - 1/2}^2 + \frac{\Delta z_{c, k + 1/2}}{4}\newtilde{\rho_{c, k}}^{(k + 1/2)}w_{c, k + 1/2}^2\right)A_{c}\nonumber\\ &= \sum_{c, k}\newtilde{\rho_{c, k}}^{(k + 1/2)}\frac{w_{c, k + 1/2}^2}{2}V_{c, k + 1/2} \end{align} \]
Here,
\[ \begin{align} V_{e, k} = \frac{l_ed_e\Delta z_{e, k}}{2} \end{align} \] and
\[ \begin{align} V_{c, k + 1/2} = A_c\Delta z_{c, k + 1/2} \end{align} \]
were substituted. This implies
\[ \begin{align} \left(\frac{\delta H_\text{kin}^\star}{\delta\rho}\right)_{c, k} = \frac{1}{2}\left(\mathbf{v}\cdot\mathbf{v}\right)_{c, k}, \end{align} \]
since the functional derivative can also be evaluated at the boundaries of the boxes and then interpolated back to the center points using a scalar product.
From Eq. (34.13) follows
\[ \begin{align} \frac{H_\text{kin}^{(u)}\left(u + \epsilon h\right) - H_\text{kin}^{(u)}\left(u\right)}{\epsilon} &= \sum_{c, k}\sum_{e\in c}\frac{\rho_{c, k}l_ed_e\Delta z_{e, k}}{2}u_{e, k}h_{e, k} + O\left(h^2\right) = \sum_{c, k}\sum_{e\in c}\rho_{c, k}u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right).\tag{34.27}\label{eq:shallow_func_deriv_deriv_0} \end{align} \]
The sum over cells can be rewritten as a sum over edges
\[ \begin{align} \frac{H_\text{kin}^{(u)}\left(u + \epsilon h\right) - H_\text{kin}^{(u)}\left(u\right)}{\epsilon} &= \sum_{e, k}\left(\rho_{c_1(e), k} + \rho_{c_2(e), k}\right)u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right)\nonumber\\ &= 2\sum_{e, k}\frac{\rho_{c_1(e), k} + \rho_{c_2(e), k}}{2}u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right) \end{align} \]
The prefactor two stands for the two spatial directions under consideration. The two cells bordering the edge $e$ were labeled $c_0$ and $c_1$. Thus one has
\[ \begin{align} \left(\frac{\delta H}{\delta\mathbf{v}}\right)_{e, k} = \frac{\rho_{c_1(e), k} + \rho_{c_2(e), k}}{2}u_{e, k}. \end{align} \]
Analogously it follows from Eq. (34.14)
\[ \begin{align} \frac{H_\text{kin}^{(w)}\left(w + \epsilon h\right) - H_\text{kin}^{(w)}\left(w\right)}{\epsilon} &= \sum_{c, k}\left(\frac{\rho_{c, k}\Delta z_{c, k - 1/2}}{2}w_{c, k - 1/2}h_{c, k - 1/2} + \frac{\rho_{c, k}\Delta z_{c, k + 1/2}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\right)A_{c}\nonumber\\ &= \sum_{c, k}\left(\frac{\rho_{c, k}}{2}w_{c, k + 1/2}h_{c, k + 1/2} + \frac{\rho_{c, k + 1}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\right)\Delta z_{c, k + 1/2}A_{c}.\tag{34.30}\label{eq:shallow_func_deriv_deriv_1} \end{align} \]
An index shift was made. With $V_{c, k + 1/2} = A_c\Delta z_{c, k + 1/2}$ one obtains
\[ \begin{align} \frac{H_\text{kin}^{(w)}\left(w + \epsilon h\right) - H_\text{kin}^{(w)}\left(w\right)}{\epsilon} &= \sum_{c, k}\frac{\rho_{c, k} + \rho_{c, k + 1}}{2}w_{c, k + 1/2}h_{c, k + 1/2}V_{c, k + 1/2}. \end{align} \]
It thus follows
\[ \begin{align} \left(\frac{\delta H}{\delta\mathbf{v}}\right)_{c, k + 1/2} = \frac{\rho_{c, k} + \rho_{c, k + 1}}{2}w_{c, k + 1/2}. \end{align} \]
Using the modified operator $H_\text{kin}^\star\left(\mathbf{v}\right)$ from Eq. (34.15), one obtains
\[ \begin{align} \left(\frac{\delta H^\star}{\delta\mathbf{v}}\right)_{e, k} &= \newtilde{\rho_{c, k}}^{(e)}u_{e, k},\tag{34.33}\label{eq:h_mod_func_deriv_0}\\ \left(\frac{\delta H^\star}{\delta\mathbf{v}}\right)_{c, k + 1/2} &= \newtilde{\rho_{c, k}}^{(k + 1/2)}w_{c, k + 1/2}.\tag{34.34}\label{eq:h_mod_func_deriv_1} \end{align} \]
Eq. (34.27) generalizes with
\[ \begin{align} V_{e, k} \coloneqq A_{e,k}\frac{d_{e, k}}{2} \end{align} \]
to the deep atmosphere as follows:
\[ \begin{align} \frac{H_\text{kin}^{(u)}\left(u + \epsilon h\right) - H_\text{kin}^{(u)}\left(u\right)}{\epsilon} &= \sum_{c, k}\sum_{e\in c}\frac{\rho_{c, k}d_eA_{e, k}}{2}u_{e, k}h_{e, k} + O\left(h^2\right) = \sum_{c, k}\sum_{e\in c}\rho_{c, k}u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right) \end{align} \]
Now one again writes the sum over cells as a sum over edges:
\[ \begin{align} \frac{H_\text{kin}^{(u)}\left(u + \epsilon h\right) - H_\text{kin}^{(u)}\left(u\right)}{\epsilon} &= \sum_{e, k}\left(\rho_{c_1(e), k} + \rho_{c_2(e), k}\right)u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right)\nonumber\\ &= 2\sum_{e, k}\frac{\rho_{c_1(e), k} + \rho_{c_2(e), k}}{2}u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right) \end{align} \]
Thus, in the deep atmosphere as well, one has
\[ \begin{align} \left(\frac{\delta H}{\delta\mathbf{v}}\right)_{e, k} = \frac{\rho_{c_1(e), k} + \rho_{c_2(e), k}}{2}u_{e, k}. \end{align} \]
One proceeds analogously with Eq. (34.27):
\[ \begin{align} \frac{H_\text{kin}^{(w)}\left(w + \epsilon h\right) - H_\text{kin}^{(w)}\left(w\right)}{\epsilon} &= \sum_{c, k}\left(\frac{\rho_{c, k}A_{c, k - 1/2}\Delta z_{c, k - 1/2}}{2V_{c, k}}w_{c, k - 1/2}h_{c, k - 1/2} + \frac{\rho_{c, k}A_{c, k + 1/2}\Delta z_{c, k + 1/2}}{2V_{c, k}}w_{c, k + 1/2}h_{c, k + 1/2}\right)V_{c, k}\nonumber\\ &= \sum_{c, k}\frac{\rho_{c, k}A_{c, k - 1/2}\Delta z_{c, k - 1/2}}{2}w_{c, k - 1/2}h_{c, k - 1/2} + \frac{\rho_{c, k}A_{c, k + 1/2}\Delta z_{c, k + 1/2}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\nonumber\\ &= \sum_{c, k}\frac{\rho_{c, k + 1}A_{c, k + 1/2}\Delta z_{c, k + 1/2}}{2}w_{c, k + 1/2}h_{c, k + 1/2} + \frac{\rho_{c, k}A_{c, k + 1/2}\Delta z_{c, k + 1/2}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\nonumber\\ &= \sum_{c, k}\left(\frac{\rho_{c, k + 1}}{2}w_{c, k + 1/2}h_{c, k + 1/2} + \frac{\rho_{c, k}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\right)\Delta z_{c, k + 1/2}A_{c, k + 1/2} \end{align} \]
With $V_{c, k + 1/2} = A_c\Delta z_{c, k + 1/2}$ one obtains
\[ \begin{align} \frac{H_\text{kin}^{(w)}\left(w + \epsilon h\right) - H_\text{kin}^{(w)}\left(w\right)}{\epsilon} &= \sum_{c, k}\frac{\rho_{c, k} + \rho_{c, k + 1}}{2}w_{c, k + 1/2}h_{c, k + 1/2}V_{c, k + 1/2}. \end{align} \]
It thus follows
\[ \begin{align} \left(\frac{\delta H}{\delta\mathbf{v}}\right)_{c, k + 1/2} = \frac{\rho_{c, k} + \rho_{c, k + 1}}{2}w_{c, k + 1/2}. \end{align} \]
The equations (34.33) - (34.34) also apply in the deep atmosphere.
In Sect. 26.7 the requirement $\mathbf{v}\cdot\nabla\psi \hastobe \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v}$ was taken as the starting point in order to derive a discretization of the scalar product; here the discretizations of „scalar field x vector field“, gradient and divergence were presupposed. In Sect. 33.4.2 the scalar product was generalized to terrain-following coordinates. From this, the vertical contravariant measure could be reconstructed in Sect. 33.4.3. This was used in Sect. 33.4.4 to derive the shape of the divergence in terrain-following coordinates. Now the gradient still has to be generalized to terrain-following coordinates. For this, one uses
\[ \begin{align} \mathbf{v}\cdot\nabla\psi \equiv \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v} \end{align} \]
as the equation defining the gradient.
From Eq. (33.69) follow
\[ \begin{align} \psi\nabla\cdot\mathbf{v} &= \left[\sum_{e\in c}\left(l_e\frac{\Delta z_{e, k}\psi_{c, k}u_e}{A_c\Delta z_{c, k}}\right) + \frac{\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}}{\Delta z_{c, k}}\right]\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4A_c\Delta z_{c, k}}\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1} - \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\tag{34.43}\label{eq:grad_c-grid_terrain_deriv_0} \end{align} \]
and
\[ \begin{align} \nabla\cdot\left(\psi\mathbf{v}\right) &= \left[\sum_{e\in c}\left(l_e\frac{\Delta z_{e, k}\left(\psi_{o(e), k} + \psi_{c, k}\right)u_e}{2A_c\Delta z_{c, k}}\right) + \frac{\left(\psi_{c, k - 1} + \psi_{c, k}\right)w_{c, k - 1/2} - \left(\psi_{c, k + 1} + \psi_{c, k}\right)w_{c, k + 1/2}}{2\Delta z_{c, k}}\right]\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4A_c\Delta z_{c, k}}\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}\frac{\psi_{c, k - 1} + \psi_{o(e), k - 1}}{2}u_{e, k - 1} - \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}\frac{\psi_{c, k + 1} + \psi_{o(e), k + 1}}{2}u_{e, k + 1}\right).\tag{34.44}\label{eq:grad_c-grid_terrain_deriv_1} \end{align} \]
Subtracting Eq. (34.43) from Eq. (34.44), one obtains
\[ \begin{align} &\mathbf{v}\cdot\nabla\psi = \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v}\nonumber\\ &= \left[\sum_{e\in c}\left(l_e\frac{\Delta z_{e, k}\left(\psi_{o(e), k} - \psi_{c, k}\right)u_e}{2A_c\Delta z_{c, k}}\right) + \frac{\left(\psi_{c, k - 1} - \psi_{c, k}\right)w_{c, k - 1/2} - \left(\psi_{c, k + 1} - \psi_{c, k}\right)w_{c, k + 1/2}}{2\Delta z_{c, k}}\right]\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4A_c\Delta z_{c, k}}\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}\frac{\psi_{c, k - 1} + \psi_{o(e), k - 1} - 2\psi_{c, k}}{2}\textcolor{red}{u_{e, k - 1}} - \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}\frac{\psi_{c, k + 1} + \psi_{o(e), k + 1} - 2\psi_{c, k}}{2}\textcolor{red}{u_{e, k + 1}}\right)\nonumber\\ &\stackrel{\text{Eq. }\href{ch-32-kinematics.html#eq:inner_c-grid_3d_shallow_terrain}{(33.47)}}{\hastobe} \sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Delta z_{e, k}\left(\nabla\psi\right)_eu_e + \frac{\Delta z_{c, k - 1/2}}{2\Delta z_{c, k}}\left(\nabla\psi\right)_{c, k - 1/2}w_{c, k - 1/2} - \frac{\Delta z_{c, k + 1/2}}{2\Delta z_{c, k}}\left(\nabla\psi\right)_{c, k + 1/2}w_{c, k + 1/2}. \end{align} \]
The terms marked in red appear only on the left-hand side, but not on the right. $\nabla\psi$ can therefore no longer be defined in terrain-following coordinates in such a way that this equation is satisfied. One could eliminate these terms by calculating the contravariant flux densities only via a one-sided vertical reconstruction, which would be a contradiction to the requirement of mass conservation, since in this case the flux through a surface would depend on which side the surface is viewed from. So the problem cannot be solved exactly.
Nevertheless, conservation of total energy can still be achieved. For this it is necessary that the identity
\[ \begin{align} \int_A\psi\nabla\cdot\mathbf{v}d^3r = -\int_A\mathbf{v}\cdot\nabla\psi d^3r \end{align} \]
carries over to the discretization. For this calculation, the orientations $\gamma_{e(c)}$ of the vectors must also be taken into account. One defines
\[ \begin{align} \gamma_{e(c)} \coloneqq \begin{cases} 1,\text{ if the vector at }e\text{ points outward with respect to }c\text{,}\\ -1,\text{ otherwise} \end{cases} \end{align} \]
Thus one obtains
\[ \begin{align} &\int_A\psi\nabla\cdot\mathbf{v}d^3r \to \sum_{c, k}\left(\psi\nabla\cdot\mathbf{v}\right)_{c, k}A_c\Delta z_{c, k}\nonumber\\ &= \sum_{c, k}\Bigg[\sum_{e\in c}\left(l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e\right) + A_c\left(\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}\right)\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1} - \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\Bigg]. \end{align} \]
An index shift yields
\[ \begin{align} &\sum_{c, k}\left(\psi\nabla\cdot\mathbf{v}\right)_{c, k}A_c\Delta z_{c, k} = \sum_{c, k}\Bigg[\sum_{e\in c}\left(l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e\right) + A_c\left(\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}\right)\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k + 1}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k - 1}u_{e, k}\right)\Bigg]\nonumber\\ &= \sum_{c, k}\Bigg[\sum_{e\in c}\left(l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e\right) + A_c\left(\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}\right)\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k + 1}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k - 1}u_{e, k} + \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k}u_{e, k}\right)\Bigg]\nonumber \end{align} \] \[ \begin{align} &= \sum_{c, k}\sum_{e\in c}l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e + \sum_{c, k + 1/2}A_c\left(\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}\right)\nonumber\\ &- \sum_{e, k}\sum_{c\in e}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k + 1}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k - 1}u_{e, k} + \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k}u_{e, k}\right)\nonumber\\ &= \sum_{e, k}\sum_{c\in e}l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\frac{\psi_{c, k} - \psi_{c, k + 1}}{\Delta z_{c, k + 1/2}}\nonumber\\ &- \sum_{e, k}\sum_{c\in e}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k}}{\Delta z_{c, k}}\Delta z_{c, k + 1/2}J_{e, k}\frac{\psi_{c, k + 1} - \psi_{c, k}}{\Delta z_{c, k + 1/2}}u_{e, k} + \frac{\Delta z_{e, k}}{\Delta z_{c, k}}\Delta z_{c, k - 1/2}J_{e, k}u_{e, k}\frac{\psi_{c, k} - \psi_{c, k - 1}}{\Delta z_{c, k - 1/2}}\right)\nonumber\\ &= \sum_{e, k}-l_ed_e\Delta z_{e, k}u_e\sum_{c\in e}-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e} - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\frac{\psi_{c, k} - \psi_{c, k + 1}}{\Delta z_{c, k + 1/2}}\nonumber \end{align} \] \[ \begin{align} &- \sum_{e, k}l_ed_e\Delta z_{e, k}u_{e, k}J_{e, k}\sum_{c\in e}\frac{1}{4}\left(\frac{\Delta z_{c, k + 1/2}}{\Delta z_{c, k}}\frac{\psi_{c, k + 1} - \psi_{c, k}}{\Delta z_{c, k + 1/2}} + \frac{\Delta z_{c, k - 1/2}}{\Delta z_{c, k}}\frac{\psi_{c, k} - \psi_{c, k - 1}}{\Delta z_{c, k - 1/2}}\right)\nonumber\\ &= \sum_{e, k}-l_ed_e\Delta z_{e, k}u_e\sum_{c\in e}-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e} - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\left(\nabla_z\psi\right)_{c, k + 1/2}\nonumber\\ &+ \sum_{e, k}l_ed_e\Delta z_{e, k}u_{e, k}J_{e, k}\frac{1}{2}\newoverline{\left(\frac{\Delta z_{c, k - 1/2}}{\Delta z_{c, k}}\left(\nabla_z\psi\right)_{c, k - 1/2} + \frac{\Delta z_{c, k + 1/2}}{\Delta z_{c, k}}\left(\nabla_z\psi\right)_{c, k + 1/2}\right)}^{(e)}\nonumber \end{align} \] \[ \begin{align} &= \sum_{e, k}-l_ed_e\Delta z_{e, k}u_e\sum_{c\in e}-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e} - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\left(\nabla_z\psi\right)_{c, k + 1/2}\nonumber\\ &+ \sum_{e, k}l_ed_e\Delta z_{e, k}u_{e, k}J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}\nonumber\\ &= \sum_{e, k}-l_ed_e\Delta z_{e, k}u_e\left[\sum_{c\in e}\left(-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e}\right) - J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}\right] - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\left(\nabla_z\psi\right)_{c, k + 1/2}\nonumber\\ &= -\left\{2\sum_{e, k}V_{e, k}u_e\left[\sum_{c\in e}\left(-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e}\right) - J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}\right] + \sum_{c, k + 1/2}V_{c, k + 1/2}w_{c, k + 1/2}\left(\nabla_z\psi\right)_{c, k + 1/2}\right\}. \end{align} \]
The factor 2 stands for the two spatial directions in the horizontal. This results in the discretization
\[ \begin{align} \left(\nabla\psi\right)_{e, k} = \frac{1}{d_e}\sum_{c\in e}\left(-\psi_{c, k}\gamma_{e(c)}\right) - J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)} \end{align} \]
of the gradient in terrain-following coordinates.
Even in the deep atmosphere, the discretized product rule of the Nabla operator is not locally valid in coordinates following the terrain. Here, according to Eq. (33.71)
\[ \begin{align} \int_A\psi\nabla\cdot\mathbf{v}d^3r&\to\sum_{c, k}\left(\psi\nabla\cdot\mathbf{v}\right)_{c, k}V_{c, k} = \sum_{c, k}\Bigg\lbrace\left[\sum_{e\in c}\left(A_{e, k}\psi_{c, k}u_{e, k}\right) + A_{c, k - 1/2}\psi_{c, k}w_{c, k - 1/2} - A_{c, k + 1/2}\psi_{c, k}w_{c, k + 1/2}\right]\nonumber\\ &+ \sum_{e\in c}\Bigg[-A_{c, k - 1/2}\left(\frac{A_{e, k}d_{e, k}}{4V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{A_{e, k - 1}d_{e, k - 1}}{4V_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1}\right)\nonumber\\ &+ A_{c, k + 1/2}\left(\frac{A_{e, k}d_{e, k}}{4V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{A_{e, k + 1}d_{e, k + 1}}{4V_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\Bigg]\Bigg\rbrace. \end{align} \]
Except for the correction terms for the coordinates following the terrain, the product rule of the Nabla operator is also locally valid here; the derivation in Sect. 26.7.5 still applies. In addition, a global integral disappears over a divergence here too, the reason for this is the same as in the shallow atmosphere. One can therefore limit the partial integration to the terms of the terrain-following correction:
\[ \begin{align} &\sum_{c, k}\sum_{e\in c}\Bigg[-A_{c, k - 1/2}\left(\frac{A_{e, k}d_{e, k}}{4V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{A_{e, k - 1}d_{e, k - 1}}{4V_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1}\right)\nonumber\\ &+ A_{c, k + 1/2}\left(\frac{A_{e, k}d_{e, k}}{4V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{A_{e, k + 1}d_{e, k + 1}}{4V_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\Bigg]\nonumber\\ &= \sum_{c, k}\sum_{e\in c}\Bigg[-A_{c, k - 1/2}\left(\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{V_{e, k - 1}}{2V_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1}\right)\nonumber\\ &+ A_{c, k + 1/2}\left(\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{V_{e, k + 1}}{2V_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\Bigg]\nonumber\\ &= \sum_{c, k}\sum_{e\in c}\Bigg[-A_{c, k - 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} - A_{c, k - 1/2}\frac{V_{e, k - 1}}{2V_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1}\nonumber\\ &+ A_{c, k + 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + A_{c, k + 1/2}\frac{V_{e, k + 1}}{2V_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\Bigg]\nonumber \end{align} \] \[ \begin{align} &= \sum_{c, k}\sum_{e\in c}\Bigg(-A_{c, k - 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} - A_{c, k + 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k + 1}u_{e, k}\nonumber\\ &+ A_{c, k + 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + A_{c, k - 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k - 1}u_{e, k}\Bigg)\nonumber\\ &= \sum_{c, k}\sum_{e\in c}J_{e, k}u_{e, k}\left(-\frac{A_{c, k - 1/2}}{2V_{c, k}}\psi_{c, k} - \frac{A_{c, k + 1/2}}{2V_{c, k}}\psi_{c, k + 1} + \frac{A_{c, k + 1/2}}{2V_{c, k}}\psi_{c, k} + \frac{A_{c, k - 1/2}}{2V_{c, k}}\psi_{c, k - 1}\right)V_{e, k}\nonumber\\ &= \sum_{c, k}\sum_{e\in c}J_{e, k}u_{e, k}\left(\frac{A_{c, k - 1/2}}{2V_{c, k}}\Delta z_{c, k - 1/2}\frac{\psi_{c, k - 1} - \psi_{c, k}}{\Delta z_{c, k - 1/2}} + \frac{A_{c, k + 1/2}}{2V_{c, k}}\Delta z_{c, k + 1/2}\frac{\psi_{c, k} - \psi_{c, k + 1}}{\Delta z_{c, k + 1/2}}\right)V_{e, k}\nonumber \end{align} \] \[ \begin{align} &= \sum_{c, k}\sum_{e\in c}J_{e, k}u_{e, k}\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}V_{e, k} = -\sum_{c, k}\sum_{e\in c}-J_{e, k}u_{e, k}\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}V_{e, k}\nonumber\\ &= -\sum_{c, k}\sum_{e\in c}-J_{e, k}2u_{e, k}\frac{1}{2}\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}V_{e, k} = -\sum_{e, k}J_{e, k}u_{e, k}V_{e, k}\left(-2\sum_{c\in e}\frac{1}{2}\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}\right)\nonumber\\ &= -2\sum_{e, k}u_{e, k}V_{e, k}\left(-J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}\right) \end{align} \]
Therefore, in the deep atmosphere as well, the following holds:
\[ \begin{align} \left(\nabla\psi\right)_{e, k} = \frac{1}{d_{e, k}}\sum_{c\in e}\left(-\psi_{c, k}\gamma_{e(c)}\right) - J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}, \end{align} \]
where only this order of the averaging operators is permitted here.
\[ \begin{align} \sum_{e, k}\newoverline{\rho}_c^{(e)}v_e\frac{\partial v_{e, k}}{\partial t} &= \sum_{e, k}\newoverline{\rho}_c^{(e)}v_e\newoverline{\newoverline{\rho w\cdot\eta_t}^{(e)}}^{(k)} \end{align} \]
Since the Hamilton function does not explicitly depend on time, it is equal to the total energy. Substituting $F = H$ into Eq. (10.66), one obtains
\[ \begin{align} \frac{dH}{dt} = \lbrace H, H\rbrace = \lbrace H, Z_a, H\rbrace + \lbrace H, H\rbrace_\rho + \lbrace H, H\rbrace_{\newtilde{s}}. \end{align} \]
Furthermore, the following hold:
\[ \begin{align} \lbrace H, Z_a, H\rbrace \stackrel{\href{ch-09-application-of-the-hamilton-formalism-to.html#eq:poisson_bracket_atm_gen_deriv_5}{\text{Eq. (10.54)}}}{=} 0, & {} & \lbrace H, H\rbrace_\rho \stackrel{\href{ch-09-application-of-the-hamilton-formalism-to.html#eq:poisson_bracket_atm_gen_deriv_3}{\text{Eq. (10.50)}}}{=} 0, & {} & \lbrace H, H\rbrace_{\newtilde{s}} \stackrel{\href{ch-09-application-of-the-hamilton-formalism-to.html#eq:poisson_bracket_atm_gen_deriv_4}{\text{Eq. (10.51)}}}{=} 0. \end{align} \]
This implies conservation of energy
\[ \begin{align} H = \text{const.} \end{align} \]
In order to transfer energy conservation to spatial discretization, three requirements must be taken into account:
The helicity bracket $\lbrace H, Z_a, H\rbrace = \lbrace K, Z_a, H\rbrace$ is energetically inactive.
The following applies: $\lbrace H, H\rbrace_\rho = \lbrace K, H\rbrace_\rho + \lbrace I', H\rbrace_\rho = 0$.
We have $\lbrace H, H\rbrace_{\newtilde{s}} = \lbrace K, H\rbrace_{\newtilde{s}} + \lbrace I', H\rbrace_{\newtilde{s}} = 0$.
The first point corresponds to the fact that the generalized Coriolis force does no work, which was already demonstrated in Chap. 28. The second and third statements have the same computational prerequisite. To this end, one adds Eqs. (10.78) - (10.79), which leads to
\[ \begin{align} \int_A\underbrace{-\rho\mathbf{v}\cdot\nabla\phi}_{\text{A}}\underbrace{ - c_d^{(p)}\rho\theta\mathbf{v}\cdot\nabla\Pi}_{\text{B}}\underbrace{ - \phi\nabla\cdot\left(\rho\mathbf{v}\right)}_{\text{A}}\underbrace{ - c_d^{(p)}\Pi\nabla\cdot\left(\rho\theta\mathbf{v}\right)}_{\text{B}}d^3r = 0\tag{34.58}\label{eq:poisson_energy_deriv} \end{align} \]
The integrals over the terms A and B each cancel out. The following two facts were used:
The calculation rule $\nabla\cdot\left(\psi\mathbf{v}\right) = \mathbf{v}\cdot\nabla\psi + \psi\nabla\cdot\mathbf{v}$.
Global integrals over $\nabla\cdot\left(\psi\mathbf{v}\right)$ vanish with Gauss' theorem due to kinematic boundary conditions.
Both facts also hold in the discretization. The proofs of this have already been given:
The definition of the gradient in Sect. 34.2.2 ensures that $\nabla\cdot\left(\psi\mathbf{v}\right) = \mathbf{v}\cdot\nabla\psi + \psi\nabla\cdot\mathbf{v}$ also holds in the discretization.
The fact that global integrals over divergences vanish has already been proven in Sect. 26.5.
Another aspect should be pointed out here: if one discretizes Eq. (34.58), $\theta$ occurs only at the edges of the grid, not in the cell centers. One can therefore freely choose the value of $\theta$ at the edge, as long as one also does so in the pressure-gradient acceleration. This will be used in Sect. 34.4 to compute $\theta$ at the edge in such a way that the accuracy of the flux operator $-\nabla\cdot\left(\rho\theta\mathbf{v}\right)$ is increased. The $\theta$-conservation itself is not violated by this, since what flows out of one grid box still flows into the neighboring one.
Pressure-gradient and gravitational acceleration as well as the generalized Coriolis force are thereby covered. The gradient of kinetic energy, however, cancels out in the Poisson brackets because, although it transports kinetic energy, it does not create or destroy any energy. This will be retraced here once more. To see this, one writes the momentum equation in the form
\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} = -\nabla k + \mathbf{l}, \end{align} \]
where pressure gradient and gravitational acceleration were combined in the term $\mathbf{l}$. The Coriolis force can be ignored here. Scalar multiplication by $\mathbf{v}$ yields
\[ \begin{align} \mathbf{v}\cdot\frac{\partial\mathbf{v}}{\partial t} = -\mathbf{v}\cdot\nabla k + \mathbf{v}\cdot\mathbf{l}\nonumber & {} & \Leftrightarrow\frac{\partial k}{\partial t} = -\mathbf{v}\cdot\nabla k + \mathbf{v}\cdot\mathbf{l}. \end{align} \]
Multiplying this by the mass density $\rho$ gives
\[ \begin{align} \rho\frac{\partial k}{\partial t} &= -\rho\mathbf{v}\cdot\nabla k + \rho\mathbf{v}\cdot\mathbf{l}.\tag{34.60}\label{eq:no-split-explicit_deriv_0} \end{align} \]
The continuity equation reads
\[ \begin{align} \frac{\partial\rho}{\partial t} = -\nabla\cdot\left(\rho\mathbf{v}\right). \end{align} \]
Multiplying this by the specific kinetic energy $k$ gives
\[ \begin{align} k\frac{\partial\rho}{\partial t} = -k\nabla\cdot\left(\rho\mathbf{v}\right).\tag{34.62}\label{eq:no-split-explicit_deriv_1} \end{align} \]
Adding Eqs. (34.60) and (34.62), one obtains
\[ \begin{align} \frac{\partial\left(\rho k\right)}{\partial t} = -\nabla\cdot\left(\rho k\mathbf{v}\right) + \mathbf{v}\cdot\mathbf{l}.\tag{34.63}\label{eq:no-split-explicit_deriv_2} \end{align} \]
This means that only the global integral over $\mathbf{v}\cdot\mathbf{l}$ contributes to the evolution of the total kinetic energy of the atmosphere. The derivation only includes calculation rules for $\frac{\partial}{\partial t}$ and the rule $\nabla\cdot\left(\psi\mathbf{v}\right) = \mathbf{v}\cdot\nabla\psi + \psi\nabla\cdot\mathbf{v}$, which is already considered to be fulfilled. The acceleration $-\nabla k$ does not influence the energetic integrals even in the discretization.
Thus, the spatial discretization developed in this chapter is energetically self-consistent.
A hydrostatic background state as described in Sect. 13.7.5 can be used without destroying energy conservation. It must be ensured that
\[ \begin{align} \mathbf{0} = -c^{(p)}\newoverline{\theta}\nabla\newoverline{\Pi} - \nabla\phi \end{align} \]
also holds discretely. The corresponding terms then cancel each other out in the discretization.
However, this is only the case in the absence of orography. If mountains are present, the horizontal components of the gradients do not cancel exactly, giving rise to a small error in energy conservation. In the interplay of spatial and temporal discretization, however, the advantages of the hydrostatic background state outweigh the disadvantages. The most important advantage here is that the horizontal orographic correction components of $-c^{(p)}\theta\nabla\Pi$ and $-\nabla\phi$ are two quantities that almost completely cancel each other out, i.e.
\[ \begin{align} -c^{(p)}\newoverline{\theta_{c, k}}^{(e)}J_{e, k}\newoverline{\left(\nabla_z\Pi\right)_{c, k + 1/2}}^{(k)} - J_{e, k}\newoverline{\left(\nabla_z\phi\right)_{c, k + 1/2}}^{(k)} \approx 0. \end{align} \]
The remaining difference is largely a discretization error. This becomes smaller if one uses a hydrostatic background state, since then only
\[ \begin{align} -c^{(p)}\newoverline{\theta_{c, k}}^{(e)}J_{e, k}\newoverline{\left(\nabla_z\Pi'\right)_{c, k + 1/2}}^{(k)} - c^{(p)}\newoverline{\theta'_{c, k}}^{(e)}J_{e, k}\newoverline{\left(\nabla_z\newoverline{\Pi}\right)_{c, k + 1/2}}^{(k)} \end{align} \]
must be calculated.
The kinematic boundary condition $\mathbf{v}\cdot\mathbf{n} = 0$ is used as the upper boundary condition in the context of reversible dynamics. In the model, this corresponds to
\[ \begin{align} w_{c, 1/2} = 0 \end{align} \]
for all cells $c$.
The lower boundary condition in the context of reversible dynamics is the no-slip boundary condition $\mathbf{v} = \mathbf{0}$. Since only vertical vector components are located there, this corresponds to
\[ \begin{align} w_{c, N_L + \frac{1}{2}} = 0 \end{align} \]
for all cells $c$.
Within the framework of the formalism developed so far, there are three possible approaches to increasing the convergence order:
Modification of the kinetic energy functional including interpolation of the density $\newtilde{\rho}^{(e)}_{c, k}$ onto the edge.
Interpolation of $\theta$ onto the edge (see Sect. 34.4.1).
Increasing the number of edges taken into account when calculating the vorticity flux term.
These remaining freedoms are not sufficient to increase the convergence order of the overall model.
The accuracy of the discretization can be increased by interpolating $\theta$ onto the edge not via the arithmetic mean
\[ \begin{align} \newtilde{\theta}^{(e)} = \frac{\theta_{c_1(e)} + \theta_{c_2(e)}}{2}, \end{align} \]
which is of second order, but in a more complex way. The vertical index $k$ is omitted here because two-dimensionality is assumed in this section, and $c_1(e)$, $c_2(e)$ are, as usual, the two cells that border the edge $e$.
First, this is considered one-dimensionally. To do this, one assumes the true spatial dependence $\theta = \theta\left(x\right)$ and assumes that the values of $\theta$ are known at the scalar grid points $i$, which are positioned at locations $x_i = i\Delta x$, i.e.
\[ \begin{align} \theta_i = \theta\left(x_i\right). \end{align} \]
A suitable approximation for $\newtilde{\theta}^{(e)}$ is now sought, i.e. $\theta$ at position $x = \frac{\Delta x}{2}$. To do this, one first replaces the spatial dependence $\theta\left(x\right)$ by $\mu^{(n)}\left(x\right)$, where $n$ stands for the convergence order.
For $n = 2$ one makes the ansatz
\[ \begin{align} \mu^{(2)}\left(0\right) &= \theta\left(0\right),\tag{34.71}\label{eq:skam_gass_deriv_1}\\ \mu^{(2)}\left(\Delta x\right) &= \theta\left(\Delta x\right).\tag{34.72}\label{eq:skam_gass_deriv_2} \end{align} \]
Since one has two linear equations for $\mu^{(2)}$, one makes an ansatz with two coefficients, i.e.
\[ \begin{align} \mu^{(2)}\left(x\right) &= a + bx. \end{align} \]
Substituting this into Eqs. (34.71) - (34.72), one obtains
\[ \begin{align} a &= \theta\left(0\right),\\ a + b\Delta x &= \theta\left(\Delta x\right). \end{align} \]
This leads to
\[ \begin{align} a &= \theta_{c_1(e)},\\ b &= \frac{\theta\left(\Delta x\right) - \theta\left(0\right)}{\Delta x} = \frac{\theta_{c_2(e)} - \theta_{c_1(e)}}{\Delta x}. \end{align} \]
In order to calculate $\newtilde{\theta}^{(e,2)}$ from this, where the superscript two stands for the order, an integration must be carried out over the time interval $\left[0, \Delta t\right]$, where $\Delta t$ is the length of a time step. One calculates:
\[ \begin{align} \newtilde{\theta}^{(e,2)} &= a + b\frac{\Delta x}{2} = \theta_{c_1(e)} + \frac{\theta\left(\Delta x\right) - \theta\left(0\right)}{2} = \frac{\theta\left(0\right) + \theta\left(\Delta x\right)}{2} \end{align} \]
To obtain a third-order approximation ($n = 3$), one includes the second derivative of $\theta$ in the grid cell $c_1\left(e\right)$ in the reconstruction of $\theta$, i.e
\[ \begin{align} \mu^{(3)}\left(x\right) &= a + bx + cx^2. \end{align} \]
From this follow
\[ \begin{align} f_0 &= a,\\ f_1 &= a + b\Delta x + c\Delta x^2,\\ f_{-1} &= a - b\Delta x + c\Delta x^2, \end{align} \]
which can be rearranged to
\[ \begin{align} a &= f_0,\\ f_{-1} + f_1 &= 2a + 2c\Delta x^2 = 2\left(f_0 + c\Delta x^2\right) \Rightarrow \frac{f_{-1} + f_1}{2} - f_0 = c\Delta x^2 \Rightarrow c = \frac{f_{-1} - 2f_0 + f_1}{2\Delta x^2},\\ f_1 - f_{-1} &= 2b\Delta x \Rightarrow b = \frac{f_1 - f_{-1}}{2\Delta x} \end{align} \]
From this it follows
\[ \begin{align} \mu^{(3)}\left(x\right) &= f_0 + \frac{f_1 - f_{-1}}{2\Delta x}x + \frac{f_{-1} - 2f_0 + f_1}{2\Delta x^2}x^2. \end{align} \]
The value of this function on the edge ($x = \frac{\Delta x}{2}$) is therefore
\[ \begin{align} \mu^{(3)}\left(\frac{\Delta x}{2}\right) &= f_0 + \frac{f_1 - f_{-1}}{4} + \frac{f_{-1} - 2f_0 + f_1}{8} = \frac{8f_0 + 2f_1 - 2f_{-1} + f_{-1} - 2f_0 + f_1}{8}\nonumber\\ &= \frac{6f_0 + 3f_1 - f_{-1}}{8} = \frac{4f_0 + 4f_1 + 2f_0 - f_1 - f_{-1}}{8} = \frac{f_0 + f_1}{2} + \frac{2f_0 - f_1 - f_{-1}}{8}\nonumber\\ &= \frac{f_0 + f_1}{2} - \frac{f_1 - 2f_0 + f_{-1}}{8} = \frac{f_0 + f_1}{2} - \Delta x^2\frac{f_1 - 2f_0 + f_{-1}}{8\Delta x^2}\nonumber\\ &= \frac{f_0 + f_1}{2} - \frac{\Delta x^2}{8}\delta_x^2f \end{align} \]
with
\[ \begin{align} \delta_x^2f \coloneqq \frac{f_1 - 2f_0 + f_{-1}}{\Delta x^2} \end{align} \]
as an approximation for the second derivative of $\theta$ in the $x-$direction at the center of the cell $c_1\left(e\right)$.
In two dimensions, the one-dimensional formulas are simply applied in every spatial direction. This is not a problem on quadrangular grids where all vectors are parallel to a global coordinate axis. On grids that do not meet this criterion, this is more complicated, which is what the next section will cover.
On the hexagonal grid, the vectors are not aligned along coordinate lines, which complicates the calculation of the second derivatives at the cell centers. Therefore, one makes the ansatz for $\theta$
\[ \begin{align} g = g\left(x,y\right) = c_1 + c_2x + c_3y + c_4x^2 + c_5xy + c_6y^2. \end{align} \]
$x$ and $y$ are calculated on the tangent plane to the sphere in the center of the cell under consideration. The $x-$axis points in the direction of the edge on which $\newtilde{\theta}^{(e)}$ is to be calculated. Now define the coefficient vector
\[ \begin{align} \mathbf{f} \coloneqq \left(\begin{array}{c} c_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6 \end{array}\right). \end{align} \]
Only $c_4$ is relevant for the second derivative in the $x-$direction at the center of the cell, one has
\[ \begin{align} \frac{\partial^2 g}{\partial x^2}\left(0,0\right) = 2c_4.\tag{34.91}\label{eq:adv_hex_deriv} \end{align} \]
The goal now is to calculate the vector $\mathbf{f}$. To do this, one first defines the vector
\[ \begin{align} \mathbf{s} \coloneqq \left(\begin{array}{c} \theta\left(x_1,y_1\right)\\ \vdots\\ \theta\left(x_m,y_m\right) \end{array}\right) \end{align} \]
of the actual values of $\theta$ in the cell centers. Here, $m \geq 6$ is the number of points used to compute $\mathbf{f}$; these are the cell in whose center the second derivative is to be computed and all neighboring cells. One now defines an error function
\[ \begin{align} \psi\left(\mathbf{f}\right) \coloneqq \sum_i^m\left[g\left(x_i,y_i\right) - s_i\right]^2. \end{align} \]
Now one defines the polynomial matrix
\[ \begin{align} P \coloneqq \left(\begin{array}{cccccc} 1 & x_1 & y_1 & x_1^2 & x_1y_1 & y_1^2\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & x_m & y_m & x_m^2 & x_my_m & y_m^2 \end{array}\right). \end{align} \]
This allows $\psi$ to be written in the form
\[ \begin{align} \psi\left(\mathbf{f}\right) &= \left(P\mathbf{f} - \mathbf{s}\right)^2 = \left(P\mathbf{f}\right)\cdot\left(P\mathbf{f}\right) - 2\left(P\mathbf{f}\right)\cdot\mathbf{s} + \mathbf{s}^2 = \sum_{i=1}^m\left(P\mathbf{f}\right)_i^2 - 2\sum_{i=1}^m\left(P\mathbf{f}\right)_is_i + \sum_{i=1}^ms_i^2\nonumber\\ &= \sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}f_j\right)^2 - 2\sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}f_j\right)s_i + \sum_{i=1}^ms_i^2 \end{align} \]
From this it follows
\[ \begin{align} \frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= 2\sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}\delta_{j,k}\right)^2\left(\sum_{j=1}^nP_{i, j}f_j\right) - 2\sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}\delta_{j,k}\right)s_i = 0\nonumber\\ \Leftrightarrow\frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= \sum_{i=1}^m\left(\sum_{j=1}^nP_{i,j}\delta_{j,k}\right)\left(\sum_{j=1}^nP_{i, j}f_j\right) - \sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}\delta_{j,k}\right)s_i = 0\nonumber\\ \Leftrightarrow\frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= \sum_{i=1}^mP_{i,k}\left(\sum_{j=1}^nP_{i, j}f_j\right) - \sum_{i=1}^mP_{i,k}s_i = \sum_{i=1}^m\sum_{j=1}^nP_{i,k}P_{i, j}f_j - \sum_{i=1}^mP_{i, k}s_i = 0\nonumber\\ \Leftrightarrow\frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= \sum_{j=1}^n\sum_{i=1}^mP^T_{k,i}P_{i,j}f_j - \sum_{i=1}^mP^T_{k,i}s_i = \sum_{j=1}^n\left(P^TP\right)_{k,j}f_j - \sum_{i=1}^mP^T_{k,i}s_i = 0\nonumber\\ \Leftrightarrow\frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= \left[\left(P^TP\right)\cdot\mathbf{f}\right]_k - \left(P^T\cdot\mathbf{s}\right)_k = 0. \end{align} \]
Thus the equivalence holds:
\[ \begin{align} \nabla_\mathbf{k}\psi = 0 \Leftrightarrow \left(P^TP\right)\cdot\mathbf{f} = P^T\cdot\mathbf{s}, \end{align} \]
which one can transform into
\[ \begin{align} \mathbf{f} = \left[\left(P^TP\right)^{-1}P^T\right]\cdot\mathbf{s} = B\cdot\mathbf{s} \end{align} \]
with
\[ \begin{align} B \coloneqq \left[\left(P^TP\right)^{-1}P^T\right] \end{align} \]
From this, the second derivatives in the cell centers in the x-direction can be calculated using Eq. (34.91).
There is the following conflict between explicit and implicit methods:
Explicit methods are easy to implement and fast. However, the time step resulting from the CFL criterion is usually significantly smaller than the physically necessary one. This significantly reduces the performance of explicit methods. Explicit multi-step procedures also require filters, which can destroy mimetic properties.
Implicit methods are more complex to implement. Their speed depends on the complexity of the problem. They are more stable than explicit methods and allow a large time step, which contributes positively to their performance.
It follows that a fast implicit method would be ideal for a dynamic core. However, since the equations to be solved are nonlinear, a three-dimensional implicit method is only possible via an iterative approach. In this section, the compromise chosen is the HEVI concept, which stands for horizontally explicit, vertically implicit.
A significant limitation with regard to the time step and thus also with regard to the efficiency of a model is the CFL criterion associated with sound waves. Therefore, the system of equations is often divided into so-called fast modes or divergent modes on the one hand and the slow modes, advective modes or non-divergent modes on the other hand. The fast modes include all terms associated with sound waves, while the slow modes include the remaining terms. With split-explicit techniques, the fast modes are now integrated with a time step that is approximately three times smaller. In the case of the equation system Equations (34.4) - (34.7), only the advection of the momentum for the larger time step would come into question. The advections of density and entropy, on the other hand, should be counted among the fast modes in order not to have to break down the divergences of the flux densities into advection and compression components. However, if one places high demands on the energetic self-consistency of a time step method, splitting off momentum advection does not make sense.
In order to see this, one must first realize that the advection of kinetic energy is connected to the continuity equation and is therefore one of the fast modes. This has already been seen in Sect. 34.2.4. Therefore, both terms or equations must be integrated with the same time step. In order to preserve the discretized analogue of the Lamb transformation $-\left(\mathbf{v}\cdot\nabla\right)\mathbf{v} = -\nabla k - \left(\nabla\times\mathbf{v}\right)\times\mathbf{v}$, the generalized Coriolis term and thus the complete system of equations must then also be integrated with the smaller time step.
Due to the disadvantages discussed here, no split-explicit option is built into the dynamic core being developed.
Vertically propagating waves are treated partially implicitly for reasons of stability. For this, one first writes
\[ \begin{align} \theta'' &\stackrel{\theta = \newtilde{\theta}/\rho}{=} \frac{\partial\theta}{\partial\rho}\rho' + \frac{\partial\theta}{\partial\newtilde{\theta}}\newtilde{\theta}' = -\frac{\newtilde{\theta}}{\rho^2}\rho' + \frac{1}{\rho}\newtilde{\theta}',\\ \Pi'' &\stackrel{\href{ch-08-first-law-in-the-atmosphere.html#eq:exner_pressure_diag}{\text{Eq. (9.71)}}}{=} \frac{R_s}{c^{(V)}\newtilde{\theta}}\left(\frac{R_s\newtilde{\theta}}{p_0}\right)^{R_s/c^{(V)}}\newtilde{\theta}'. \end{align} \]
The two primes stand for linear expansions relative to the previous time step and not for differences from the hydrostatic background state, which are denoted by $\theta',\Pi'$. The derivatives are abbreviated by
\[ \begin{align} \alpha \coloneqq -\frac{\newtilde{\theta}}{\rho^2},& \beta \coloneqq \frac{1}{\rho},\\ \gamma \coloneqq \frac{R_s}{c^{(V)}\newtilde{\theta}}\Pi \end{align} \]
These are calculated in the cell centers, using the old time step (time step $n$) for the predictor step and the arithmetic mean of the steps $n$ and $n + 1$ for the corrector step. The vertical momentum equation is written down in the form
\[ \begin{align} w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left[\left(\newoverline{\newoverline{\theta}_k}^{(k + 1/2)} + \newoverline{\theta_k^{'(n + 1)^\star}}^{(k + 1/2)}\right)\frac{\Pi^{'(n + 1)}_{k} - \Pi^{'(n + 1)}_{k + 1}}{z_{k} - z_{k + 1}} + \newoverline{\theta_k^{'(n + 1)}}^{(k + 1/2)}\frac{\newoverline{\Pi}_{k} - \newoverline{\Pi}_{k + 1}}{z_{k} - z_{k + 1}}\right]\nonumber\\ &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(n + 1)}_{k} - \Pi^{'(n + 1)}_{k + 1}}{z_{k} - z_{k + 1}} + \newoverline{\theta_k^{'(n + 1)}}^{(k + 1/2)}\frac{\newoverline{\Pi}_{k} - \newoverline{\Pi}_{k + 1}}{z_{k} - z_{k + 1}}\right)\nonumber\\ &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(n + 1)}_{k} - \Pi^{'(n + 1)}_{k + 1}}{z_{k} - z_{k + 1}} + \newoverline{\theta_k^{'(n + 1)}}^{(k + 1/2)}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\right).\tag{34.104}\label{eq:wave_vertical_deriv_0} \end{align} \]
Here $\epsilon$ is the implicit weight, for which $\frac{c^{(v)}}{c^{(p)}}$ is usually used based on the findings in Sect. 25.5.1.3. $\left(n + 1\right)^\star$ is a preliminary value for the new time step. For the implicit quantities, one now writes
\[ \begin{align} \theta^{'(n + 1)}_k &= \theta^{'(\text{expl.})}_k + \theta''_k = \theta^{'(\text{expl.})}_k + \alpha_k\rho'_k + \beta_k\newtilde{\theta}'_k,\\ \Pi^{'(n + 1)}_k &= \Pi^{'(\text{expl.})}_k + \Pi''_k = \Pi^{'(\text{expl.})}_k + \gamma_k\newtilde{\theta}'_{k}. \end{align} \]
The primed quantities are the deviations from the explicit components, i.e.
\[ \begin{align} \psi'_k \coloneqq \psi_k^{(n + 1)} - \psi_k^{(\text{expl.})} \end{align} \]
for $\psi = \rho, \newtilde{\theta}.$ For these, the following hold:
\[ \begin{align} \rho_k' = -\Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{V_k},\\ \newtilde{\theta}_k' = -\Delta t\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}}{V_k}. \end{align} \]
Here,
\[ \begin{align} M_{k + 1/2} \coloneqq A_{k + 1/2}\frac{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}w_{k + 1/2}^{(n + 1)} + \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}w_{k + 1/2}^{(n)}}{2}\tag{34.110}\label{eq:def_ml} \end{align} \]
for $1 \leq k \leq N_L - 1$ defines the vertical covariant mass flux density multiplied by the boundary surfaces $A_{k + 1/2}$ of the grid boxes; this is the vector of unknowns. Substituting this into Eq. (34.104), one obtains
\[ \begin{align} w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \Pi''_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \Pi''_{k + 1}}{z_{k} - z_{k + 1}} + \newoverline{\theta_k^{'(n + 1)}}^{(k + 1/2)}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\right)\nonumber\\ \Leftrightarrow w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \Pi''_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \Pi''_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\theta_{k}^{'(n + 1)}+\theta_{k + 1}^{'(n + 1)}}{2}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\right)\nonumber\\ \Leftrightarrow w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \gamma_{k}\newtilde{\theta}'_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\theta_{k}^{'(n + 1)}+\theta_{k + 1}^{'(n + 1)}}{2}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\right)\nonumber\\ \Leftrightarrow w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \gamma_{k}\newtilde{\theta}'_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\theta_{k}^{'(\text{expl.})}+\theta''_{k}+\theta_{k + 1}^{'(\text{expl.})}+\theta''_{k + 1}}{2}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\right)\\ \Leftrightarrow w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\cdot\Bigg(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \gamma_{k}\newtilde{\theta}'_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}}\nonumber\\ &+ \frac{\theta_{k}^{'(\text{expl.})}+\alpha_{k}\rho'_{k}+\beta_{k}\newtilde{\theta}'_{k}+\theta_{k + 1}^{'(\text{expl.})}+\alpha_{k + 1}\rho'_{k + 1}+\beta_{k + 1}\newtilde{\theta}'_{k + 1}}{2}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\Bigg). \end{align} \]
The primed quantities contain the vector of unknowns, so one moves them to the left-hand side:
\[ \begin{align} &\epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\gamma_{k}\newtilde{\theta}'_{k} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\alpha_{k}\rho'_{k}+\beta_{k}\newtilde{\theta}'_{k}+\alpha_{k + 1}\rho'_{k + 1}+\beta_{k + 1}\newtilde{\theta}'_{k + 1}}{2}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\right)\nonumber\\ &= w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)} - \epsilon\Delta tc^{(p)}\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{z_{k} - z_{k + 1}} - \epsilon\Delta tc^{(p)}\frac{\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}}{2}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\nonumber\\ &\Leftrightarrow\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\gamma_{k}\newtilde{\theta}'_{k} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\alpha_{k}\rho'_{k}+\beta_{k}\newtilde{\theta}'_{k}+\alpha_{k + 1}\rho'_{k + 1}+\beta_{k + 1}\newtilde{\theta}'_{k + 1}}{2}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\nonumber\\ &= \frac{w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)}}{\epsilon\Delta tc^{(p)}} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{z_{k} - z_{k + 1}} - \frac{\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}}{2}\frac{d\newoverline{\Pi}}{dz}\vert_{k+1/2}\nonumber\\ &\Leftrightarrow\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\left(\gamma_{k}\newtilde{\theta}'_{k} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}\right) + \frac{d\newoverline{\Pi}}{2dz}\vert_{k+1/2}\left(z_{k} - z_{k + 1}\right)\left(\alpha_{k}\rho'_{k}+\beta_{k}\newtilde{\theta}'_{k}+\alpha_{k + 1}\rho'_{k + 1}+\beta_{k + 1}\newtilde{\theta}'_{k + 1}\right)\nonumber\\ &= \frac{\left(w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)}\right)\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta tc^{(p)}} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\left(\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}\right) - \left(z_{k} - z_{k + 1}\right)\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vert_{k+1/2}\nonumber\\ &\Leftrightarrow\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg(-\gamma_{k}\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}}{V_{k}}\nonumber\\ &+ \gamma_{k + 1}\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}}{V_{k + 1}}\Bigg)\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vert_{k+1/2}\left(z_{k} - z_{k + 1}\right)\Bigg(-\alpha_{k}\frac{M_{k - 1/2} - M_{k + 1/2}}{V_{k}}-\beta_{k}\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}}{V_{k}}\nonumber\\ &-\alpha_{k + 1}\frac{M_{k + 1/2} - M_{k + 3/2}}{V_{k + 1}}-\beta_{k + 1}\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}}{V_{k + 1}}\Bigg)\nonumber\\ &= \frac{\left(w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)}\right)\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} - \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vert_{k+1/2}. \end{align} \]
One now makes the substitutions
\[ \begin{align} \alpha_k \to \frac{\alpha_k}{V_k}, \beta_k \to \frac{\beta_k}{V_k}, \gamma_k \to \frac{\gamma_k}{V_k} \end{align} \]
It follows
\[ \begin{align} &\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg[-\gamma_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &+ \gamma_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vert_{k+1/2}\left(z_{k} - z_{k + 1}\right)\Bigg[-\alpha_{k}\left(M_{k - 1/2} - M_{k + 1/2}\right)-\beta_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &-\alpha_{k + 1}\left(M_{k + 1/2} - M_{k + 3/2}\right)-\beta_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &= \frac{\left(w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)}\right)\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} - \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vert_{k+1/2}\nonumber\\ &\Leftrightarrow\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg[\gamma_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &- \gamma_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vert_{k+1/2}\left(z_{k} - z_{k + 1}\right)\Bigg[\alpha_{k}\left(M_{k - 1/2} - M_{k + 1/2}\right)+\beta_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &+\alpha_{k + 1}\left(M_{k + 1/2} - M_{k + 3/2}\right)+\beta_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg] - \frac{w_{k + 1/2}^{(n + 1)}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}}\nonumber\\ &= -\frac{w_{k + 1/2}^{(\text{expl.})}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} + \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vert_{k + 1/2}.\tag{34.115}\label{eq:vert_solver_pre_1} \end{align} \]
$w_{k + 1/2}^{(n + 1)}$ also depends on $M_{k + 1/2}$. From Eq. (34.110) follows
\[ \begin{align} w^{(n + 1)}_{k + 1/2} &= \frac{\frac{2M_{k + 1/2}}{A_{k + 1/2}} - \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}w^{(n)}_{k + 1/2}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}.\tag{34.116}\label{eq:w_from_ml} \end{align} \]
Putting this into Eq. (34.115), one obtains
\[ \begin{align} &\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg[\gamma_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &- \gamma_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vert_{k+1/2}\left(z_{k} - z_{k + 1}\right)\Bigg[\alpha_{k}\left(M_{k - 1/2} - M_{k + 1/2}\right)+\beta_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &+\alpha_{k + 1}\left(M_{k + 1/2} - M_{k + 3/2}\right)+\beta_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &- \frac{\frac{2M_{k + 1/2}}{A_{k + 1/2}} - \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}w^{(n)}_{k + 1/2}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}\frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}}\nonumber\\ &= -\frac{w_{k + 1/2}^{(\text{expl.})}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} + \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vert_{k + 1/2}. \end{align} \]
For the density $\newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}$, one has
\[ \begin{align} \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)} &= \frac{1}{2}\rho_{k}^{(n + 1)} + \frac{1}{2}\rho_{k + 1}^{(n + 1)}\nonumber\\ &= \frac{1}{2}\rho_{k}^{(\text{expl.})} - \Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{2V_{k}} + \frac{1}{2}\rho_{k + 1}^{(\text{expl.})} - \Delta t\frac{M_{k + 1/2} - M_{k + 3/2}}{2V_{k + 1}}\nonumber\\ &= \newoverline{\rho_{k}^{(\text{expl.})}}^{(k + 1/2)} - \Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{2V_{k}} - \Delta t\frac{M_{k + 1/2} - M_{k + 3/2}}{2V_{k + 1}}. \end{align} \]
From this it follows
\[ \begin{align} &\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg[\gamma_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &- \gamma_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vert_{k + 1/2}\left(z_{k} - z_{k + 1}\right)\Bigg[\alpha_{k}\left(M_{k - 1/2} - M_{k + 1/2}\right) + \beta_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &+\alpha_{k + 1}\left(M_{k + 1/2} - M_{k + 3/2}\right)+\beta_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &- \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}}\frac{\frac{2M_{k + 1/2}}{A_{k + 1/2}} - w^{(n)}_{k + 1/2}\newoverline{\rho_k^{(\text{expl.})}}^{(k + 1/2)} + w^{(n)}_{k + 1/2}\Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{2V_{k}} + w^{(n)}_{k + 1/2}\Delta t\frac{M_{k + 1/2} - M_{k + 3/2}}{2V_{k + 1}}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}\nonumber\\ &= -\frac{w_{k + 1/2}^{(\text{expl.})}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} + \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vert_{k + 1/2}. \end{align} \]
This is the equation to be solved. It is of the form
\[ \begin{align} c_{k - 1/2}M_{k - 1/2} + d_{k + 1/2}M_{k + 1/2} + e_{k + 3/2}M_{k + 3/2} = r_{k + 1/2}. \end{align} \]
For the coefficients, one has
\[ \begin{align} c_{k - 1/2} &= \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\gamma_{k}\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)} + \frac{d\newoverline{\Pi}}{2dz}\vert_{k + 1/2}\left(z_{k} - z_{k + 1}\right)\left(\alpha_{k} + \beta_{k}\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}\right)\nonumber\\ &- \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta tc^{(p)}}\frac{w^{(n)}_{k + 1/2}}{2V_{k}\newoverline{\rho_k^{(n)}}^{(k + 1/2)}},\\ d_{k + 1/2} &= -\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\right)^2\left(\gamma_{k} + \gamma_{k + 1}\right) + \frac{d\newoverline{\Pi}}{2dz}\vert_{k + 1/2}\left(z_{k} - z_{k + 1}\right)\left[\alpha_{k + 1} - \alpha_{k} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\left(\beta_{k + 1} - \beta_{k}\right)\right]\nonumber\\ &- \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}\left[\frac{2}{A_{k + 1/2}} + \frac{\Delta tw^{(n)}_{k + 1/2}}{2}\left(-\frac{1}{V_{k}} + \frac{1}{V_{k + 1}}\right)\right],\\ e_{k + 3/2} &= \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\gamma_{k + 1}\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)} - \frac{d\newoverline{\Pi}}{2dz}\vert_{k + 1/2}\left(z_{k} - z_{k + 1}\right)\left(\alpha_{k + 1} + \beta_{k + 1}\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}\right)\nonumber\\ &+ \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta tc^{(p)}}\frac{w^{(n)}_{k + 1/2}}{2V_{k + 1}\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}, \end{align} \] \[ \begin{align} r_{k + 1/2} &= -\frac{w_{k + 1/2}^{(\text{expl.})}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} + \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vert_{k + 1/2}\nonumber\\ &- \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}}\frac{w^{(n)}_{k + 1/2}\newoverline{\rho_k^{(\text{expl.})}}^{(k + 1/2)}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}. \end{align} \]
Linear systems of equations with tridiagonal matrices can be solved using the Thomas algorithm. From the result for $M_{k + 1/2}$ the other unknowns can be determined as follows:
\[ \begin{align} \rho^{(n + 1)}_k &= \rho^{(\text{expl.})}_k - \Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{V_k},\\ \newtilde{\theta}^{(n + 1)}_k &= \newtilde{\theta}^{(\text{expl.})}_k - \Delta t\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}}{V_k},\\ \Pi_k^{'(n+1)} &= \Pi_k^{'(\text{expl.})} + \gamma_kV_k\left(\newtilde{\theta}_k^{(n + 1)} - \newtilde{\theta}_k^{(\text{expl.})}\right),\\ w^{(n + 1)}_{k + 1/2} &\stackrel{\href{#eq:w_from_ml}{\text{Eq. (34.116)}}}{=} \frac{\frac{2M_{k + 1/2}}{A_{k + 1/2}} - \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}w^{(n)}_{k + 1/2}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}. \end{align} \]
If one does not want to use a hydrostatic background state, one only has to make the substitutions
\[ \begin{align} \newoverline{\theta}_k &= 0,\\ \newoverline{\Pi}_k &= 0,\\ w_{k + 1/2}^{(\text{expl.})} &\to w_{k + 1/2}^{(\text{expl.})} - \Delta t\frac{d\newoverline{\phi}}{dz}\vert_{k + 1/2} \end{align} \]
The forward-backward procedure was presented in Sect. 25.4.1. It does not increase the accuracy, but it does increase the stability of the model. In the dynamic core to be formulated here, the horizontal momentum equation is first integrated forward, then the scalar equations are solved with the new values of the velocities.
The performance of a HEVI method with separation of horizontal fast (acoustic) and slow (advective) components can be generally calculated. Let $\tau$ be the acoustic time step and $T$ be the advective time step. We have $T\geq\tau$. However, when processing a complete time step, higher-order methods calculate the corresponding terms several times. $n_\tau$ is the number of calculations of the acoustic components and $n_T$ is the number of calculations of the advective components. This allows effective time steps
\[ \begin{align} \tau_\text{eff} &\coloneqq \frac{\tau}{n_\tau},\\ T_\text{eff} &\coloneqq \frac{T}{n_T} \end{align} \]
to be defined. The larger these time steps are, the better for the performance of the model. The speed $c$ of the model, i.e. the integrated time $t$ divided by the time required for it, is calculated by
\[ \begin{align} c = \frac{t}{\frac{t}{T_\text{eff}}t_T + \frac{t}{\tau_\text{eff}}t_\tau} = \frac{1}{\frac{1}{T_\text{eff}}t_T + \frac{1}{\tau_\text{eff}}t_\tau} = \frac{\tau_\text{eff}T_\text{eff}}{\tau_\text{eff}t_T + T_\text{eff}t_\tau}. \end{align} \]
Here $t_T$ is the time required for an advective time step and $t_\tau$ is the time required for an acoustic time step. With the notation $t_T\equiv\alpha t_\tau$ this becomes
\[ \begin{align} c = \frac{\tau_\text{eff}T_\text{eff}}{\tau_\text{eff}\alpha t_\tau + T_\text{eff}t_\tau} = \frac{\tau_\text{eff}}{t_\tau}\frac{T_\text{eff}}{\alpha\tau_\text{eff} + T_\text{eff}}. \end{align} \]
The basic time-stepping procedure of the model should be explicit. The explicit Euler method is only first order and therefore not precise enough. A multi-step process produces unphysical modes, which is also disadvantageous. Therefore, a predictor-corrector method is used.
The basic architecture of the time step method is that of the predictor-corrector method. In addition, some modifications are made to stabilize and increase efficiency. The basic procedure for each of the three sub-steps is as follows:
First, the velocity tendencies of the explicit part are calculated. In the second RK step, the final tendency is calculated as the arithmetic mean of the two steps.
Now the horizontal divergences are calculated using the scalar values of the predictor step and the new values of the horizontal velocity components; diabatic effects are also taken into account. In the second RK step, the final tendency is calculated as the arithmetic mean of the two steps.
Using the procedure described in Sect. 34.5.2.1, the new values of the vertical velocity and the scalar fields are calculated.
Finally, the new values of the mass densities of the tracers are calculated using the continuity equations, see Sect. 36.1.
The atmosphere is global. A regional model is still practical for conducting experiments and operational runs with increased resolution. The reason for choosing the hexagonal grid for the global model developed in this book is its numerical properties, which are similar to the quadrangular grid, combined with quasi-uniformity and orthogonality. Therefore, one could also implement a local model on a hexagonal grid. The disadvantage, however, is the effort that the grid generation process entails as well as the difficulty of parallelization and the possible additional effort that is necessary to achieve satisfactory performance. Therefore, in this section, the previous findings in this chapter will be generalized to a regional longitude-latitude grid, since the pole problem does not occur in a regional model anyway.
The term nesting describes the generation of boundary conditions for a regional model. Let
\[ \begin{align} \psi = \psi\left(\mathbf{x}, z, t\right) \end{align} \]
be a model state, where $\mathbf{x}$ denotes the horizontal coordinates. A state is now generated from a model with a larger domain, which is referred to here as the boundary model, for example a global model
\[ \begin{align} \phi = \phi\left(\mathbf{x}, z, t\right), \end{align} \]
which is obtained by spatial and temporal interpolation of the boundary model onto the grid of the target model. Here, linear interpolation in time between the time steps of the boundary model is usually used, i.e.
\[ \begin{align} \phi\left(\mathbf{x}, z, t\right) = \frac{t_{n+1} - t}{t_{n+1} - t_n}\phi\left(\mathbf{x}, z, t_n\right) + \frac{t - t_n}{t_{n+1} - t_n}\phi\left(\mathbf{x}, z, t_{n + 1}\right). \end{align} \]
At each time step of the target model, the model state $\psi$ is moved a little towards $\phi$ at the edges of the domain, i.e.
\[ \begin{align} \psi_\text{neu}\left(\mathbf{x}, z, t\right) = \left(1 - f\left(\mathbf{x}\right)\right)\psi_\text{alt}\left(\mathbf{x}, z, t\right) + f\left(\mathbf{x}\right)\phi\left(\mathbf{x}, z, t\right). \end{align} \]
Here $f = f\left(\mathbf{x}\right)$ is a function that describes the admixture of the boundary state. The simplest approach would be to set this function to one at the very edge of the domain and zero otherwise. However, this can create relatively sharp gradients, which can emit unphysical waves into the model area. A more sensible choice, for example, is the ansatz
\[ \begin{align} f = f\left[d\left(\mathbf{x}\right)\right] = \cos\left[\text{min}\left(\frac{d\pi}{2n\Delta x},\frac{\pi}{2}\right)\right]^2. \end{align} \]
Here, $d\left(\mathbf{x}\right)$ is the distance of the location $\mathbf{x}$ from the edge of the domain and $n$ is the number of grid cells within which the weight of the boundary conditions drops to zero. These procedures are also known as Nudging.