Ziel dieses Kapitels ist es, herzuleiten, wie differenzielle und algebraische Operatoren und Randbedingungen zu diskretisieren sind. Hierzu betrachtet man die Atmosphäre $A \subseteq \mathbb{R}^3$ als offene Menge und verwendet als Gleichungssystem die reversiblen Gleichungen trockener Luft mit den Randbedingungen
\[ \begin{align} \mathbf{v}\cdot\mathbf{n} = 0 \end{align} \]
am oberen Rand und
\[ \begin{align} \mathbf{v} = \mathbf{0} \end{align} \]
am unteren Rand.
Als prognostische Variablen werden
gewählt. Die Hamilton-Funktion $H$ schreibt man dann als Funktion von $\rho, \newtilde{s}$ und $\mathbf{v}$,
\[ \begin{align} H = H\left(\rho, \newtilde{\theta}, \mathbf{v}\right). \end{align} \]
Die prognostischen Gleichungen für diesen Fall wurden in Abschn. 10.1.2 hergeleitet, sie lauten
\[ \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} \]
Die Funktion $f$ in Glg. (34.7) ist die thermische Zustandsgleichung idealer Gase, $\etabi\coloneqq\zetabi + \mathbf{f}$ ist die absolute Vorticity.
Die Poisson-Klammern sind globale Integrale. Man definiert daher einen diskreten Integraloperator gemäß der Ersetzung
\[ \begin{align} \int_A\psi\left(\mathbf{r}\right)d^3r \to \sum_{i = 1}^N\psi_iV_i, \end{align} \]
hierbei sind $\psi\left(\mathbf{r}\right)$ ein integrierbares Skalar- oder Vektorfeld, $\psi_i$ der Wert der diskretisierten Version von $\psi$ in Gitterbox $i$, $V_i$ das Volumen der Gitterbox $i$ sowie $N$ die Anzahl der Gitterboxen. Die Poisson-Klammern werden gemäß dem Grundansatz
\[ \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} \]
diskretisiert. Zunächst wird das Hamilton-Funktional benötigt. Man diskretisiert Glg. (10.81) in der 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} \]
Man teilt dies auf in zwei Anteile
\[ \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} \]
Den Anteil der kinetischen Energie kann man noch weiter aufteilen:
\[ \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} \]
Eine alternative Formulierung der kinetischen Energie ist
\[ \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} \]
Dabei sind $\newtilde{\rho_{c, k}}^{(e)}$ ein Mittelungsoperator an die Kanten und $\newtilde{\rho_{c, k}}^{(k + 1/2)}$ ein Mittelungsoperator an die Punkte, an denen $w$ platziert ist. Diese Operatoren müssen
\[ \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} \]
erfüllen, um die Massenerhaltung nicht zu gefährden.
Glg. (A.122) diskretisiert man in der 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} \]
Die Funktionalableitung ist diejenige diskrete Funktion, welche dies für alle Hilfsfunktionen $\mathbf{h}_{c, k}$ erfüllt. $\mathbf{h}$ ist hierbei ein Tupel von Zustandsgrößen. Es kann auch über andere Volumina als $V_{c, k}$ summiert werden.
Für die Funktionalableitungen nach den thermodynamischen Größen erhält man
\[ \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} \]
Die Glg.en (34.16) - (34.17) kann man in die 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} \]
umschreiben. Dabei wurden
\[ \begin{align} V_{e, k} = \frac{l_ed_e\Delta z_{e, k}}{2} \end{align} \] und
\[ \begin{align} V_{c, k + 1/2} = A_c\Delta z_{c, k + 1/2} \end{align} \]
eingesetzt. Dies impliziert
\[ \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} \]
da die Funktionalableitung auch an den Grenzen der Boxen ausgewertet und anschließend mittels Skalarproduktes auf die Mittelpunkte zurückinterpoliert werden kann.
Aus Glg. (34.13) folgt
\[ \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} \]
Die Summe über Zellen kann als eine Summe über Kanten umgeschrieben werden
\[ \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} \]
Der Vorfaktor Zwei steht für die beiden betrachteten Raumrichtungen. Die beiden Zellen, welche an die Kante $e$ grenzen, wurden mit $c_0$ und $c_1$ bezeichnet. Somit gilt
\[ \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} \]
Analog folgt aus Glg. (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} \]
Dabei wurde eine Indexverschiebung vorgenommen. Mit $V_{c, k + 1/2} = A_c\Delta z_{c, k + 1/2}$ erhält man
\[ \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} \]
Somit folgt
\[ \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} \]
Verwendet man den modifizierten Operator $H_\text{kin}^\star\left(\mathbf{v}\right)$ aus Glg. (34.15), erhält man
\[ \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} \]
Glg. (34.27) verallgemeinert man mit
\[ \begin{align} V_{e, k} \coloneqq A_{e,k}\frac{d_{e, k}}{2} \end{align} \]
wie folgt auf die deep atmosphere:
\[ \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} \]
Nun schreibt man die Summe über Zellen wieder als Summe über Kanten auf:
\[ \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} \]
Somit gilt auch in der tiefen Atmosphäre
\[ \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} \]
Analog verfährt man mit Glg. (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} \]
Mit $V_{c, k + 1/2} = A_c\Delta z_{c, k + 1/2}$ erhält man
\[ \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} \]
Somit folgt
\[ \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} \]
Die Glg.en (34.33) - (34.34) gelten auch in der deep atmosphere.
In Abschn. 26.7 wurde von der Forderung $\mathbf{v}\cdot\nabla\psi \hastobe \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v}$ ausgegangen, um eine Diskretisierung des Skalarproduktes herzuleiten, dabei wurden die Diskretisierungen von „ Skalarfeld x Vektorfeld“, Gradient und Divergenz vorausgesetzt. In Abschn. 33.4.2 wurde das Skalarprodukt auf geländefolgende Koordinaten verallgemeinert. Hieraus konnte in Abschn. 33.4.3 die vertikale kontravariante Maßzahl rekonstruiert werden. Dies wurde in Abschn. 33.4.4 genutzt, um die Form der Divergenz in geländefolgenden Koordinaten abzuleiten. Nun muss noch der Gradient auf geländefolgende Koordinaten verallgemeinert werden. Hierzu verwendet man
\[ \begin{align} \mathbf{v}\cdot\nabla\psi \equiv \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v} \end{align} \]
als die den Gradienten festlegende Gleichung.
Aus Glg. (33.69) folgen
\[ \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} \]
und
\[ \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} \]
Zieht man Glg. (34.43) von Glg. (34.44) ab, erhält man
\[ \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{Glg. }\href{ch-32-kinematik.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} \]
Die rot markierten Terme treten nur auf der linken Seite auf, jedoch nicht auf der rechten. $\nabla\psi$ lässt sich also in geländefolgenden Koordinaten nicht mehr so definieren, dass diese Gleichung erfüllt ist. Man könnte diese Terme eliminieren, indem man die kontravarianten Flussdichten nur über eine einseitige vertikale Rekonstruktion berechnet, was ein Widerspruch zur Forderung der Massenerhaltung wäre, da in diesem Fall der Fluss durch eine Fläche davon abhängen würde, von welcher Seite man die Fläche betrachtet. Das Problem kann also nicht exakt gelöst werden.
Nichtsdestotrotz kann Erhaltung der Gesamtenergie immernoch erreicht werden. Hierfür ist es notwendig, dass sich die Identität
\[ \begin{align} \int_A\psi\nabla\cdot\mathbf{v}d^3r = -\int_A\mathbf{v}\cdot\nabla\psi d^3r \end{align} \]
auf die Diskretisierung überträgt. Für diese Rechnung müssen auch die Orientierungen $\gamma_{e(c)}$ der Vektoren berücksichtigt werden. Es wird
\[ \begin{align} \gamma_{e(c)} \coloneqq \begin{cases} 1,\text{ falls der Vektor bei }e\text{ in Bezug auf }c\text{ nach außen zeigt,}\\ -1,\text{ sonst} \end{cases} \end{align} \]
definiert. Somit erhält man
\[ \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} \]
Mit einer Indexverschiebung folgt
\[ \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} \]
Der Faktor 2 steht für die beiden Raumrichtungen in der Horizontalen. Hieraus ergibt sich die Diskretisierung
\[ \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} \]
des Gradienten in geländefolgenden Koordinaten.
Auch in der deep atmosphere ist in geländefolgenden Koordinaten die diskretisierte Produktregel des Nabla-Operators lokal nicht gültig. Hier gilt nach Glg. (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} \]
Bis auf die Korrekturterme der geländefolgenden Koordinaten ist die Produktregel des Nabla-Operators auch hier lokal gültig, es gilt weiterhin die Herleitung in Abschn. 26.7.5. Darüber hinaus verschwindet auch hier ein globales Integral über eine Divergenz, die Begründung hierfür ist dieselbe wie in der shallow atmosphere. Man kann also die partielle Integration auf die Terme der geländefolgende Korrektur beschränken:
\[ \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} \]
Daher gilt auch in der deep atmosphere
\[ \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} \]
wobei hier nur diese Reihenfolge der Mittelungsoperatoren zulässig ist.
\[ \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} \]
Da die Hamilton-Funktion nicht explizit von der Zeit abhängt, ist sie gleich der Gesamtenergie. Setzt man in Glg. (10.66) $F = H$ ein, erhält man
\[ \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} \]
Weiterhin gelten
\[ \begin{align} \lbrace H, Z_a, H\rbrace \stackrel{\href{ch-09-anwendung-des-hamilton-formalismus-auf-d.html#eq:poisson_bracket_atm_gen_deriv_5}{\text{Glg. (10.54)}}}{=} 0, & {} & \lbrace H, H\rbrace_\rho \stackrel{\href{ch-09-anwendung-des-hamilton-formalismus-auf-d.html#eq:poisson_bracket_atm_gen_deriv_3}{\text{Glg. (10.50)}}}{=} 0, & {} & \lbrace H, H\rbrace_{\newtilde{s}} \stackrel{\href{ch-09-anwendung-des-hamilton-formalismus-auf-d.html#eq:poisson_bracket_atm_gen_deriv_4}{\text{Glg. (10.51)}}}{=} 0. \end{align} \]
Dies impliziert Energieerhaltung
\[ \begin{align} H = \text{const.} \end{align} \]
Um Energieerhaltung auf die räumliche Diskretisierung zu übertragen muss man also drei Voraussetzungen beachten:
Der erste Punkt entspricht der Tatsache, dass die generalisierte Coriolis-Kraft keine Arbeit leistet, was bereits in Kap. 28 nachgewiesen wurde. Die zweite und dritte Aussage haben rechnerisch die gleiche Voraussetzung. Hierzu addiere man die Glg.en (10.78) - (10.79), was auf
\[ \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} \]
führt. Die Integrale über die Terme A und B heben sich jeweils gegenseitig auf. Dabei wurden die folgenden beiden Tatsachen benutzt:
Beide Tatsachen gelten auch in der Diskretisierung. Die Nachweise dafür wurden bereits erbracht:
Auf einen weiteren Aspekt sei hier hingewiesen: diskretisiert man Glg. (34.58), kommt $\theta$ nur an den Kanten des Gitters vor, nicht in den Zellzentren. Man kann den Werte von $\theta$ an der Kante also frei wählen, solange man dies auch in der Druckgradientbeschleunigung tut. Dies wird in Abschn. 34.4 genutzt werden, um $\theta$ an der Kante so zu berechnen, dass die Genaugikeit des Flussoperators $-\nabla\cdot\left(\rho\theta\mathbf{v}\right)$ erhöht wird. Auch die $\theta-$Erhaltung selbst wird hierdurch nicht verletzt, da weiterhin das, was aus einer Gitterbox herausfließt, in die benachbarte hineinfließt.
Druckgradient- und Schwerebeschleunigung sowie die generalisierte Corioliskraft sind hierduch abgedeckt. Der Gradient der kinetischen Energie jedoch kürzt sich in den Poisson-Klammern heraus, da er zwar kinetische Energie transportiert, jedoch keine erzeugt oder vernichtet. Dies soll hier noch einmal nachvollzogen werden. Um dies einzusehen, notiert man die Impulsgleichung in der Form
\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} = -\nabla k + \mathbf{l}, \end{align} \]
wobei Druckgradient- und Schwerebeschleunigung im Term $\mathbf{l}$ zusammengefasst wurden. Die Corioliskraft kann hier ignoriert werden. Skalare Multiplikation mit $\mathbf{v}$ ergibt
\[ \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} \]
Multipliziert man dies mit der Massendichte $\rho$, erhält man
\[ \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} \]
Die Kontinuitätsgleichung lautet
\[ \begin{align} \frac{\partial\rho}{\partial t} = -\nabla\cdot\left(\rho\mathbf{v}\right). \end{align} \]
Multipliziert man dies mit der spezifischen kinetischen Energie $k$, erhält man
\[ \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} \]
Addiert man die Glg.en (34.60) und (34.62), erhält man
\[ \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} \]
Hieraus ergibt sich, dass nur das globale Integral über $\mathbf{v}\cdot\mathbf{l}$ zur Evolution der gesamten kinetischen Energie der Atmosphäre beiträgt. In die Herleitung sind nur Rechenregeln für $\frac{\partial}{\partial t}$ sowie die bereits als erfüllt eingesehene Regel $\nabla\cdot\left(\psi\mathbf{v}\right) = \mathbf{v}\cdot\nabla\psi + \psi\nabla\cdot\mathbf{v}$ eingegangen. Die Beschleunigung $-\nabla k$ beeinflusst die energetischen Integrale also auch in der Diskretisierung nicht.
Somit ist die in diesem Kapitel entwickelte Ortsdiskretisierung energetisch selbstkonsistent.
Ein hydrostatischer Hintergrundzustand wie in Abschn. 13.7.5 beschrieben kann verwendet werden, ohne die Energieerhaltung zu zerstören. Dabei muss sichergestellt werden, dass
\[ \begin{align} \mathbf{0} = -c^{(p)}\newoverline{\theta}\nabla\newoverline{\Pi} - \nabla\phi \end{align} \]
auch diskret gilt. Die entsprechenden Terme haben sich dann in der Diskretisierung gegenseitig auf.
Dies ist so jedoch nur ohne Orographie der Fall. Sind jedoch Berge vorhanden, heben sich die horizontalen Anteile der Gradienten nicht exakt gegenseitig auf, sodass ein kleiner Fehler in der Energieerhaltung entsteht. Im Zusammenspiel von Orts- und Zeitdiskretisierung überwiegen jedoch die Vorteile des hydrostatischen Hintergrundzustandes. Der wesentlichste Vorteil hierbei ist, dass die horizontalen orogrphischen Korrekturanteile von $-c^{(p)}\theta\nabla\Pi$ und $-\nabla\phi$ zwei Größen sind, die sich fast komplett gegenseitig aufheben, also
\[ \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} \]
Die übrigbleibende Differenz ist größtenteils ein Diskretisierungsfehler. Dieser wird kleiner, wenn man einen hydrostatischen Hintergrundzustand verwendet, da dann lediglich
\[ \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} \]
berechnet werden muss.
Als obere Randbedingung im Rahmen der reversiblen Dynamik verwendet man die kinematische Randbedingung $\mathbf{v}\cdot\mathbf{n} = 0$. Im Modell entspricht dies
\[ \begin{align} w_{c, 1/2} = 0 \end{align} \]
für alle Zellen $c$.
Die untere Randbedingung im Rahmen der reversiblen Dynamik ist die No-slip-Randbedingung $\mathbf{v} = \mathbf{0}$. Da sich dort nur vertikale Vektorkomponenten befinden, entspricht dies
\[ \begin{align} w_{c, N_L + \frac{1}{2}} = 0 \end{align} \]
für alle Zellen $c$.
Im Rahmen des bisher entwickelten Formalismus gibt es drei Möglichkeiten für Ansätze zur Erhöhung der Konvergenzordnung:
Diese verbleibenden Freiheiten sind nicht ausreichend, um die Konvergenzordnung des Gesamtmodells zu erhöhen.
Die Genauigkeit der Diskretisierung kann erhöht werden, indem man $\theta$ nicht über das arithmetische Mittel
\[ \begin{align} \newtilde{\theta}^{(e)} = \frac{\theta_{c_1(e)} + \theta_{c_2(e)}}{2}, \end{align} \]
was zweiter Ordnung ist, an die Kante interpoliert, sondern in komplexerer Art und Weise. Der vertikale Index $k$ wird hier weggelassen, da in diesem Abschnitt von Zweidimensionalität ausgegangen wird und $c_1(e)$, $c_2(e)$ sind wie üblich die beiden Zeilen, welche an die Kante $e$ grenzen.
Zunächst wird dies eindimensional betrachtet. Hierzu geht man von der wahren Ortsabhängigkeit $\theta = \theta\left(x\right)$ aus und nimmt an, dass in den skalaren Gitterpunkten $i$, welche an Orten $x_i = i\Delta x$ positioniert sind, die Werte von $\theta$ bekannt sind, also
\[ \begin{align} \theta_i = \theta\left(x_i\right). \end{align} \]
Es wird nun eine geeignete Approximation für $\newtilde{\theta}^{(e)}$ gesucht, also $\theta$ an der Position $x = \frac{\Delta x}{2}$. Hierzu ersetzt man zunächst die Ortsabhängigkeit $\theta\left(x\right)$ durch $\mu^{(n)}\left(x\right)$, wobei $n$ für die Konvergenzordnung steht.
Für $n = 2$ macht man den 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} \]
Da man zwei lineare Gleichungen für $\mu^{(2)}$ hat, macht man einen Ansatz mit zwei Koeffizienten, also
\[ \begin{align} \mu^{(2)}\left(x\right) &= a + bx. \end{align} \]
Setzt man dies in die Glg.en (34.71) - (34.72) ein, erhält man
\[ \begin{align} a &= \theta\left(0\right),\\ a + b\Delta x &= \theta\left(\Delta x\right). \end{align} \]
Dies führt auf
\[ \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} \]
Um hieraus $\newtilde{\theta}^{(e,2)}$ zu berechnen, wobei die hochgestellte Zwei für die Ordnung steht, muss noch eine Integration über das Zeitintervall $\left[0, \Delta t\right]$ durchgeführt werden, wobei $\Delta t$ die Länge eines Zeitschrittes ist. Man rechnet:
\[ \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} \]
Um eine Approximation dritter Ordnung zu erhalten ($n = 3$), nimmt man die zweite Ableitung von $\theta$ in der Gitterzelle $c_1\left(e\right)$ mit in die Rekonstruktion von $\theta$ auf, also
\[ \begin{align} \mu^{(3)}\left(x\right) &= a + bx + cx^2. \end{align} \]
Hieraus folgen
\[ \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} \]
was sich zu
\[ \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} \]
umformen lässt. Hieraus folgt
\[ \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} \]
Der Wert dieser Funktion an der Kante ($x = \frac{\Delta x}{2}$) lautet somit
\[ \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} \]
mit
\[ \begin{align} \delta_x^2f \coloneqq \frac{f_1 - 2f_0 + f_{-1}}{\Delta x^2} \end{align} \]
als Approximation für die zweite Ableitung von $\theta$ in $x-$Richtung im Zentrum der Zelle $c_1\left(e\right)$.
In zwei Dimensionen werden die eindimensionalen Formeln einfach in jeder Raumrichtung angewandt. Auf viereckigen Gittern, bei denen alle Vektoren parallel zu einer globalen Koordinatenachse liegen, ist dies kein Problem. Auf Gittern, die dieses Kriterium nicht erfüllen, ist dies komplizierter, worum es im nächsten Abschnitt geht.
Auf dem hexagonalen Gitter sind die Vektoren nicht entlang von Koordinatenlinien ausgerichtet, was die Berechnung der zweiten Ableitungen in den Zellzentren verkompliziert. Daher macht man für $\theta$ den Ansatz
\[ \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$ und $y$ werden hierbei auf der Tangentialebene an die Kugel im Zentrum der betrachteten Zelle berechnet. Die $x-$Achse zeige dabei in Richtung der Kante, an der $\newtilde{\theta}^{(e)}$ berechnet werden soll. Man definiert nun den Koeffizientenvektor
\[ \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} \]
Nur $c_4$ ist für die zweite Ableitung in $x-$Richtung im Zentrum der Zelle relevant, es gilt
\[ \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} \]
Das Ziel ist nun, den Vektor $\mathbf{f}$ zu berechnen. Dazu definiert man zunächst den Vektor
\[ \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} \]
der tatsächlichen Werte von $\theta$ in den Zellzentren. Dabei ist $m \geq 6$ die Anzahl der Punkte, die zur Berechnung von $\mathbf{f}$ berechnet werden, dies sind die Zelle, in deren Zentrum die zweite Ableitung berechnet werden soll und alle benachbarten Zellen. Nun definiert man eine Fehlerfunktion
\[ \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} \]
Man definiert nun die Polynomialmatrix
\[ \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} \]
Damit kann man $\psi$ in der 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} \]
notieren. Hieraus folgt
\[ \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} \]
Es gilt also die Äquivalenz
\[ \begin{align} \nabla_\mathbf{k}\psi = 0 \Leftrightarrow \left(P^TP\right)\cdot\mathbf{f} = P^T\cdot\mathbf{s}, \end{align} \]
was man zu
\[ \begin{align} \mathbf{f} = \left[\left(P^TP\right)^{-1}P^T\right]\cdot\mathbf{s} = B\cdot\mathbf{s} \end{align} \]
mit
\[ \begin{align} B \coloneqq \left[\left(P^TP\right)^{-1}P^T\right] \end{align} \]
umformen kann. Hieraus lassen mithilfe von Glg. (34.91) die zweiten Ableitungen in den Zellzentern in $x-$Richtung berechnen.
Es besteht folgender Konflikt zwischen expliziten und impliziten Verfahren:
Hieraus folgt, dass ein schnelles implizites Verfahren ideal für einen dynamischen Kern wäre. Da die zu lösenden Gleichungen jedoch nichtlinear sind, ist ein dreidimensional implizites Verfahren nur über ein iteratives Vorgehen möglich. In diesem Abschnitt wählt man als Kompromiss das HEVI-Konzept, was für horizontally explicit, vertically implicit steht.
Eine wesentliche Einschränkung in Bezug auf den Zeitschritt und somit auch in Bezug auf die Effizienz eines Modells ist das mit Schallwellen verbundene CFL-Kriterium. Daher teilt man das Gleichungssystem häufig in sogenannte fast modes oder divergent modes einerseits und die slow modes, advective modes oder non-divergent modes andererseits ein. Die fast modes umfassen dabei alle mit Schallwellen verbundenen Terme, während die slow modes die restlichen Terme umfassen. Bei Split-explicit-Techniken integriert man nun die fast modes mit einem ca. dreimal kleineren Zeitschritt. Im Falle des Gleichungssystems Glg.en (34.4) - (34.7) käme nur die Advektion des Impulses für den größeren Zeitschritt in Frage. Die Advektionen von Dichte und Entropie sollten hingegen zu den fast modes gezählt werden, um die Divergenzen der Flussidchten nicht in Advektions- und Kompressionsanteile zerlegen zu müssen. Wenn man hohe Anforderungen an die energetische Selbstkonsistenz eines Zeitschrittverfahrens stellt, ist jedoch auch die Abspaltung der Impulsadvektion nicht sinnvoll.
Um dies einzusehen, muss man sich zunächst einmal klarmachen, dass die Advektion kinetischer Energie mit der Kontinuitätsgleichung in Verbindung steht und daher auch zu den fast modes zählt. Dies wurde bereits in Abschn. 34.2.4 eingesehen. Somit müssen auch beide Terme bzw. Gleichungen mit dem gleichen Zeitschritt integriert werden. Zur Wahrung des diskretisierten Analogons der Lamb-Transformation $-\left(\mathbf{v}\cdot\nabla\right)\mathbf{v} = -\nabla k - \left(\nabla\times\mathbf{v}\right)\times\mathbf{v}$ muss dann auch der generalisierte Coriolis-Term und somit das komplette Gleichungssystem mit dem kleineren Zeitschritt integriert werden.
Aufgrund der hier besprochenen Nachteile wird keine Split-explicit-Option in den zu entwickelnden dynamischen Kern eingebaut.
Vertikal propagierende Wellen werden aus Stabilitätsgründen anteilig implizit behandelt. Hierfür notiert man zunächst
\[ \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-erster-hauptsatz-in-der-atmosphäre.html#eq:exner_pressure_diag}{\text{Glg. (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} \]
Die zwei Striche stehen für lineare Entwicklungen gegenüber dem vorherigen Zeitschritt und nicht für Differenzen zum hydrostatischen Hintergrundzustand, welche mit $\theta',\Pi'$ bezeichet werden. Die Ableitungen kürzt man mit
\[ \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} \]
ab. Diese berechnet man in den Zellzentren, wobei man beim Predictor-Schritt den alten Zeitschritt (Zeitschritt $n$) und beim Corrector-Schritt das arithmetische Mittel der Schritte $n$ und $n + 1$ verwendet. Die vertikale Impulsgleichung notiert man in der 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}\vline_{k+1/2}\right).\tag{34.104}\label{eq:wave_vertical_deriv_0} \end{align} \]
Hierbei ist $\epsilon$ das implizite Gewicht, für welches aufgrund der Feststellungen in Abschn. 25.5.1.3 meist $\frac{c^{(v)}}{c^{(p)}}$ eingesetzt wird. $\left(n + 1\right)^\star$ ist ein vorläufiger Wert für den neuen Zeitschritt. Für die impliziten Größen notiert man nun
\[ \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} \]
Die gestrichenen Größen sind die Abweichungen von den expliziten Komponenten, also
\[ \begin{align} \psi'_k \coloneqq \psi_k^{(n + 1)} - \psi_k^{(\text{expl.})} \end{align} \]
für $\psi = \rho, \newtilde{\theta}.$ Hierfür gelten
\[ \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} \]
Hierbei definiert
\[ \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} \]
für $1 \leq k \leq N_L - 1$ die vertikale kovariante mit den Begrenzungsflächen $A_{k + 1/2}$ der Gitterboxen multiplizierte Massenflussdichte, dies ist der Unbekanntenvektor. Setzt man dies in Glg. (34.104) ein, erhält man
\[ \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}\vline_{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}\vline_{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}\vline_{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}\vline_{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}\vline_{k+1/2}\Bigg). \end{align} \]
Die gestrichenen Größen enthalten den Unbekanntenvektor, daher bringt man sie auf die linke Seite:
\[ \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}\vline_{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}\vline_{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}\vline_{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}\vline_{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}\vline_{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}\vline_{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}\vline_{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}\vline_{k+1/2}. \end{align} \]
Nun nimmt man die Ersetzungen
\[ \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} \]
vor. Es folgt
\[ \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}\vline_{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}\vline_{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}\vline_{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}\vline_{k + 1/2}.\tag{34.115}\label{eq:vert_solver_pre_1} \end{align} \]
Auch $w_{k + 1/2}^{(n + 1)}$ hängt von $M_{k + 1/2}$ ab. Aus Glg. (34.110) folgt
\[ \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} \]
Setzt man dies in Glg. (34.115) ein, erhält man
\[ \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}\vline_{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}\vline_{k + 1/2}. \end{align} \]
Für die Dichte $ \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}$ gilt
\[ \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} \]
Hieraus folgt
\[ \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}\vline_{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}\vline_{k + 1/2}. \end{align} \]
Dies ist die zu lösende Gleichung. Diese ist von der 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} \]
Für die Koeffizienten gilt
\[ \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}\vline_{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}\vline_{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}\vline_{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}\vline_{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} \]
Lineare Gleichungssysteme mit Drei-Band-Matrizen können mit dem Thomas-Algorithmus gelöst werden. Aus dem Ergebnis für $M_{k + 1/2}$ lassen sich die anderen Unbekannten wie folgt bestimmen:
\[ \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{Glg. (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} \]
Möchte man keinen hydrostatischen Hintergrundzustand verwenden, muss man nur die Ersetzungen
\[ \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}\vline_{k + 1/2} \end{align} \]
vornehmen.
Das Forward-Backward-Verfahren wurde in Abschn. 25.4.1 dargestellt. Es erhöht nicht die Genauigkeit, jedoch die Stabilität des Modells. Im hier zu formulierenden dynamischen Kern wird die horizontale Impulsgleichung zunächst vorwärts integriert, anschließend werden die skalaren Gleichungen mit den neuen Werten der Geschwindigkeiten gelöst.
Die Performance eines HEVI-Verfahrens mit Separierung von horizontalen schnellen (akutischen) und langsamen (advektiven) Anteilen kann allgemein berechnet werden. Es seien $\tau$ der akustische Zeitschritt, und $T$ der advektive Zeitschritt. Es gilt $T\geq\tau$. Innerhalb der Abarbeitung eines vollständigen Zeitschrittes werden bei Verfahren höherer Ordnung allerdings die entpsprechenden Terme mehrmals berechnet. $n_\tau$ sie die Anzahl der Berechnungen der akustischen Anteile und $n_T$ die Anzahl der Berechnungen der advektiven Anteile. Damit lassen sich effektive Zeitschritte
\[ \begin{align} \tau_\text{eff} &\coloneqq \frac{\tau}{n_\tau},\\ T_\text{eff} &\coloneqq \frac{T}{n_T} \end{align} \]
definieren. Desto größer diese Zeitschritte sind, desto besser für die Performance des Modells. Die Geschwindigkeit $c$ des Modells, also die integrierte Zeit $t$ dividiert durch die dafür benötigte Zeit, berechnet sich durch
\[ \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} \]
Hierbei ist $t_T$ die Zeit, die für einen advektiven Zeitschritt benötigt wird und $t_\tau$ die Zeit, die für einen akustischen Zeitschritt benötigt wird. Mit der Notation $t_T\equiv\alpha t_\tau$ wird dies zu
\[ \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} \]
Das grundlegende Zeitschrittverfahren des Modells soll explizit sein. Das explizite Euler-Verfahren ist nur erster Ordnung und daher nicht genau genug. Ein Multi-Schritt-Verfahren erzeugt unphysikalische Moden, was ebenfalls unvorteilhaft ist. Daher verwendet man ein Predictor-Corrector-Verfahren.
Die grundlegende Architektur des Zeitschrittverfahrens ist die des Predictor-Corrector-Verfahrens. Zusätzlich werden einige Modifikationen zur Stabilisierung und Effizienzerhöhung vorgenommen. Das prinzipielle Vorgehen bei jedem der drei Unterschritte ist folgendes:
Die Atmosphäre ist global. Ein Regionalmodell ist trotzdem praktisch, um Experimente und operationelle Läufe mit erhöhter Auflösung durchführen zu können. Die Begründung für die Auswahl des hexagonalen Gitters für das in diesem Buch entwickelte Globalmodell sind die im Vergleich zum viereckigen Gitter ähnlich guten numerischen Eigenschaften in Verbindung mit Quasi-Uniformität und Orthogonalität. Daher könnte man auch ein lokales Modell auf einem hexagonalen Gitter implementieren. Nachteil ist jedoch der Aufwand, den der Gittergenerierungsprozess mit sich bringt sowie die schwierige Parallelisierbarkeit und der eventuelle Mehraufwand, der nötig ist, um eine zufriendenstellende Performance zu erzielen. Daher sollen in diesem Abschnitt die bisherigen Erkenntnisse dieses Kapitels auf ein regionales Länge-Breite-Gitter verallgemeinert werden, da das Polproblem bei einem regionalen Modell ohnehin nicht auftritt.
Der Begriff Nesting beschreibt die Generierung von Randbedingungen für ein regionales Modell. Sei
\[ \begin{align} \psi = \psi\left(\mathbf{x}, z, t\right) \end{align} \]
ein Modellzustand, wobei $\mathbf{x}$ die horizontalen Koordinaten bezeichnet. Aus einem Modell mit größerer Domain, welches hier als Randmodell bezeichnet wird, beispielsweise einem Globalmodell, generiert man nun einen Zustand
\[ \begin{align} \phi = \phi\left(\mathbf{x}, z, t\right), \end{align} \]
der durch örtliche und zeitliche Interpolation des Randodells auf das Gitter des Zielmodells gewonnen wird. Dabei wird zeitlich üblicherweise linear zwischen den Zeitschritten des Randmodells linearisiert, also
\[ \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} \]
Bei jedem Zeitschritt des Zielmodells wird nun der Modellzustand $\psi$ an den Rändern der Domain ein Stück in Richtung $\phi$ bewegt, also
\[ \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} \]
Hierbei ist $f = f\left(\mathbf{x}\right)$ eine Funktion, die die Beimengung des Randzustandes beschreibt. Der einfachst Ansatz wäre, diese Funktion am äußersten Rand der Domain auf Eins zu setzen und sonst auf Null. Dies kann jedoch relativ scharfe Gradienten erzeugen, welche unphysikalische Wellen ins Modellgebiet emittieren können. Sinnvoller ist beispielsweise ein 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} \]
Hierbei ist ist $d\left(\mathbf{x}\right)$ der Abstand des Ortes $\mathbf{x}$ vom Rand der Domain und $n$ die Anzahl der Gitterzellen, innerhalb derer die Gewichtung der Randbedingungen auf Null abfällt. Diese Verfahren bezeichnet man auch als Nudging.