The term dynamic core is often not precisely defined, but this is what will be done here. For this purpose, various subcomponents are defined:
The reversible part of the dynamic core solves the equations of a dry adiabatic atmosphere without molecular and numerical diffusion terms. For models in which this part of the dynamic core does not exactly comply with the second law of thermodynamics, it is better to call this part primary part of the dynamic core.
The irreversible part of the dynamic core deals with molecular and numerical diffusive flows. These can be heat, mass or momentum flows. For models in which the previously defined part of the dynamic core does not exactly comply with the second law of thermodynamics, it is better to call this part secondary part of the dynamic core. Together with the previously defined part, this part forms the dry dynamic core.
The extended dynamic core handles the advection of tracers and ensures that negative tracer densities do not occur. These can be water components, aerosols and/or trace gases.
First, one inscribes an icosahedron in a sphere. It consists of 20 triangles, so the number $E$ of vertices is given by
\[ \begin{align} E = \frac{3\cdot 20}{5} = 12, \end{align} \]
since each vertex is touched by five polygons. The number $K$ of edges is
\[ \begin{align} K = \frac{3\cdot 20}{2} = 30. \end{align} \]
Now each triangle is subdivided $n$ times into four triangles each, where $n\in \mathbb{N}$. There then exist
\[ \begin{align} D = 20\cdot 4^n \end{align} \]
triangles. The dual grid of the resulting grid is the hexagonal grid. It consists of
\[ \begin{align} N = N_5 + N_6 \end{align} \]
polygons, where $N_5$ denotes the number of pentagons and $N_6$ denotes the number of hexagons. Pentagons only form at the twelve vertices of the icosahedron, so $N_5 = 12$ is independent of $n$. The number of vertices of this grid is
\[ \begin{align} \frac{5N_5 + 6N_6}{3} = D = 20\cdot 4^n, \end{align} \]
since each vertex is touched by three polygons; this is equal to the number of triangles. It thus follows
\[ \begin{align} 5N_5 + 6N_6 = 60 + 6N_6 = 60\cdot 4^n\Rightarrow N_6 = 10\left(4^n - 1\right)\Rightarrow N = 10\left(4^n - 1\right) + 12. \end{align} \]
This is the number of scalar degrees of freedom per model level $N_S^{(H)}$. The number $N_V^{(H)}$ of vector degrees of freedom per model level is given by
\[ \begin{align} N_V^{(H)} = \frac{5N_5 + 6N_6}{2} = \frac{60\cdot 4^n}{2} = 30\cdot 4^n. \end{align} \]
The number of layers $N_L$ is arbitrarily set to
\[ \begin{align} N_L = 2 + 6\cdot n \end{align} \]
The total number of horizontal vectors is $N_L\cdot N_V^{(H)}$, to which the number of vertical vectors $N_S^{(H)}\cdot\left(N_L + 1\right)$ must be added to get the total number
\[ \begin{align} N_V = N_L\cdot N_V^{(H)} + N_S^{(H)}\cdot\left(N_L + 1\right) \end{align} \]
of vectors. For the total number of scalars $N_S$, one has
\[ \begin{align} N_S = N_L\cdot N_S^{(H)} = N_L\left[10\left(4^n - 1\right) + 12\right]. \end{align} \]
On the dual grid (the triangular grid), one has
\[ \begin{align} \prescript{{(D)}}{}{N_S^{(H)}} = 20\cdot 4^n, & {} & \prescript{{(D)}}{}{N_S} = 20\cdot 4^n\left(N_L + 1\right),\\ \prescript{{(D)}}{}{N_V^{(H)}} = 30\cdot 4^n, & {} & \prescript{{(D)}}{}{N_V} = 30\cdot 4^n\left(N_L + 1\right) + 20\cdot 4^n\cdot N_L. \end{align} \]
| index | geographical latitude | geographical longitude |
|---|---|---|
| $0$ | $\frac{\pi}{2}$ | $0$ |
| $1$ | $\arctan\left(\frac{1}{2}\right)$ | $0$ |
| $2$ | $\arctan\left(\frac{1}{2}\right)$ | $\frac{2\pi}{5}$ |
| $3$ | $\arctan\left(\frac{1}{2}\right)$ | $2\frac{2\pi}{5}$ |
| $4$ | $\arctan\left(\frac{1}{2}\right)$ | $3\frac{2\pi}{5}$ |
| $5$ | $\arctan\left(\frac{1}{2}\right)$ | $4\frac{2\pi}{5}$ |
| $6$ | $-\arctan\left(\frac{1}{2}\right)$ | $\frac{2\pi}{10}$ |
| $7$ | $-\arctan\left(\frac{1}{2}\right)$ | $\frac{2\pi}{10} + \frac{2\pi}{5}$ |
| $8$ | $-\arctan\left(\frac{1}{2}\right)$ | $\frac{2\pi}{10} + 2\frac{2\pi}{5}$ |
| $9$ | $-\arctan\left(\frac{1}{2}\right)$ | $\frac{2\pi}{10} + 3\frac{2\pi}{5}$ |
| $10$ | $-\arctan\left(\frac{1}{2}\right)$ | $\frac{2\pi}{10} + 4\frac{2\pi}{5}$ |
| $11$ | $-\frac{\pi}{2}$ | $0$ |
In order to make communication easier, the following conventions are made:
A Layer is a layer of grid cells.
A Level is the edge of a layer.
The numbering starts at zero at the top of the atmosphere. If there are $n$ layers, then there are $n + 1$ levels.
All variables are placed in the middle of the layers, vertical components of vector fields in the levels. Such a vertical lattice structure is called a Lorenz lattice.
The radial component of a vector field points positively upwards.
Let a grid with $n$ layers be given. According to Sect. 24.2, central difference quotients are of second order, while other types of difference quotients are only of first order. Therefore, one specifies for the z-coordinates of the levels that they should lie in the middle of the scalar grid points of the neighboring layers, i.e.
\[ \begin{align} z_\text{Level} = \frac{z_\text{Layer above} + z_\text{Layer below}}{2}.\tag{33.13}\label{eq:pos_level} \end{align} \]
The z coordinates of the scalar grid points ($z_\text{Layer}$) must first be determined. It makes sense to choose all the layers to be as equally thick as possible, so that the resolution is as homogeneous as possible. Near the earth's surface, the grid boxes must follow the orography. Since this leads to additional numerical effort, it is useful not to do this above a certain height $\newtilde{z}$. This quantity is a compromise between accuracy and computing time. The height $\newtilde{z}$ is set indirectly by the user. The user sets the height $H$ of the model atmosphere and the number $N_L$ of layers. In addition, the user specifies a number $1 \leq N_{L, O} \leq N_L$ of layers that should follow the orography. $\newtilde{z}$ is then a horizontal scalar field, which is
\[ \begin{align} \newtilde{z}\left(\phi, \lambda\right) = h\left(\phi, \lambda\right) + \frac{N_{L, O}}{N_L}\left(H - h\left(\phi, \lambda\right)\right) \end{align} \]
defined in terms of the orography $h = h\left(\phi, \lambda\right)$.The idea that $h$ is a function of the horizontal coordinates is incorrect when one thinks of buildings, plants and canyons. However, these structures are small-scale. Only when the horizontal model resolution reaches the order of ten meters does one have to take them into account. From this, one calculates preliminary positions $z_{w, \text{pre}, i}$ of vertical vectors, where $i$ is the layer index:
\[ \begin{align} z_{w, \text{pre}, i} = \begin{cases} H - i\frac{H - \newtilde{z}}{N_L - N_{L, O}}, \text{ if } i < N_{L, O},\\ h + \left(N_L - i\right)\frac{\newtilde{z} - h}{N_{L, O}}, \text{ otherwise.} \end{cases} \end{align} \]
The scalar grid points are then placed in the middle of two adjacent vertical vectors before Eq. (33.13) is applied to finally determine the position of the levels. For the horizontal vector points it is specified that they are also placed vertically in the middle of the two adjacent boxes, i.e.
\[ \begin{align} z_\text{vector, h} = \frac{z_{\text{scalar, source}} + z_{\text{scalar, target}}}{2}. \end{align} \]
The same procedure is followed for dual scalar fields; here the determining quantity is the position $z_{j, i}$ of the dual scalar grid points, where $j$ is the horizontal index and $i$ is again the layer. For this quantity one sets
\[ \begin{align} z_{j, i} = \frac{1}{3}\sum_{j' \in \text{APC}(j)}z_{j', w, i} \end{align} \]
fixed. So it is the mean of the z-positions of the vertical vectors of the three surrounding cells.
The starting point is spherical coordinates with a generalized vertical coordinate $\zeta$ that depends on the horizontal coordinates. The transformation to global coordinates reads
\[ \begin{align} x = r\left(\zeta, \phi, \lambda\right)\cos\left(\phi\right)\cos\left(\lambda\right), & {} & y = r\left(\zeta, \phi, \lambda\right)\cos\left(\phi\right)\sin\left(\lambda\right), & {} & z = r\left(\zeta, \phi, \lambda\right)\sin\left(\phi\right). \end{align} \]
The covariant basis elements of this coordinate system
\[ \begin{align} \mathbf{j}_\eta &= \frac{\partial r}{\partial\zeta}\mathbf{e}_r, & {} & \mathbf{j}_\phi &= r\mathbf{e}_\phi + \frac{\partial r}{\partial\phi}\mathbf{e}_r, & {} & \mathbf{j}_\lambda &= r\cos\left(\phi\right)\mathbf{e}_\lambda + \frac{\partial r}{\partial\lambda}\mathbf{e}_r. \end{align} \]
From this one obtains for the six relevant scalar products of the basis elements
\[ \begin{align} \mathbf{j}_\eta\cdot\mathbf{j}_\eta = \left(\frac{\partial r}{\partial\zeta}\right)^2, & {} & \mathbf{j}_\eta\cdot\mathbf{j}_\phi = \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\phi}, & {} & \mathbf{j}_\eta\cdot\mathbf{j}_\lambda = \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\lambda},\\ \mathbf{j}_\phi\cdot\mathbf{j}_\phi = r^2 + \left(\frac{\partial r}{\partial\phi}\right)^2, & {} & \mathbf{j}_\phi\cdot\mathbf{j}_\lambda = \frac{\partial r}{\partial\phi}\frac{\partial r}{\partial\lambda}, & {} & \mathbf{j}_\lambda\cdot\mathbf{j}_\lambda = r^2\cos^2\left(\phi\right) + \left(\frac{\partial r}{\partial\lambda}\right)^2. \end{align} \]
For the determinant of the metric tensor, expansion along the first row yields
\[ \begin{align} g &= \left|\begin{array}{ccc} \left(\frac{\partial r}{\partial\zeta}\right)^2 & \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\phi} & \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\lambda}\\ \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\phi} & r^2 + \left(\frac{\partial r}{\partial\phi}\right)^2 & \frac{\partial r}{\partial\phi}\frac{\partial r}{\partial\lambda}\\ \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\lambda} & \frac{\partial r}{\partial\phi}\frac{\partial r}{\partial\lambda} & r^2\cos^2\left(\phi\right) + \left(\frac{\partial r}{\partial\lambda}\right)^2\\ \end{array}\right|\nonumber\\ &= \left(\frac{\partial r}{\partial\zeta}\right)^2\left[r^4\cos^2\left(\phi\right) + \textcolor{blue}{r^2\left(\frac{\partial r}{\partial\lambda}\right)^2} + \textcolor{red}{\left(\frac{\partial r}{\partial\phi}\right)^2r^2\cos^2\left(\phi\right)} + \textcolor{green}{\left(\frac{\partial r}{\partial\phi}\right)^2\left(\frac{\partial r}{\partial\lambda}\right)^2 - \left(\frac{\partial r}{\partial\phi}\right)^2\left(\frac{\partial r}{\partial\lambda}\right)^2}\right]\nonumber\\ & - \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\phi}\left[\textcolor{red}{\frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\phi}r^2\cos^2\left(\phi\right)} + \textcolor{magenta}{\frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\phi}\left(\frac{\partial r}{\partial\lambda}\right)^2 - \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\lambda}\frac{\partial r}{\partial\phi}\frac{\partial r}{\partial\lambda}}\right]\nonumber\\ & + \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\lambda}\left[\textcolor{cyan}{\frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\phi}\frac{\partial r}{\partial\phi}\frac{\partial r}{\partial\lambda}} \textcolor{blue}{- \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\lambda}r^2} \textcolor{cyan}{- \frac{\partial r}{\partial\zeta}\frac{\partial r}{\partial\lambda}\left(\frac{\partial r}{\partial\phi}\right)^2}\right]\nonumber\\ &= \left(\frac{\partial r}{\partial\zeta}\right)^2r^4\cos^2\left(\phi\right). \end{align} \]
The terms marked in color cancel each other out. For the functional determinant one obtains
\[ \begin{align} \sqrt{g} = \left|\frac{\partial r}{\partial\zeta}\right|r^2\cos\left(\phi\right). \end{align} \]
If one uses the linearly interpolated level index for $\eta$, it follows
\[ \begin{align} \sqrt{g} \stackrel{\text{model}}{=} \Delta zr^2\cos\left(\phi\right). \end{align} \]
The shallow-atmosphere approximation assumes $r = a = \text{const.}$, where $a$ is the radius of the Earth. As a constant, this factor can be neglected, and one obtains
Let $\psi: A \to \mathbb{R}^n$ with $n \in \lbrace 1, 2, 3\rbrace$ be a doubly continuously differentiable scalar or vector field. The discretization of $\psi$ is written in the form
\[ \begin{align} \psi \to \psi_{h, v}^{(\tau)}, \end{align} \]
here $h \in \mathbb {N}$ is the horizontal index, $v \in \mathbb {N}$ is the vertical index and $\tau \in \mathbb {N}$ is the index defining the time step. If one or more of the indices are not important for consideration, they will be ignored in the notation. If no different definitions are introduced at the respective location, $c$ denotes cells, $e$ edges, $dc$ dual cells and $de$ dual edges. In the vertical, the layers are numbered from $k = 1$ to $k = N_L$, the levels are numbered from $k = \frac{1}{2}$ to $k = N_L + \frac{1}{2}$ (with a distance of one).
The following sets of points play a special role:
the centers of the cells (index $c$)
the midpoints of the dual cells (index $d$)
the midpoints of the dual edges from the primary grid perspective (index $e$)
the midpoints of the dual edges from the perspective of the dual grid (index $de$)
Both the layers and the levels are important. Each element of these eight point sets is assigned a volume $V_{l, k}$ with $l \in \lbrace c, d, e, de\rbrace$. Here one has to distinguish between the primary and dual grid at the edges, even though they are the same points. Normal vectors are denoted by $\mathbf{n}_e$ and $\mathbf{n}_{de}$, respectively.
Let $\psi$ be a scalar field and $\mathbf{v} = \mathbf{u} + w\mathbf{k}$ a vector field. The following types of fields are defined:
scalar field:
\[ \begin{align} \psi_{c, k} \coloneqq \frac{1}{V_{c, k}}\int_{V_{c, k}}\psi d^3r \end{align} \]
horizontal vector field:
\[ \begin{align} u_{e, k} \coloneqq \frac{1}{V_{e, k}}\int_{V_{e, k}}\mathbf{v}\cdot\mathbf{n}_ed^3r \end{align} \]
vertical vector field:
\[ \begin{align} w_{c, k + 1/2} \coloneqq \frac{1}{V_{c, k + 1/2}}\int_{V_{c, k + 1/2}}\mathbf{v}\cdot\mathbf{k}d^3r \end{align} \]
horizontal dual vector field:
\[ \begin{align} u_{e, k + 1/2} \coloneqq \frac{1}{V_{e, k + 1/2}}\int_{V_{e, k + 1/2}}\mathbf{v}\cdot\mathbf{n}_{de}d^3r \end{align} \]
vertical dual vector field:
\[ \begin{align} w_{d, k} \coloneqq \frac{1}{V_{d, k}}\int_{V_{d, k}}\mathbf{v}\cdot\mathbf{k}d^3r \end{align} \]
horizontal rotation field: identical to horizontal dual vector fields
vertical rotation field:
\[ \begin{align} w_{e, k} \coloneqq \frac{1}{V_{e, k}}\int_{V_{e, k}}\mathbf{v}\cdot\mathbf{k}d^3r \end{align} \]
According to Sect. 27.6.1, rotation fields are necessary.
The notation did not distinguish between the set $V \in \mathbb{R}^3$ and its volume $\left|V\right|$, and this will continue to be the case whenever it is clear what is meant.
The slope of a coordinate surface is denoted by
\[ \begin{align} \mathbf{J} \coloneqq \nabla_hz\vert_{\zeta}. \end{align} \]
This is a horizontal vector field. $\mathbf{J}$ is called Slope.
Let $x, a \in \lbrace c, d, e, de\rbrace$ and $\psi_{a, k}$ be a discretized field. Then one defines
\[ \begin{align} \newoverline{\psi_{a, k}}^{(x)} \coloneqq \frac{1}{V_{x, k}}\sum_{y \in a, k'}\psi_{y, k'}\left|V_{x, k}\cap V_{y, k'}\right| \end{align} \]
as an averaging operator. The weighting factors $\frac{\left|V_{x, k}\cap V_{y, k'}\right|}{V_{x, k}}$ correspond to the relative overlap of the corresponding sets. One has
\[ \begin{align} \sum_{y \in a, k'}\frac{\left|V_{x, k}\cap V_{y, k'}\right|}{V_{x, k}} = 1. \end{align} \]
Vertical averaging is also defined analogously. Thus, a single averaging operator is either exclusively horizontal or vertical.
In the vertical, the covariant measures are prognostic variables and therefore do not need to be reconstructed. In the horizontal, however, one needs
\[ \begin{align} \mathbf{v}\cdot\mathbf{e}_h = u_h, \end{align} \]
where $\mathbf{e}_h$ is the horizontal covariant unit vector. To determine $u_h$ from $u$ and the vertical covariant components $w$, one first considers a two-dimensional xz coordinate system. For the sake of simplicity, one assumes a coordinate surface of the form
\[ \begin{align} z\left(x\right) = Jx = \tan\left(\alpha\right)x \end{align} \]
where $\alpha$ is the angle between the coordinate surface and the x-axis. In this case, one has
\[ \begin{align} \mathbf{e}_h = \left(\begin{array}{c} \cos\left(\alpha\right)\\ \sin\left(\alpha\right) \end{array}\right). \end{align} \]
Writing the vector field as
\[ \begin{align} \mathbf{v} = \left(\begin{array}{c} u\\ w \end{array}\right), \end{align} \]
it follows
\[ \begin{align} \mathbf{v}\cdot\mathbf{e}_h = u\cos\left(\alpha\right) + w\sin\left(\alpha\right). \end{align} \]
A horizontal length $L'$ along a coordinate line increases, when tilted by the angle $\alpha$, by the factor
\[ \begin{align} L' = L\cos\left(\alpha\right) \Leftrightarrow L = \frac{L'}{\cos\left(\alpha\right)}. \end{align} \]
In a numerical implementation, it is often more efficient to absorb this horizontal length magnification factor into the vector field component. This yields
\[ \begin{align} \mathbf{v}\cdot\mathbf{e}_h\vert_\text{GK} = u + w\tan\left(\alpha\right) = u + Jw, \end{align} \]
where the index GK stands for geometry-corrected. In three dimensions this relation does not change because the horizontal wind component perpendicular to $u$ does not interfere with these calculations. On a C-grid, one has to average the vertical wind component onto the edge:
\[ \begin{align} v_{h, e, \text{GK}} = u_{e} + J_e\newoverline{w_{c}}^{(e)}. \end{align} \]
In the three-dimensional case, further vertical interpolation of the vertical wind is necessary:
\[ \begin{align} v_{h, e, k, \text{GK}} = u_{e, k} + J_{e, k}\newoverline{\newoverline{w_{c, k + \frac{1}{2}}}^{(e)}}^{(k)} \end{align} \]
One has
\[ \begin{align} \newoverline{\newoverline{w_{c, k + \frac{1}{2}}}^{(e)}}^{(k)} &= \newoverline{\frac{1}{2}\left(w_{c, k + \frac{1}{2}} + w_{o(e), k + \frac{1}{2}}\right)}^{(k)} = \frac{1}{2}\newoverline{\left(w_{c, k + \frac{1}{2}} + w_{o(e), k + \frac{1}{2}}\right)}^{(k)} = \frac{1}{2}\left(\newoverline{w_{c, k + \frac{1}{2}}}^{(k)} + \newoverline{w_{o(e), k + \frac{1}{2}}}^{(k)}\right)\nonumber\\ &= \newoverline{\newoverline{w_{c, k + \frac{1}{2}}}^{(k)}}^{(e)}. \end{align} \]
Thus, the order of the averagings plays no role here.
In the deep atmosphere, only the averaging order $\newoverline{\newoverline{w_{c, k + \frac{1}{2}}}^{(k)}}^{(e)}$ is permitted. The averaging onto the edge $(e)$ is still the arithmetic mean of the two values averaged onto the cell centers. For the vertical averaging onto the centers, the weights from Eq. (26.72) are used:
\[ \begin{align} \newoverline{w_{c, k + \frac{1}{2}}}^{(k)} = \frac{A_{c, k - 1/2}\Delta z_{c, k - 1/2}w_{c, k - \frac{1}{2}} + A_{c, k + 1/2}\Delta z_{c, k + 1/2}w_{c, k + \frac{1}{2}}}{2V_{c, k}} \end{align} \]
Here the vertical areas also depend on the vertical index.
The scalar product turns two vector fields into a scalar field and should therefore also be considered in the context of the averaging operators.
In Sect. 26.7 it was assumed that the discretizations of gradient and divergence are fixed, and the discretization of the three-dimensional scalar product on the C-grid, Eq. (26.66), was derived. Here one proceeds the other way around and presupposes the discretization of the scalar product. One generalizes Eq. (26.66) to terrain-following coordinates as follows:
\[ \begin{align} \left(\mathbf{u}^{(1)}\cdot\mathbf{u}^{(2)}\right)_{c, k} = \sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Delta z_{e, k}u_{e, k}^{(1)}u_{e, k}^{(2)} + \frac{\Delta z_{c, k - 1/2}}{2\Delta z_{c, k}}w_u^{(1)}w_u^{(2)} + \frac{\Delta z_{c, k + 1/2}}{2\Delta z_{c, k}}w_l^{(1)}w_l^{(2)}\tag{33.47}\label{eq:inner_c-grid_3d_shallow_terrain} \end{align} \]
This can be viewed as an averaging operator from the edges to the cell centers:
\[ \begin{align} \newoverline{\mathbf{u}^{(1)}\cdot\mathbf{u}^{(2)}}^{(c)} = \sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Delta z_{e, k}u_{e, k}^{(1)}u_{e, k}^{(2)}\tag{33.48}\label{eq:edge_to_cell_shallow} \end{align} \]
In the deep atmosphere, one makes the substitutions
\[ \begin{align} A_c\Delta z_{c, k} &\to V_{c, k},\\ \Delta z_{e, k}l_{e, k} &\to A_{e, k},\\ d_e &\to d_{e, k} \end{align} \]
The vertical component is based on Eq. (26.72). This leads to
\[ \begin{align} \left(\mathbf{u}^{(1)}\cdot\mathbf{u}^{(2)}\right)_{c, k} = \sum_{e\in c}\frac{A_{e, k}d_e}{2V_{c, k}}u_{e, k}^{(1)}u_{e, k}^{(2)} + \frac{A_{c, k - 1/2}\Delta z_{c, k - 1/2}}{2V_{c, k}}w_u^{(1)}w_u^{(2)} + \frac{A_{c, k + 1/2}\Delta z_{c, k + 1/2}}{2V_{c, k}}w_l^{(1)}w_l^{(2)}. \end{align} \]
Analogously, the averaging operator from the edges to the cell centers reads
\[ \begin{align} \newoverline{\mathbf{u}^{(1)}\mathbf{u}^{(2)}}^{(c, k)} = \sum_{e\in c}\frac{A_{e, k}d_{e, k}}{2V_{c, k}}u_{e, k}^{(1)}u_{e, k}^{(2)}.\tag{33.53}\label{eq:edges2cells_deep} \end{align} \]
In the horizontal, the contravariant measures are prognostic variables and therefore do not need to be reconstructed. In the vertical, however, one needs
\[ \begin{align} \mathbf{v}\cdot\mathbf{e}^{(v)} = v^{(v)}, \end{align} \]
where $\mathbf{e}^{(v)}$ is the vertical contravariant unit vector. To determine $v^{(v)}$ from $w$ and the horizontal contravariant components $u$, one first considers a two-dimensional xz coordinate system; in the y direction the system is homogeneous. For the sake of simplicity, one assumes a coordinate surface of the form
\[ \begin{align} z\left(x\right) = Jx = \tan\left(\alpha\right)x \end{align} \]
where $\alpha$ is the angle between the coordinate surface and the x-axis. In this case, one has
\[ \begin{align} \mathbf{e}^{(v)} = \left(\begin{array}{c} -\sin\left(\alpha\right)\\ \cos\left(\alpha\right) \end{array}\right). \end{align} \]
Writing the vector field as
\[ \begin{align} \mathbf{v} = \left(\begin{array}{c} u\\ w \end{array}\right), \end{align} \]
it follows
\[ \begin{align} \mathbf{v}\cdot\mathbf{e}^{(v)} = -u\sin\left(\alpha\right) + w\cos\left(\alpha\right). \end{align} \]
A horizontal coordinate surface $A'$ increases, when tilted by the angle $\alpha$, by the factor
\[ \begin{align} A' = A\cos\left(\alpha\right) \Leftrightarrow A = \frac{A'}{\cos\left(\alpha\right)}. \end{align} \]
In a numerical implementation, it is often more efficient to absorb this horizontal area magnification factor into the vector field component. This yields
\[ \begin{align} \mathbf{v}\cdot\mathbf{e}^{(v)}\vert_\text{GK} = w - u\tan\left(\alpha\right) = w - Ju, \end{align} \]
where the index GK stands for geometry-corrected. In the three-dimensional case, the product $Ju$ is generalized to a scalar product:
\[ \begin{align} \mathbf{v}\cdot\mathbf{e}^{(v)}\vert_\text{GK} = w - \mathbf{J}\cdot\mathbf{v}_h \end{align} \]
On a C-grid this becomes
\[ \begin{align} v^{(v)}_{c, \text{GK}} = w_{c} - \newoverline{J_eu_e}^{(c)}. \end{align} \]
If one uses an L-grid in the vertical, a further vertical interpolation is necessary.
One has
\[ \begin{align} &\newoverline{\newoverline{J_{e, k}u_{e, k}}^{(c)}}^{(k + \frac{1}{2})} = \newoverline{\sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Delta z_{e, k}J_{e, k}u_{e, k}}^{(k + \frac{1}{2})} = \sum_{e\in c}\frac{l_ed_e}{2A_c}\newoverline{\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}}^{(k + \frac{1}{2})}\nonumber\\ &= \sum_{e\in c}\frac{l_ed_e}{2A_c}\frac{1}{2\Delta z_{c, k + 1/2}}\left(\Delta z_{c, k}\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k} + \Delta z_{c, k + 1}\frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}u_{e, k + 1}\right)\nonumber\\ &= \sum_{e\in c}\frac{l_ed_e}{2A_c}\frac{1}{2\Delta z_{c, k + 1/2}}\left(\Delta z_{e, k}J_{e, k}u_{e, k} + \Delta z_{e, k + 1}J_{e, k + 1}u_{e, k + 1}\right)\nonumber\\ &= \sum_{e\in c}\frac{l_ed_e}{2A_c}\frac{\Delta z_{e, k + 1/2}}{\Delta z_{c, k + 1/2}}\frac{\Delta z_{e, k}J_{e, k}u_{e, k} + \Delta z_{e, k + 1}J_{e, k + 1}u_{e, k + 1}}{2\Delta z_{e, k + 1/2}}\nonumber\\ &= \sum_{e\in c}\frac{l_ed_e}{2A_c}\frac{\Delta z_{e, k + 1/2}}{\Delta z_{c, k + 1/2}}\newoverline{J_{e, k}u_{e, k}}^{(k + 1/2)} = \newoverline{\newoverline{J_{e, k}u_{e, k}}^{(k + 1/2)}}^{(c)}. \end{align} \]
Here too, the order of the averaging operators does not play any role.
In the deep atmosphere,
is defined. Only this order of averaging operators is permitted. Writing this out, one obtains with Eq. (33.53)
\[ \begin{align} \newoverline{\newoverline{J_{e, k}u_{e, k}}^{(c)}}^{(k + \frac{1}{2})} = \frac{1}{2}\left(\sum_{e\in c}\frac{A_{e, k}d_{e, k}}{2V_{c, k}}J_{e, k}^{(1)}u_{e, k}^{(2)} + \sum_{e\in c}\frac{A_{e, k + 1}d_{e, k + 1}}{2V_{c, k + 1}}J_{e, k + 1}^{(1)}u_{e, k + 1}^{(2)}\right). \end{align} \]
Since one now knows the vertical contravariant mass flux density, one can generalize Eq. (26.63) to
\[ \begin{align} \left(\nabla\cdot\mathbf{v}\right)_{c, k} &= \frac{1}{A_c\Delta z_{c, k}}\left(\sum_{e\in c}l_e\Delta z_{e, k}u_{e, k} + A_cv^{(v)}_{c, k - 1/2, \text{GK}} - A_cv^{(v)}_{c, k + 1/2, \text{GK}}\right)\tag{33.67}\label{eq:div_c-grid_shallow_terrain} \end{align} \]
Here, the replacement
\[ \begin{align} w_{k + 1/2} \to v^{(v)}_{c, k + 1/2, \text{GK}} \end{align} \]
was made, and the layer thickness was made dependent on the horizontal position. Substituting Eq. (33.63) into Eq. (33.67), one obtains
\[ \begin{align} \left(\nabla\cdot\mathbf{v}\right)_{c, k} &= \frac{1}{A_c\Delta z_{c, k}}\left[\sum_{e\in c}l_e\Delta z_{e, k}u_{e, k} + A_c\left(w_{c, k - 1/2} - \newoverline{\newoverline{J_{e, k}u_{e, k}}^{(c)}}^{(k - \frac{1}{2})}\right) - A_c\left(w_{c, k + 1/2} - \newoverline{\newoverline{J_{e, k}u_{e, k}}^{(c)}}^{(k + \frac{1}{2})}\right)\right]\nonumber\\ &= \sum_{e\in c}\frac{l_e\Delta z_{e, k}}{A_c\Delta z_{c, k}}u_{e, k} + \frac{w_{c, k - 1/2}}{\Delta z_{c, k}} - \frac{\newoverline{\newoverline{J_{e, k}u_{e, k}}^{(c)}}^{(k - \frac{1}{2})}}{\Delta z_{c, k}} - \frac{w_{c, k + 1/2}}{\Delta z_{c, k}} + \frac{\newoverline{\newoverline{J_{e, k}u_{e, k}}^{(c)}}^{(k + \frac{1}{2})}}{\Delta z_{c, k}}\nonumber\\ &\stackrel{\href{#eq:edge_to_cell_shallow}{\text{Eq. (33.48)}}}{=} \frac{w_{c, k - 1/2} - w_{c, k + 1/2}}{\Delta z_{c, k}} + \sum_{e\in c}\frac{l_e\Delta z_{e, k}}{A_c\Delta z_{c, k}}u_{e, k}\nonumber\\ &- \frac{\newoverline{\sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Delta z_{e, k}J_{e, k}u_{e, k}}^{(k - \frac{1}{2})} - \newoverline{\sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Delta z_{e, k}J_{e, k}u_{e, k}}^{(k + \frac{1}{2})}}{\Delta z_{c, k}}\nonumber\\ &= \frac{w_{c, k - 1/2} - w_{c, k + 1/2}}{\Delta z_{c, k}} + \sum_{e\in c}\frac{l_e\Delta z_{e, k}}{A_c\Delta z_{c, k}}u_{e, k} - \frac{\sum_{e\in c}\frac{l_ed_e}{2A_c}\newoverline{\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}}^{(k - \frac{1}{2})} - \sum_{e\in c}\frac{l_ed_e}{2A_c}\newoverline{\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}}^{(k + \frac{1}{2})}}{\Delta z_{c, k}}\nonumber \end{align} \] \[ \begin{align} &= \frac{w_{c, k - 1/2} - w_{c, k + 1/2}}{\Delta z_{c, k}} + \sum_{e\in c}\frac{l_e\Delta z_{e, k}u_{e, k} - \frac{l_ed_e}{2}\newoverline{\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}}^{(k - \frac{1}{2})} + \frac{l_ed_e}{2}\newoverline{\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}}^{(k + \frac{1}{2})}}{A_c\Delta z_{c, k}}\nonumber\\ &= \frac{w_{c, k - 1/2} - w_{c, k + 1/2}}{\Delta z_{c, k}} + \sum_{e\in c}l_e\frac{\Delta z_{e, k}u_{e, k} - \frac{d_e}{2}\newoverline{\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}}^{(k - \frac{1}{2})} + \frac{d_e}{2}\newoverline{\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}}^{(k + \frac{1}{2})}}{A_c\Delta z_{c, k}}\nonumber\\ &= \frac{w_{c, k - 1/2} - w_{c, k + 1/2}}{\Delta z_{c, k}}\nonumber\\ & + \sum_{e\in c}l_e\frac{\Delta z_{e, k}u_{e, k} - \frac{d_e}{4}\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}u_{e, k - 1} + \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}\right) + \frac{d_e}{4}\left(\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k} + \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}u_{e, k + 1}\right)}{A_c\Delta z_{c, k}}\nonumber \end{align} \]
\[ \begin{align} \Rightarrow\left(\nabla\cdot\mathbf{v}\right)_{c, k} &= \left[\sum_{e\in c}\left(l_e\frac{\Delta z_{e, k}u_{e, k}}{A_c\Delta z_{c, k}}\right) + \frac{w_{c, k - 1/2} - w_{c, k + 1/2}}{\Delta z_{c, k}}\right]\nonumber\\ &- \frac{1}{2}\sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Bigg[\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}u_{e, k - 1} + \textcolor{blue}{\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}}\right)\nonumber\\ &- \left(\textcolor{blue}{\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}u_{e, k}} + \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}u_{e, k + 1}\right)\Bigg].\tag{33.69}\label{eq:div_c-grid_shallow_terrain_concrete} \end{align} \]
The terms in the first line correspond to Eq. (26.63), the terms in the last two equations are the orographic correction terms. These vanish in the case $\mathbf{J} = \mathbf{0}$. The terms marked in blue cancel each other out and can be ignored.
In the deep atmosphere, Eq. (33.67) generalizes to
\[ \begin{align} \left(\nabla\cdot\mathbf{v}\right)_{c, k} &= \frac{1}{V_{c, k}}\left(\sum_{e\in c}A_{e, k}u_{e, k} + A_{c, k - 1/2}v^{(v)}_{c, k - 1/2, \text{GK}} - A_{c, k + 1/2}v^{(v)}_{c, k + 1/2, \text{GK}}\right). \end{align} \]
Substituting Eq. (33.63) here, one obtains
\[ \begin{align} \left(\nabla\cdot\mathbf{v}\right)_{c, k} &= \frac{1}{V_{c, k}}\left(\sum_{e\in c}A_{e, k}u_{e, k} + A_{c, k - 1/2}v^{(v)}_{c, k - 1/2, \text{GK}} - A_{c, k + 1/2}v^{(v)}_{c, k + 1/2, \text{GK}}\right)\nonumber\\ &= \frac{1}{V_{c, k}}\left[\sum_{e\in c}\left(A_{e, k}u_{e, k}\right) + A_{c, k - 1/2}w_{c, k - 1/2} - A_{c, k + 1/2}w_{c, k + 1/2}\right]\nonumber\\ &- \frac{1}{V_{c, k}}\sum_{e\in c}\Bigg[A_{c, k - 1/2}\frac{1}{2}\left(\frac{A_{e, k}d_{e, k}}{2V_{c, k}}J_{e, k}u_{e, k} + \frac{A_{e, k - 1}d_{e, k - 1}}{2V_{c, k - 1}}J_{e, k - 1}u_{e, k - 1}\right)\nonumber\\ &- A_{c, k + 1/2}\frac{1}{2}\left(\frac{A_{e, k}d_{e, k}}{2V_{c, k}}J_{e, k}u_{e, k} + \frac{A_{e, k + 1}d_{e, k + 1}}{2V_{c, k + 1}}J_{e, k + 1}u_{e, k + 1}\right)\Bigg]\nonumber \end{align} \]
\[ \begin{align} \Rightarrow\left(\nabla\cdot\mathbf{v}\right)_{c, k} &= \frac{1}{V_{c, k}}\left[\sum_{e\in c}\left(A_{e, k}u_{e, k}\right) + A_{c, k - 1/2}w_{c, k - 1/2} - A_{c, k + 1/2}w_{c, k + 1/2}\right]\nonumber\\ &+ \frac{1}{V_{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}u_{e, k} + \frac{A_{e, k - 1}d_{e, k - 1}}{4V_{c, k - 1}}J_{e, k - 1}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}u_{e, k} + \frac{A_{e, k + 1}d_{e, k + 1}}{4V_{c, k + 1}}J_{e, k + 1}u_{e, k + 1}\right)\Bigg].\tag{33.71}\label{eq:div_deep} \end{align} \]
The simplification possible in Eq. (33.69) is not applicable here.