the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The upperatmosphere extension of the ICON general circulation model (version: uaicon1.0)
Sebastian Borchert
Guidi Zhou
Michael Baldauf
Hauke Schmidt
Günther Zängl
Daniel Reinert
How the upperatmosphere branch of the circulation contributes to and interacts with the circulation of the middle and lower atmosphere is a research area with many open questions. Inertia–gravity waves, for instance, have moved in the focus of research as they are suspected to be key features in driving and shaping the circulation. Numerical atmospheric models are an important pillar for this research. We use the ICOsahedral Nonhydrostatic (ICON) general circulation model, which is a joint development of the Max Planck Institute for Meteorology (MPIM) and the German Weather Service (DWD), and provides, e.g., local mass conservation, a flexible grid nesting option, and a nonhydrostatic dynamical core formulated on an icosahedral–triangular grid. We extended ICON to the upper atmosphere and present here the two main components of this new configuration named UAICON: an extension of the dynamical core from shallow to deepatmosphere dynamics and the implementation of an upperatmosphere physics package. A series of idealized test cases and climatological simulations is performed in order to evaluate the upperatmosphere extension of ICON.
 Article
(7618 KB)  Fulltext XML

Supplement
(4455 KB)  BibTeX
 EndNote
In climate simulations and numerical weather prediction (NWP), there are ongoing efforts to raise the upper model lid, acknowledging possible influences of middle and upperatmosphere dynamics on tropospheric weather and climate (e.g., Thompson et al., 2002; Scaife et al., 2012; CharltonPerez et al., 2013). The dynamics of the largescale flow in the middle and upper atmosphere is determined, for instance, by the interaction with smallscale gravity waves. These waves, predominantly forced in the troposphere, can propagate vertically until they become unstable and break. As a result of this, a drag is exerted on the atmospheric background flow. This is an important route of momentum flux from the lower atmosphere to the middle and upper atmosphere, which shapes the meridional circulation in the latter regions (e.g., Fritts and Alexander, 2003; Kim et al., 2003). To have a model at hand that allows to study such processes on a wide range of spatial and temporal scales was one of our central motivations for the upperatmosphere extension of the ICON model, which we present in the following.
The ICOsahedral Nonhydrostatic (ICON) general circulation model (Zängl et al., 2015; Dipankar et al., 2015; Heinze et al., 2017; Giorgetta et al., 2018; Crueger et al., 2018), a joint development of the Max Planck Institute for Meteorology (MPIM) and the German Weather Service (DWD), is designed as a unified modeling system to allow simulations from the turbulent scales in largeeddy simulations (LESs) up to climatological timescales. Its triangular horizontal grid provides an almost uniform horizontal resolution, and the nonhydrostatic dynamical core makes it applicable on a wide range of spatial scales. The discretization of the governing equations makes use of a combination of finitedifference and finitevolume methods, e.g., to ensure mass conservation, which is an important property for climate simulations (e.g., Staniforth and Wood, 2008). In addition to a limitedarea mode, ICON offers a gridnesting option (one or twoway nesting) for local refinement up to potentially very high horizontal resolutions. With three distinct packages of physics parameterizations, ICON meets the different needs of climate simulations, NWP and LES. To prepare ICON for simulations with model tops in the lower thermosphere, some extensions of the dynamical core and the physics parameterizations are necessary, as presented in this work.
Among the most important approximations which are applied to the dynamical core of ICON in its standard configuration are the shallowatmosphere approximation and traditional approximation (e.g., Phillips, 1966; White and Bromley, 1995; White et al., 2005; Staniforth and Wood, 2008). These approximations are applied to the mapping of the budget equations on spherical coordinates relative to the center of the Earth. The shallowatmosphere approximation is basically associated with the neglect of terms related to the spherical curvature of the atmosphere as well as variations of the gravitational field by assuming the field strength to be constant (see Thuburn and White, 2013, for a detailed examination of the metrical implications of this approximation). The traditional approximation refers to neglecting the contribution of the horizontal component of the Earth's angular velocity to the Coriolis acceleration. Following the usual terminology, we will call the system of equations with and without the two approximations the shallowatmosphere equations and deepatmosphere equations, respectively. Both approximations are generally applied together in order to satisfy the conservation of the energy, the axial component of angular momentum, and the potential vorticity (e.g., Phillips, 1966; White and Bromley, 1995; Staniforth and Wood, 2003). However, as Tort and Dubos (2014) have shown, it is possible to extend the shallowatmosphere equations in such a way that the full Coriolis acceleration can be retained without violating the conservation principles.
The accuracy of the shallowatmosphere and traditional approximations can be estimated by comparing the magnitude of the terms neglected in the shallowatmosphere equations to the magnitude of the terms that are present in both systems. Such scale analysis has been used, for instance, by White and Bromley (1995), to show that for diabatically driven flows in the tropics and planetaryscale flows the neglected terms of the Coriolis acceleration might reach magnitudes up to about 10 % of the magnitude of key terms of the shallowatmosphere momentum budget. On the other side, normalmode analyses, done by Thuburn et al. (2002a) for the deepatmosphere equations, and by Kasahara (2003) with focus on a Boussinesq model featuring the full Coriolis acceleration, show that the differences in the spatial structure and the frequencies of the energetically most significant modes between the shallow and the deepatmosphere equations are relatively small (with differences in the frequency magnitude being typically less than about 1 %; Thuburn et al., 2002b; Kasahara, 2003). Both, the scale analysis and the normalmode analysis are important tools to figure out the differences between the deep and shallowatmosphere equations. However, the applicability of the results on longterm integrations might be limited. The systematic errors introduced by the approximations, albeit small in magnitude, could accumulate over time and lead to significantly different flow patterns of the model atmosphere. This might be especially important for the largescale circulations of the middle and upper atmosphere. Furthermore, in view of the everincreasing computational power, some approximations that were meaningful under the restrictions of past computer architectures might nowadays lose their justification. Therefore, we decided to expand the dynamical core of ICON by a deepatmosphere option.
Examples for other models that use a deepatmosphere formulation, or offer the option to do so, are the Met Office's Unified Model (UM) (Davies et al., 2005; Wood et al., 2014), the Nonhydrostatic Icosahedral Atmospheric Model (NICAM) (e.g., Tomita and Satoh, 2004), the Ocean Land Atmosphere Model (OLAM) (e.g., Walko and Avissar, 2008), the MCore model by Ullrich and Jablonowski (2012), and the Finite Volume Model (FVM) of the Integrated Forecasting System (IFS) developed at the European Centre for MediumRange Weather Forecasts (ECMWF) (e.g., Smolarkiewicz et al., 2016). An overview of some of these models can be found in Ullrich et al. (2017).
Apart from the dynamics, the physics parameterizations are the second important model pillar that has to be extended for applications including part of the upper atmosphere. The ICON model offers basically three different physics packages: one which has been largely adopted from the ECHAM model intended for climate simulations (e.g., Stevens et al., 2013; Giorgetta et al., 2018; Crueger et al., 2018), another one which is used for NWP (some aspects of which can be found in Zängl et al., 2015), and a third one for LES (Dipankar et al., 2015). The upperatmospherespecific physics parameterizations have been integrated into the ECHAM and NWP packages, but we will focus on the extended ECHAM physics package in the remainder of this work. To avoid confusion in the following, a side note on our terminology may be in order: if we discuss the physics parameterizations, we typically use the attribute “upper atmosphere” to denote that their effects become significant in and above, say, the uppermesosphere–lowerthermosphere region. In contrast, the attribute “deep atmosphere” is used in the context of the modifications of the dynamical core, since they apply to the entire air column, and no more or less welldefined vertical significance threshold can be made out for them. If we address both extensions as a whole, we use again the attribute “upper atmosphere”.
While the model lid is raised and the deepatmospheric dynamics is applied, the model cannot produce physically reasonable results without the extended physics parameterizations. In fact, our experience shows that when the model lid is raised to above ∼85 km, the model without the extended physics would suffer from strong numerical instability, which could only be suppressed by using extremely large and unphysical numerical damping. The reason is related to the characteristics of the upper atmosphere, most importantly the rarefied air and the broader spectrum of incoming solar irradiance, which give rise to some physical phenomena that are negligible in the lower atmosphere and thus not parameterized in the current model but become crucial in maintaining the upperatmospheric dynamics and thermodynamics.
One of such physical phenomena is molecular diffusion of momentum and heat, which is of negligible magnitude in the lower atmosphere compared to turbulent diffusion. In the upper atmosphere, as turbulence dies away and the air molecules are capable of traveling a long distance, molecular diffusion becomes dominant. Besides molecules, the upper atmosphere is also abundant with atoms and radicals produced as a result of photolysis. Chemical heating is the release of heat by recombination reactions between atoms or radicals, and is of particular importance in the upper atmosphere where photolysis products can travel large distances before recombining. Moreover, the higher altitude also means that the upper atmosphere would receive and absorb more solar irradiance in higher frequencies than at lower levels and on the surface. This brings the need of parameterizing the ultraviolet radiation including the Schumann–Runge bands and continuum, and the extreme ultraviolet bands. The solar radiation also acts to ionize the atmosphere, establishing the ionosphere from about 60 km. The electrically charged ions in the ionosphere are then aligned with the magnetic field of the Earth, thus creating a drag and a heating source to the neutral mass flow. Further, the usual assumption of local thermodynamical equilibrium (LTE) does not hold in the upper atmosphere as a consequence of the low collision frequency between particles; therefore, some modification must be made to the longwave radiation parameterization. The heating and the eddy diffusion of momentum and heat generated by the breaking of gravity waves have to be taken into account near the mesopause, too.
In recent decades, several general circulation models (GCMs) have been extended to the middle and upper atmosphere. An early development in this direction is the spectral model described by Miyahara et al. (1993). With an upper boundary at an altitude of about 165 km, it extends into the lower thermosphere, and has been used to study thermal tides, for instance. The influence of gravity wave parameterizations on the circulation in the middle atmosphere has been studied with the MA/ECHAM4 model, an extension of the hydrostatic spectral ECHAM4 up to the mesopause region (Manzini et al., 1997; Manzini and McFarlane, 1998). Later versions of MA/ECHAM provided the basis for the Hamburg Model of the Neutral and Ionized Atmosphere (HAMMONIA), which extends into the thermosphere and includes a package of upperatmospherespecific physics parameterizations, together with an interactive chemistry module (Schmidt et al., 2006). The upperatmosphere physics package implemented into UAICON has been largely adopted from HAMMONIA. A good overview of the physical processes that are typically parameterized can also be found in the description of an extension of the Canadian Middle Atmosphere Model (CMAM) from an upper boundary at an altitude of about 95 to 210 km by Fomichev et al. (2002). In addition, the authors present an examination of how much these processes contribute to the energy and momentum budgets. Another spectral hydrostatic GCM that includes the uppermesosphere–lowerthermosphere region is the Kühlungsborn Mechanistic general Circulation Model (KMCM, e.g., Becker, 2009). It has been used recently to study secondary gravity wave generation in the mesosphere as a result of the breaking of gravity waves originating from the troposphere (Becker and Vadas, 2018). In contrast to the aforementioned spectral GCMs, WACCM (Whole Atmosphere Community Climate Model; see, e.g., Richter et al., 2008 for its version 3), which is based on the Community Atmosphere Model and Community Climate Model (CAM and CCM) of the National Center for Atmospheric Research (NCAR), offers a hydrostatic finitevolume dynamical core. The different versions of WACCM generally extend to about 150 km and share many parameterizations of its upperatmosphere physics package with CMAM and HAMMONIA. The offspring WACCMX (Liu et al., 2010) extends the simulations even deeper into the thermosphere up to about 500 km. Recent developments, such as selfconsistent electrodynamics, the transport of O^{+}, and a modification of the dynamical core to account for the variation of the specific heats and mean molecular weight in the heterosphere (above ∼100 km, say) led to improved simulations of space weather and space climate (WACCMX 2.0, Liu et al., 2018). As a final example, we mention the middle atmosphere GCM developed by Watanabe et al. (2008) and later complemented by physics parameterizations that allow for simulations covering the lower thermosphere up to about 150 km (the Japanese Atmospheric general circulation model for Upper Atmosphere Research (JAGUAR); Watanabe and Miyahara, 2009). It has been employed, for instance, to study the interaction of resolved gravity waves with thermal tides by global simulations of a horizontal, triangularly truncated spectral resolution of T213 and a vertical layer spacing of 500 m throughout the middle atmosphere (Watanabe and Miyahara, 2009). In a subsequent study on the dependence of the resolved gravitywaveborne vertical flux of zonal momentum on the vertical resolution, the authors were able to conduct simulations with a layer spacing down to 200 m (Watanabe et al., 2015).
The outline of the paper is as follows: in Sect. 2, we describe the modifications and extensions to the dynamical core and to the physics parameterizations, followed by a presentation of results from idealized test cases and climate simulations for the evaluation of the new implementations in Sect. 3. We close with a conclusion in Sect. 4.
2.1 Deepatmosphere dynamics
The dynamical core of the standard configuration of ICON makes use of the shallowatmosphere approximation, which mainly consists in simplifying the governing equations measured in a spherical coordinate system in the following way: the radial distance of an air parcel to the center of the Earth r is approximated by the radius of the Earth a, and metrical terms which result from the unit vectors of the coordinate system to be functions of position are neglected. In addition, the traditional approximation is applied, by which the acceleration due to the horizontal component of the angular velocity of the Earth is neglected (Phillips, 1966; Staniforth and Wood, 2003). For atmospheric models having a model top below, say, the mesopause region at an altitude of about 70 to 100 km, the shallowatmosphere approximation is likely a very good approximation (e.g., Ullrich et al., 2014), albeit some adverse impacts might exist, for instance, in the tropics where the cosine of latitude is of order 1, questioning the neglect of the nontraditional part of the Coriolis acceleration to some extent (White and Bromley, 1995; White et al., 2005). If the model top is raised into the lower thermosphere, the systematic errors introduced by the shallowatmosphere approximation might start to outweigh its benefits, especially on a climatological timescale. So we think the extension of the dynamical core from shallow to deepatmosphere dynamics is an important component for UAICON.
2.1.1 Model equations
We restrict our considerations to the dry atmosphere hereafter in order to focus on the particular aspects of the deepatmosphere dynamics. The budget equations for the momentum, mass, and heat are rearranged to the following set of prognostic equations (White et al., 2005; Zängl et al., 2015):
and the equation of state reads
where $\mathit{v}={\mathit{v}}_{\mathrm{h}}+w{\mathit{e}}_{\mathrm{r}}$ is the wind vector split in horizontal and radial components, Ω is the angular velocity vector of the Earth, and r denotes the position relative to the center of the Earth. In addition, $\mathit{\pi}=(p/{p}_{\mathrm{00}}{)}^{\mathrm{R}/{c}_{p}}$ denotes the Exner pressure relative to the reference pressure p_{00}=1000 hPa, $\mathit{\theta}=T/\mathit{\pi}$ is the potential temperature, ϕ_{g} denotes the gravitational potential, and c_{p}=1004.64 $\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{\mathrm{1}}$, c_{v}=717.6 $\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{\mathrm{1}}$, and $\mathrm{R}={c}_{p}{c}_{v}$ are the specific heat capacities at constant pressure and volume, and the gas constant for dry air, respectively. Finally, F denotes velocity tendencies due to dissipative processes or other parameterized processes, and Q denotes Exner pressure tendencies due to diabatic processes.
We have added the centrifugal acceleration $\mathbf{\Omega}\times (\mathbf{\Omega}\times \mathit{r})$ in Eq. (1), since it is of importance for a test case presented in Sect. 3.1.1. In the standard configuration of ICON, where this term is absent, we have apparent gravity $\mathbf{\nabla}{\mathit{\varphi}}_{\mathrm{g}+\mathrm{c}}$ (true gravity plus centrifugal acceleration) on the righthand side of Eq. (1), to which the spherical–geopotential approximation is applied: in the case of Eq. (1), the vertical unit vector of the coordinate system, in which the atmospheric flow is measured, is defined as the normal vector of the spherical equipotential surfaces ${\mathit{e}}_{\mathrm{r}}=\mathbf{\nabla}{\mathit{\varphi}}_{\mathrm{g}}/\left\mathbf{\nabla}{\mathit{\varphi}}_{\mathrm{g}}\right$. In the case of apparent gravity, we continue to write ${\mathit{e}}_{\mathrm{r}}=\mathbf{\nabla}{\mathit{\varphi}}_{\mathrm{g}+\mathrm{c}}/\left\mathbf{\nabla}{\mathit{\varphi}}_{\mathrm{g}+\mathrm{c}}\right$, but the equipotential surfaces are still approximated by spherical surfaces, neglecting the oblateness caused by the centrifugal acceleration (Gill, 1982; White et al., 2005; Vallis, 2006; Staniforth and Wood, 2008).
The advection term in Eq. (1) is expressed in the socalled 2D vectorinvariant formulation (e.g., Phillips, 1966; Sadourny, 1972; Vallis, 2006):
where ${\mathit{\omega}}_{\mathrm{h}}={\mathbf{\nabla}}_{\mathrm{h}}\times {\mathit{v}}_{\mathrm{h}}=\mathit{\zeta}{\mathit{e}}_{\mathrm{r}}$ is the vertical component of the relative vorticity. During the development of ICON, it turned out that it is advantageous to separate the Exner pressure into a hydrostatically balanced part and a deviation from it $\mathit{\pi}={\mathit{\pi}}_{\mathrm{0}}\left(r\right)+{\mathit{\pi}}^{\prime}$, where ${c}_{p}{\mathit{\theta}}_{\mathrm{0}}\mathbf{\nabla}{\mathit{\pi}}_{\mathrm{0}}\mathbf{\nabla}{\mathit{\varphi}}_{\mathrm{g}}=\mathrm{0}$. Thus, the righthand side of Eq. (1) reads
with ${\mathit{\theta}}^{\prime}=\mathit{\theta}{\mathit{\theta}}_{\mathrm{0}}\left(r\right)$. The background temperature profile of the hydrostatically balanced part is defined by (Zängl, 2012)
where ${z}_{\mathrm{g}}=z/(\mathrm{1}+z/a)$ denotes the geopotential height, and T_{sl}=288.15 K as well as T_{str}=213.15 K are characteristic temperature values at sea level (z_{g}=0) and in the stratosphere, respectively. H_{scal}=10 000 m is a characteristic (geopotential) temperature scale height. The geopotential height follows as “natural” height measure from the gravitational acceleration (Gill, 1982; Wood and Staniforth, 2003; Wood et al., 2014; Ullrich et al., 2017):
with the geopotential ${\mathit{\varphi}}_{\mathrm{g}}=g{a}^{\mathrm{2}}/r+{\mathit{\varphi}}_{\mathrm{g},\mathrm{0}}$, where ϕ_{g,0} denotes some constant, and g=9.80665 m s^{−2} is the mean gravitational acceleration at sea level (this is actually the magnitude of apparent gravity, but since the contribution of the centrifugal acceleration is relatively small, we use this value also for the true gravity). The values of z_{g} range from ${lim}_{z\to a}{z}_{\mathrm{g}}=\mathrm{\infty}$ to ${lim}_{z\to \mathrm{\infty}}{z}_{\mathrm{g}}=a$. The profile (7) allows to integrate the hydrostatic equilibrium analytically. Note that if the centrifugal acceleration is treated in its explicit form $\mathbf{\Omega}\times (\mathbf{\Omega}\times \mathit{r})$, it is not taken into account in the hydrostatic equilibrium outlined above.
In ICON, governing Eqs. (1) to (3) are measured in a spherical coordinate system with local unit vectors e_{𝔱}(r), e_{𝔫}(r), and e_{r}(r) (forming a righthanded system in this order), where e_{𝔱} and e_{𝔫} are the horizontal unit vectors tangential and normal to an edge separating two adjacent triangular cells, and e_{r} is defined in the center of the bottom and top surfaces of a cell (see Fig. 1)^{1}. In order to find expressions for the projection of the momentum equation (Eq. 1) onto the unit vectors of the local edge coordinate system and e_{r}, we can make use of the projections in the standard geographical coordinate system, with its horizontal unit vectors in meridional (south–north) and zonal (west–east) directions, e_{φ} and e_{λ} (e.g., Zdunkowski and Bott, 2003). The Coriolis and centrifugal accelerations are excluded from the following consideration, since their projections onto the unit vectors of the local coordinate systems of the triangular grid can be determined directly. Each edge is part of a great circle of the Earth, which can be regarded as an imaginary Equator. Rotating the geographical coordinate system accordingly and in such a way that e_{𝔱} is parallel to ${\mathit{e}}_{{\mathit{\lambda}}^{\prime}}$ at the considered location (where the prime indicates this coordinate transformation), we would find e_{𝔫} being parallel to ${\mathit{e}}_{{\mathit{\phi}}^{\prime}}$ (the radial unit vectors of both systems are parallel anyway). The governing equations for the velocity components $\mathit{v}=u{\mathit{e}}_{\mathit{\lambda}}+v{\mathit{e}}_{\mathit{\phi}}+w{\mathit{e}}_{\mathrm{r}}$ (or here, $\mathit{v}={u}^{\prime}{\mathit{e}}_{{\mathit{\lambda}}^{\prime}}+{v}^{\prime}{\mathit{e}}_{{\mathit{\phi}}^{\prime}}+w{\mathit{e}}_{\mathrm{r}}$) can be found in many textbooks (e.g., Gill, 1982; Zdunkowski and Bott, 2003; Holton, 2004; Vallis, 2006). Evaluating them at the Equator (${\mathit{\phi}}^{(\prime )}=\mathrm{0}$) yields almost immediately the components of the governing equations for $\mathit{v}={v}_{\mathfrak{t}}{\mathit{e}}_{\mathfrak{t}}+{v}_{\mathfrak{n}}{\mathit{e}}_{\mathfrak{n}}+w{\mathit{e}}_{\mathrm{r}}$ in the local edge coordinate systems if we identify v_{𝔱} with u^{′} and v_{𝔫} with v^{′}. The equation for the vertical velocity component w in the local coordinate system defined in the center of the vertical cell interfaces can be found in a similar way. Together, they read
where ${K}_{\mathrm{h}}=({v}_{\mathfrak{t}}^{\mathrm{2}}+{v}_{\mathfrak{n}}^{\mathrm{2}})/\mathrm{2}$ is the horizontal massspecific kinetic energy, and f_{r}=2Ωsin (φ), ${f}_{\mathfrak{t},\mathfrak{n}}=\mathrm{2}\mathrm{\Omega}\mathrm{cos}\left(\mathit{\phi}\right){\mathit{e}}_{\mathit{\phi}}\cdot {\mathit{e}}_{\mathfrak{t},\mathfrak{n}}$ are the Coriolis parameters. The values of the projections ${\mathit{e}}_{\mathit{\phi}}\cdot {\mathit{e}}_{\mathfrak{t},\mathfrak{n}}$ are already provided by the standard shallowatmosphere configuration of ICON, so they pose no additional problem. Here and in the following, we will formulate the deepatmosphere equations as a modification of the shallowatmosphere equations. This simplifies the comparison and corresponds actually to the way we have implemented the deepatmosphere dynamics in ICON. For instance, we made use of the expansion of the gradient in the local coordinate systems $\mathbf{\nabla}={\mathit{e}}_{\mathfrak{t}}\mathit{\{}a/r\mathit{\}}\partial /\partial \mathfrak{t}+{\mathit{e}}_{\mathfrak{n}}\mathit{\{}a/r\mathit{\}}\partial /\partial \mathfrak{n}+{\mathit{e}}_{\mathrm{r}}\partial /\partial r$ if $\mathbf{\nabla}={\mathit{e}}_{\mathfrak{t}}\partial /\partial \mathfrak{t}+{\mathit{e}}_{\mathfrak{n}}\partial /\partial \mathfrak{n}+{\mathit{e}}_{z}\partial /\partial z$ is its expansion under shallowatmosphere approximation (with $\partial /\partial \mathfrak{t}$ and $\partial /\partial \mathfrak{n}$ denoting the horizontal derivatives along great circle arcs under shallowatmosphere approximation, i.e., $r=a+z\approx a$, and e_{z}=e_{r}, $\partial /\partial z=\partial /\partial r$). So the factors {⋅} become 1 under the shallowatmosphere approximation. In addition, the underlined terms in Eqs. (9) and (10) are neglected in the shallowatmosphere equations.
The velocity component v_{𝔱} is reconstructed diagnostically from the adjacent values of v_{𝔫} by a radial basis function interpolation and subsequent projection. More details on this approach can be found in Wan (2009, Appendix D).
The dissipative and diabatic processes summarized in the terms ρ^{−1}F and RQ∕(c_{v}c_{p}ρθ) in Eqs. (1) and (3) are not modified for the deep atmosphere. Their modification for spherical coordinates would have introduced numerous additional metric terms into the governing equations (compare, e.g., Baldauf and Brdar, 2016), whose rigorous implementation into the dynamical core would cause considerable additional computational costs, which are not affordable for us for the time being.
2.1.2 Numerical implementation
A thorough description of the spatial and temporal discretization of the governing equations used in ICON can be found in Zängl et al. (2015), and for development stages of ICON in Bonaventura and Ringler (2005), Wan (2009) and Wan et al. (2013). Here, we will focus only on those elements of the spatial discretization that are affected by the deepatmosphere modifications.
First, we would like to point to an important simplification, which we have done in our implementation of the deepatmosphere dynamics into ICON. Since ICON is used for operational NWP at DWD, a key criterion for new developments to be integrated into the dynamical core is computational efficiency (of course, this criterion applies to climate models as well). For this reason, priority is given to efficiency over accuracy, where we assume this to be acceptable. Where possible the deepatmosphere dynamics are realized by modification factors to the existing terms of the shallowatmosphere dynamics. In addition, topography is not taken into account by the deepatmosphere modification, which, for instance, means that the dynamics in a cell of the first grid layer above the Himalaya experience the same deepatmosphere modification as the dynamics in a cell above the ocean surface. This simplification, albeit the most severe one, has many advantages. First, we avoid the complicated calculations of the geometric measures of a cell (e.g., its volume and surface areas) in spherical geometry if they are distorted by topography and lose the center of Earth as their center of curvature. Second, the abovementioned deepatmosphere modification factors depend only on height, which saves a considerable amount of memory and computational cost. Since topography imprints on the grid layers only up to a certain height in ICON (typically up to 16 km), the errors introduced by this measure are assumed to be relatively small. For instance, the difference of a typical grid cell volume at mean sea level and at 16 km is of the order of a few per mille in the spherical geometry. Furthermore, the massconserving property of the dynamical core of ICON is not affected by our simplification.
We start with the modification of the gravitational acceleration (Wood and Staniforth, 2003):
with the modification factor parenthesized by braces. The formula (Eq. 11) enters Eq. (10) only implicitly via the hydrostatic background state. For its computation, Eq. (11) is used only once during model initialization. Next, we consider the gradient in edgenormal and tangential directions:
where χ stands for the massspecific horizontal kinetic energy K_{h} and the Exner pressure π in Eq. (9), and the vertical velocity w in Eq. (10). $\mathrm{\Delta}\mathit{\chi}/\mathrm{\Delta}\mathfrak{t},\mathfrak{n}$ denotes the gradient as it is computed under shallowatmosphere approximation. The factor a∕r derives from the modification of the length of horizontal distances $l\to l\mathit{\{}r/a\mathit{\}}$. Like all other modification factors, it is precomputed during model initialization for the vertical positions, where the gradients have to be evaluated during runtime. Another quantity the computation of which has to be modified for the deep atmosphere is the vertical vorticity component ζ on the lefthand side of Eq. (9). ζ is computed at the vertices of the triangular cells and then interpolated to the edge centers, where v_{𝔫} is defined. This allows to compute the vorticity from the v_{𝔫} using the Stokes theorem:
where A_{v} is the area of the hexagonal or pentagonal cells centered at the vertices, which form the dual cells to the primal triangular cells by connecting the mass points of the triangular cells around a vertex; see Fig. 2 in Wan et al. (2013). In addition, n_{e} is the number of the dual edges of length l_{d} around the vertex, crossing the primal edges of the triangular cells perpendicularly (n_{e}=6, or 5 for the 12 pentagon points of the grid, respectively). Finally, f_{o} is an orientation factor, which is 1 or −1 according to whether e_{𝔫} is parallel or antiparallel to the cyclonic direction of the integration path. Again, the first factor on the righthand side of Eq. (13) is the vorticity under shallowatmosphere approximation; thus, the geometric measures A_{v} and l_{d} are those at r=a. The modification factor a∕r results from the quotient l_{d}∕A_{v} in spherical geometry and is identical to the modification factor for the horizontal gradient (Eq. 12). Equations. (2) and (3) are written in flux form, which guarantees mass conservation up to machine precision and the conservation of heat as far as advective fluxes are concerned. The deepatmosphere modifications of the flux divergences result from the quotient of area and volume A∕V for the five surfaces of a triangular cell. The modification of a cell volume itself reads
where $V={A}_{\mathrm{c}}({z}_{\text{top}}{z}_{\text{bot}})$ is the cell volume under shallowatmosphere approximation, with the cell surface area A_{c}, and ${r}_{\text{bot}}=a+{z}_{\text{bot}}$, ${r}_{\text{top}}=a+{z}_{\text{top}}$ are the radii of the cell's bottom surface and top surface, respectively. This modification is required, for instance, if global integrals of mass and other quantities are computed for diagnostic purposes. Together with the area modification for the three side faces and the bottom and top surfaces,
where l_{p} is the length of the primal edges, we find for the flux divergences
Here, ${\mathit{\xi}}_{i,\text{top},\text{bot}}$ denotes the scalar product of the flux density and e_{𝔫}, or e_{r} at the respective cell face. The first term on the righthand side of Eq. (17) sums the fluxes across the side faces, where the orientation factor f_{o} is either 1 or −1 according to whether e_{𝔫} is parallel or antiparallel to the normal vector of the side faces from the point of view of the cell. The second term sums the fluxes across the top and bottom surfaces.
The underlined terms in Eqs. (9) and (10) are new implementations in the dynamical core. Their discrete formulation closely follows the formulation of structurally similar terms of the discrete shallowatmosphere equations in Zängl et al. (2015).
The discretization employed in ICON necessitates the spatial interpolation of variables between cell center, edge midpoint, vertex, and the center of the vertical cell interfaces (bottom and top surface; Zängl et al., 2015), which could be affected potentially when changing from a shallow to a deepatmosphere formulation. In the horizontal, linear, bilinear, and areaweighted interpolations are used primarily. These interpolations are assumed to take place on coordinate surfaces z=const., i.e., planes in the case of the shallow atmosphere (at least practically, but see Thuburn and White, 2013, for the geometric implications of the shallowatmosphere approximation) and spherical shells in the case of the deep atmosphere. This means that no explicit deepatmosphere modification is necessary, because it would be the same for the respective geometric measures (distances, areas) on a coordinate surface, which cancel each other in the interpolation. For the vertical interpolation between cell centers and the cell interfaces, a linear interpolation is used in the case of the shallow atmosphere. We interpret this as an interpolation along the coordinate lines $\mathit{\lambda},\mathit{\phi}=\text{const.}$ in the case of the deep atmosphere, so that we can avoid deepatmosphere modifications in this case, too. If we interpret the vertical interpolation, e.g., as a volumeweighted average, a modification would be necessary. Apart from the aforementioned standard interpolations, ICON makes use of an upwindbiased reconstruction of density and potential temperature from cell centers to the centers of the (horizontal and vertical) interfaces. These reconstructions enter, for instance, the divergence of mass and heat flux (Zängl et al., 2015) and make use of a Taylor expansion up to the first order along backward trajectories, which requires modifying the corresponding gradient according to Eq. (12).
For the time integration, a twotimelevel predictor–corrector scheme is used, which is not directly affected by deepatmosphere modifications. However, certain terms on the righthand side of the governing equations, which are associated with vertical sound wave propagation, are treated implicitly in the temporal discretization in order to allow reasonable time steps. This directly affects the discretized forms of Eqs. (3) and (10), which can be combined to yield a diagnostic, elliptic equation for the vertical wind w for each vertical grid cell column of the model atmosphere (Zängl et al., 2015). Thus, some of the coefficients of this equation have to be complemented by the corresponding deepatmosphere modification factors.
2.1.3 Model initialization
One motivation to implement the deepatmosphere dynamics in ICON is to increase the accuracy of simulations with a model top ≳100 km. However, initial data are usually only available up to altitudes of about 70 to 80 km. For instance, the model top is at 75 km in the operational NWP with ICON at DWD, and the IFS model of the ECMWF, whose operational analysis data can be used to initialize ICON, has its model top at 0.01 hPa (see https://www.ecmwf.int/en/forecasts/documentationandsupport, last access: October 2018). To our knowledge, no regular and reliable measurements with global coverage are available for the (lower) thermosphere, which could be used for some kind of data assimilation. Therefore, climatological tables appear as a possible second choice. However, the momentary state at a point in the atmosphere generally deviates more or less strongly from the state given by any climatology. The model thermosphere initialized with climatological data will consequently undergo an adjustment process during some spinup phase. With this uncertainty in mind, we decided to begin with a very simple approach. The state of the lower and middle atmosphere is initialized by an IFS analysis. In order to complement this by a state for the upper atmosphere, a mean vertical temperature profile (neglecting horizontal and temporal variations) obtained from Bates (1959), Hedin (1983), and Fleming et al. (1988) is used. It provides the basis to obtain first the pressure from the hydrostatic balance, followed by a computation of the horizontal wind from the geostrophic balance. The technical details can be found in Appendix A. Later on, this simple approach could be improved on by using a climatology generated by UAICON itself.
2.2 Upperatmosphere physics
A new physics package, which parameterizes processes specific to the upper atmosphere has been developed for UAICON. This package, referred to as the UA package, can be called in combination with either the NWP package or the ECHAM package. The processes taken into consideration in the UA package are summarized in Table 1.
Huang et al. (1998)Banks and Kockarts (1973)Gill (1982)Hong and Lindzen (1976)Hines (1997a, 1997b)Strobel (1978)Richards et al. (1994)Fomichev and Blanchet (1995, 1998)Ogibalov and Fomichev (2003)Kockarts (1980)The processes are categorized into three groups (kinetics, radiation, and chemical heating), as described below. Most of the parameterizations are adopted from the HAMMONIA model, a spectral model based on ECHAM5 (Roeckner et al., 2006), covering the atmosphere up to the thermosphere ($\mathrm{1.7}\times {\mathrm{10}}^{\mathrm{7}}$ hPa, ∼250 km). A detailed description of the physics parameterizations used in HAMMONIA can be found in Schmidt et al. (2006); thus, here we keep the description brief, only noting important and differing treatments. An overall difference is that in UAICON all the parameterizations are implemented such that the computation only starts at a certain altitude above which the forcings are expected to become relevant (see Table 1). This increases computational efficiency significantly.
2.2.1 Kinetics
Above the mesopause, molecular diffusion, which is negligible at lower altitude, becomes significant. In fact, there the downward transport of heat by molecular diffusion appears as a strong cooling in balance with the strong solar heating. Hence, it is of primary importance to parameterize molecular processes in this region of the upper atmosphere. Molecular transport of heat, momentum, and tracers is parameterized in the UA package following Huang et al. (1998) and Banks and Kockarts (1973), as in HAMMONIA. The computation starts at 75 km. Besides direct transport of heat by molecular diffusion, the momentum transport also leads to energy deposition in the form of heat, known as frictional heating. In the UA package, this is parameterized following Gill (1982). The computation starts at 85 km.
In the mesosphere and lower thermosphere, unlike in lower layers, a larger number of air particles are ionized and thus aligned with the electromagnetic field of the Earth. This produces a force on the neutral mass flow, its tangential and normal components known as the ion drag and the Lorenz force, respectively. Joule heating is produced by the ion drag as well. In UAICON, as in HAMMONIA, this effect is parameterized following the simple Hong and Lindzen (1976) approach. The computation starts at 80 km.
A large portion of gravity wave momentum energy is deposited and transferred to turbulent energy near the mesopause, where turbulence is otherwise very weak. The turbulent mixing effects induced by gravity waves can be estimated using the Hines (1997a, 1997b) parameterization included in the ECHAM package. This option is switched off in standard ICON simulations. In UAICON simulations, however, we enable the calculation and pass the computed turbulent diffusion coefficient to the turbulent mixing subroutine to account for gravitywaveinduced turbulent mixing.
2.2.2 Radiation
In standard ICON, the PSrad radiation package (Pincus and Stevens, 2013) is employed, which is itself an extension to the Rapid Radiative Transfer Model for GCMs (RRTMG) (Mlawer et al., 1997; Iacono et al., 2008). The shortwave component of the radiation package covers wavelengths of the solar spectrum longer than 200 nm. This is a sufficient bandwidth only up to about the mesopause, beyond which radiative heating in the ultraviolet and extreme ultraviolet frequencies become dominant. In the UA package, starting from 50 km, ultraviolet solar forcing for the O_{2} Schumann–Runge bands (SRBs; 175 to 205 nm) and continuum (SRC; 125 to 175 nm) is calculated based on the model of Strobel (1978). Efficiency factors multiplied by the SRBC (SRB plus SRC) heating rates account for the loss of internal energy due to airglow processes and are taken from Mlynczak and Solomon (1993). For the extreme ultraviolet (EUV; 5 to 105 nm) solar forcing, starting above 90 km, a model based on Richards et al. (1994) taken from HAMMONIA is currently used. Efficiency factors multiplied by the EUV heating rates are based on Roble (1995). Their values also account for the energy loss due to radiative cooling in the 5.3 µm NO band. Since this process is explicitly calculated in our model (see below), a factor of 1.33 is multiplied by these efficiency factors (see Richards et al., 1982). Some additional adjustments to the PSrad/RRTMG shortwave radiation are necessary in UAICON due to the introduction of chemical heating. More details on this are given in Sect. 2.2.3.
The longwave component of PSrad/RRTMG covers terrestrial wavelengths shorter than 1 mm. This bandwidth is still valid at thermospheric heights, yet a few important additions have to be made:

The usual assumption of local thermodynamical equilibrium (LTE) does not hold above the mesopause; thus, nonLTE effects must be taken into account. As in HAMMONIA, nonLTE infrared cooling by O_{3} and CO_{2} is calculated from the parameterization of Fomichev and Blanchet (1995) with the modifications of Fomichev and Blanchet (1998). The calculation starts at 65 km, and the calculated values are multiplied by a scaling factor α equaling 0 at 65 km and linearly growing to 1 at 75 km. Correspondingly, the longwave radiation computed by PSrad/RRTMG is scaled with the factor 1α, effectively discarding it above 75 km.

As in HAMMONIA, a parameterization of CO_{2} nonLTE absorption in the nearinfrared following Ogibalov and Fomichev (2003) is employed. The computed values are ignored below 19.25 km and fully considered above 24.5 km.

NO cooling at 5.3 µm is calculated utilizing the parameterization from Kockarts (1980). The computation starts at 60 km.
A noteworthy difference between ICON and HAMMONIA is that ICON is a nonhydrostatic model on height levels, whereas HAMMONIA is hydrostatic on hybrid pressure levels. In HAMMONIA, for the use in the radiation computation, number densities of radiatively active tracers are calculated based on the mass of air in a given layer which is derived from pressure differences between the upper and lower surfaces of a layer. This approach is only valid under the assumption of hydrostatic balance, since otherwise pressure is not guaranteed to decrease strictly monotonically with increasing altitude. Therefore, in the UA package, the computation of number density for the radiation parameterization is utilizing the mass of air, which is a globally conserved quantity in ICON.
Moreover, HAMMONIA has the upper boundary of its top pressure level at 0 hPa, effectively covering the whole atmosphere. The height levels of ICON, on the other hand, cover a finite range and unavoidably ignore the atmospheric air mass above the model lid. The effect of this missing amount of air on radiative fluxes is ignored in UAICON.
2.2.3 Chemical heating
Chemical heating, i.e., the release of heat due to recombination reactions between atoms or radicals produced as a result of photolysis, becomes important in maintaining the upperatmospheric thermodynamic balance, where the photolysis products can travel a long distance before recombining. HAMMONIA employs a condensed version of the MOZART3 chemistry model (Kinnison et al., 2007) to explicitly compute chemical heating online. In our case, however, given the finer target resolution (from 160 km down to a few tens of kilometers) and the central goal of studying gravity waves, using a coupled chemistry model is overly expensive. Therefore, we deploy a simpler strategy of prescribing monthly zonalmean climatological chemical heating rates from a 35year HAMMONIA simulation with constant presentday boundary conditions. As an example, the chemical heating rates for January are shown in Fig. 2. Technically, below 70 km, all heating is calculated in the PSrad/RRTMG shortwave radiation code and no chemical heating rates are prescribed, whereas above 80 km full chemical heating rates are used and the radiative heating provided by the PSrad/RRTMG scheme is reduced by 23 % in order not to count solar energy twice. This approximation accounts for the energy used to break the chemical bond of ozone. Additionally, heightdependent efficiency factors from Mlynczak and Solomon (1993) are applied at these altitudes. This approach was also used in HAMMONIA (Schmidt et al., 2006). Between 70 and 80 km, the two heating sources are linearly merged.^{2}
3.1 Idealized test cases
To test the deepatmosphere implementation in the dynamical core, we used two test cases. In the first test case, the propagation of a sound wave is considered, for which an analytical solution (of the linearized equations) is available. It is aimed at testing especially the accuracy of the spherical geometry in its imprint on the grid cells and the corresponding modification factors described in Sect. 2.1.2, and at testing the metric terms and the complete Coriolis acceleration in the components of the momentum Eqs. (9) and (10). The second test case is the Jablonowski–Williamson baroclinic instability test case (Jablonowski and Williamson, 2006) in its extension for deepatmosphere dynamical cores by Ullrich et al. (2014). It reveals if the height dependence of gravity is properly implemented and maintains the hydrostatic background state of the test case atmosphere (especially if gravity enters the momentum equation only implicitly, as in Eq. 10). In addition, the performance of the entire dynamical core is tested when it comes to reproducing the development of the baroclinic wave. Both test cases make use of the smallEarth approach of Wedi and Smolarkiewicz (2009) to pronounce deepatmosphere effects.
3.1.1 Sound wave test case
The particular motivation for this first test case is that an analytical solution is available to which the numerical solution can be compared. We have developed this test case with a method originally proposed by Läuter et al. (2005) for the shallowwater equations on the sphere. The method was developed further for the shallow and deepatmosphere equations, e.g., by Staniforth and White (2008); Baldauf et al. (2014). An atmosphere at rest in the absolute frame is considered. If a nontrivial analytical solution is known for this case, it can be transformed into a rotating frame (e.g., regarded as a rotating Earth slipping through the air without exchange of tangential momentum). Depending on the solution in the absolute frame being either stationary or time dependent, potentially all aspects of a dynamical core can be tested. However, a disadvantage of this method is that the centrifugal acceleration has to be taken into account explicitly in the dynamical core (see Eq. 1). Some aspects of this transformation method are shown in Appendix B1, and a thorough mathematical description can be found in the literature cited above.
Given a solution in the absolute frame, it appears to be advected with $\mathit{v}={\mathit{v}}_{F}=\mathbf{\Omega}\times (\mathit{X}\mathit{A})$ from the perspective of the rotating frame (with the center of Earth A and an arbitrary point X not coincident with A). In practice, this means that the solution has merely to be rotated by an angle −Ωt about an axis being parallel to Ω and crossing A. Therefore, we will direct our attention to the solution to be rotated in the following.
Baldauf et al. (2014) derived analytical gravity and sound wave solutions for the linearized deepatmosphere equations. However, certain terms of the equations had to be omitted in order to allow the solution to be expanded in a system of orthonormal basis functions. Although these terms were shown to be only of second order, their effects are not controllable in a dynamical core, where the omission of most of the corresponding discretized terms is unfeasible without greater effort. Therefore, it would be desirable to find solutions for the linearized equations with the omitted terms restored, or alternatively with only those omissions retained which can be realized in a dynamical core without greater effort. The fact that the gravity $g(a/r{)}^{\mathrm{2}}(\mathit{X}\mathit{A})/r$ imprints a spherical symmetry on the equations with its distinct point A turns out to be a severe obstacle to that aim. So we decided to switch off gravity (which in model practice means to set the constant parameter g to a very small value >0, since divisions by g are used throughout the model, especially in the physics parameterizations). This greatly simplifies the problem but has the consequence that the test case makes no statement about the implementation of gravity. Under these circumstances, the atmosphere, enclosed between the spherical boundaries of the model bottom and top, is isotropic, with the constant pressure p_{0} and temperature T_{0}. As long as the sound waves, propagating with the speed of sound ${c}_{\mathrm{s}}=\sqrt{({c}_{p}/{c}_{v})\mathrm{R}{T}_{\mathrm{0}}}$, do not interact with the boundaries, they are not “aware” of the spherical shape of the atmosphere as a whole. Therefore, the challenge for the model is to properly simulate the sound wave propagation on the anisotropic spherically curved grid. For this test case, we consider a spherically symmetric acoustic wave, shown schematically in Fig. 3, which consists of an outwardpropagating part, only. The derivation is shown in Appendix B2, and the solution for the pressure perturbation ${p}^{\prime}=({c}_{p}/\mathrm{R})({p}_{\mathrm{0}}/{\mathit{\pi}}_{\mathrm{0}}){\mathit{\pi}}^{\prime}$ associated with the sound wave reads
where $\mathit{\delta}p=({c}_{p}/\mathrm{R})(\mathit{\delta}T/{T}_{\mathrm{0}}){p}_{\mathrm{0}}$ denotes the pressure amplitude of the wave, determined by a temperature amplitude δT in our implementation, $x=\mathit{X}\mathit{B}$ is the distance from the center of the spherical wave B, $x={b}_{\mathrm{1}}>\mathrm{0}$, and $x={b}_{\mathrm{2}}>{b}_{\mathrm{1}}$ are the radial boundaries within which the wave has a nonvanishing amplitude at t=0 (see Fig. 3), and n is the number of wave crests (in this work, we consider only n=1). In addition, Θ denotes the Heaviside step function defined as (e.g., Bronstein et al., 2001)
The solution (Eq. 18) is only valid until the first reflection occurs. Of course, sound waves have been thoroughly investigated in the literature and solutions to the sound wave equation are all but new (e.g., Kirchhoff, 1876), but since their propagation in combination with the method of Läuter et al. (2005) involves potentially all parts of a dynamical core (except for gravity), we found this test useful, also in view of the relative scarcity of test cases dedicated to deepatmosphere dynamical cores in the literature. In order to highlight the effect of the spherical curvature (on the model grid), the radius of Earth can be rescaled a→η_{1}a, (η_{1}<1). This is the smallEarth approach proposed by Wedi and Smolarkiewicz (2009) (see also Baldauf et al., 2014; Ullrich et al., 2014). The model time step is rescaled accordingly in order to account for the correspondingly smaller mesh size of the horizontal grid. Furthermore, it might be advantageous to rescale the angular velocity of the Earth, Ω→η_{2}Ω, in order to control the velocity $\mathit{v}={\mathit{v}}_{F}$ with which the sound wave is advected. Further details on the implementation can be found in Appendix B2. Apart from that, we followed closely the guidelines given by Baldauf et al. (2014) for the implementation into UAICON.
We envisaged two concrete test configurations: the first without rotation (η_{2}=0), to simulate the sound wave propagation in the absolute frame, and the second with rotation (η_{2}>0), to test if the dynamical core is able to maintain the balance of the background velocity in the rotating frame $\mathrm{d}{\mathit{v}}_{F}/\mathrm{d}t+\mathrm{2}\mathbf{\Omega}\times {\mathit{v}}_{F}+\mathbf{\Omega}\times [\mathbf{\Omega}\times (\mathit{X}\mathit{A}\left)\right]=\mathrm{0}$, so as to advect the sound wave in a shapeconserving way. Further parameter settings are listed in Table 2. The temperature amplitude of the sound wave δT was chosen small enough that its nonlinear dynamics as computed by the dynamical core is negligibly small (an initial amplitude of δT=0.1 K appears large when compared, e.g., to typical values used in Baldauf et al. (2014), but even for larger amplitudes, ∼1 K, the numerical solution of UAICON was in relatively good agreement with the analytical solution of the linearized equations).
A height–longitude cross section at the Equator of the numerical solution from UAICON, the analytical solution, as well as the difference between the two for both configurations are shown in Fig. 4, shortly before the periphery of the sound wave would impinge on the bottom and top boundaries. The angular velocity of the second configuration was chosen such that the center of the sound wave, B, would be advected by the zonal background wind $u\left(\mathit{B}\right)={u}_{F}\left(\mathit{B}\right)=\mathrm{\Omega}r\mathrm{cos}\left(\mathit{\phi}\right){}_{\mathit{B}}=\mathrm{100}$ m s^{−1} in the rotating frame. This is close to the value used by Baldauf et al. (2014) in their test scenario (C). We use a time step of $\mathrm{\Delta}t={\mathit{\eta}}_{\mathrm{1}}\cdot \mathrm{13.2}$ s for both configurations. It should satisfy the Courant–Friedrichs–Lewy (CFL) criterion in both configurations. The maximum propagation velocity in the first configuration is $v{}_{max}={c}_{\mathrm{s}}=\mathrm{317}$ m s^{−1}, whereas in the second configuration it is $v{}_{max}={u}_{F}{}_{max}+{c}_{\mathrm{s}}=\mathrm{134}$ m s^{−1} +317 m s^{−1} =451 m s^{−1}.
Shape and amplitude are relatively well captured by the numerical solution in both configurations, with the difference in amplitude to the analytical solution being about 1 order of magnitude smaller than the magnitude of the wave's pressure perturbation itself. However, in the second configuration, where the sound wave is advected westward while radially propagating, the magnitude of the error has increased slightly, and the symmetry of the pressure difference with respect to a vertical axis crossing the center of the sound wave is lost due to the horizontal advection. The amount by which symmetry is lost is a measure for the phase error of the horizontal advection implementation (e.g., Skamarock and Klemp, 2008). The pressure difference in the first configuration not being radially symmetric with respect to the center of the sound wave is probably due to at least three anisotropies between the vertical and the horizontal: first, the horizontal and vertical mesh sizes (Δx=597.8 m and Δz=555.6 m) slightly differ. Second, the extension of a grid cell increases with height in the spherical geometry, and third, the horizontally explicit–vertically implicit scheme employed in the dynamical core of ICON introduces an anisotropy as well.
We repeated the simulation for different grid resolutions and computed the L^{2} norm and L^{∞} norm of the pressure difference between the numerical and analytical solutions on the entire circumequatorial height–longitude cross section, of which a part is plotted in Fig. 4, and at time t=60 s, according or in analogy to the formula employed by Baldauf and Brdar (2014, p. 1983). All pressure values entering the computation of the two norms are weighted equally, i.e., no weighting with the cell volume is applied (in which case the two norms would not be with respect to the pressure difference Δp but with respect to a worklike quantity ∝ΔpV). The results and some further information on the employed grids are shown in Fig. 5. In the first configuration, without rotation, the convergence rate is dominated by a secondorder behavior, although a relatively small firstorder component seems to be present, especially in the case of the L^{∞} norm. In the second configuration, with rotation, the convergence rate seems to start with a secondorder behavior for the lower grid resolutions and changes into a firstorder behavior for the higher resolutions. This is in agreement with the findings of Baldauf et al. (2014) for their test scenario (C) (compare their Fig. 7). The reason for the firstorder convergence in the presence of a background wind is still unknown. Nevertheless, we regard the agreement between the analytical and numerical solutions in both configurations as satisfactory, as we assume a critical deficiency in the deepatmosphere modification of the dynamical core to leave a much more distinct fingerprint in the numerical solution.
3.1.2 Jablonowski–Williamson baroclinic instability test case
The previous test case focused on one particular emergent structure of the atmosphere. However, if we turn our focus to the atmospheric features on the synoptic scale, other structures, such as baroclinic waves, are much more important than sound waves. The Jablonowski–Williamson baroclinic instability test case (Jablonowski and Williamson, 2006) is a standard test to investigate the performance of atmospheric models in representing a key feature of midlatitude dynamics. It consists of a baroclinically unstable atmosphere in hydrostatic and geostrophic balance to which a perturbation is added which triggers the instability. This test case reveals, on the one hand, if the model is able to maintain the hydrostatically and geostrophically balanced background state during the first days of the wave evolution, when its amplitude is still relatively small, and on the other hand, how the model performs in reproducing the amplitude growth of the wave and its shape. However, a disadvantage of this test is that no analytical solution for the problem is known, so the evaluation has to be based on a model intercomparison. Ullrich et al. (2014) have extended this test case for deepatmosphere dynamical cores and introduced some further improvements to the original formulation. The approach of the small Earth is employed to highlight the differences between the shallow and the deepatmosphere dynamics. The rescale factors are ${\mathit{\eta}}_{\mathrm{1}}=\mathrm{1}/\mathrm{20}$ and η_{2}=20, for the Earth radius and angular velocity, respectively. For the test with UAICON, we used a R3B4 grid which provides a horizontal mesh size of $\mathrm{\Delta}\mathit{\phi}={\mathrm{0.95}}^{\circ}$. This is close to the value used for the production of the numerical benchmark solution in Ullrich et al. (2014). The vertical grid is stretched, with layer thicknesses increasing from the model bottom to the model top at 30 km. Following Ullrich et al. (2014), we use n_{lev}=30 levels; however, the vertical stretching of the ICON grid differs from their formula (28):
where z_{j} denotes the height of the jth interface separating the layer j−1 from the layer j. The index j counts from the model top to the model bottom, with z_{1}=H_{top} and ${z}_{{n}_{\text{lev}}+\mathrm{1}}=\mathrm{0}$. With a value of μ=15 for the flattening parameter, the lowermost and uppermost layers have a thickness of ${z}_{{n}_{\text{lev}}}\approx \mathrm{82}$ m and ${H}_{\text{top}}{z}_{\mathrm{2}}\approx \mathrm{1249}$ m, respectively. For ICON, the formula reads
A value of 1 is used for the stretching parameter σ, and a value of 3.16 follows for the exponent λ from Eq. (21) if we require a thickness of ${z}_{{n}_{\text{lev}}}=\mathrm{100}$ m for the lowermost layer. This setting results in a thickness of ${H}_{\text{top}}{z}_{\mathrm{2}}\approx \mathrm{1969}$ m for the uppermost layer. We assume that the differences between the vertical grids (20) and (21) are negligible, since tests with n_{lev}=60 and n_{lev}=120 revealed that the numerical solution is largely converged on the vertical grid of ICON for n_{lev}=30 (not shown). The results for days 8 and 10 of the simulation with UAICON are shown in Fig. 6a, b. In order to study the convergence behavior with regard to the horizontal grid resolution, we doubled the same twice (see Fig. 6c–f). First of all, we can state that UAICON is able to maintain the hydrostatic and geostrophic balances of the background state in the first days of the simulation relatively well. This indicates, for instance, that the vertical variation of the gravitational acceleration is adequately implemented. Second, the amplitude and shape of the baroclinic wave, as they show up in the surface pressure in Fig. 6, compare relatively well to the benchmark solution of Ullrich et al. (2014, their Fig. 9), and also to the solution of Wood et al. (2014, their Figs. 4 and 5). However, some differences can be recognized, especially in the tail of the baroclinic wave. The convergence tests revealed that the numerical solution is largely converged with regard to the horizontal resolution in the zonal range from 120 to 240^{∘} E (120^{∘} W), say. As mentioned before, the tail of the baroclinic wave in the zonal range from about 60 to 120^{∘} E shows a greater variation between the different horizontal resolutions (Fig. 6). To see if the resolution R3B6 is converged in that regard, we tested the R3B7 grid as well. However, it developed a numerical instability in the tail region of the baroclinic wave around day 8. The reason for the instability has not yet been clearly identified, but we assume that the nontraditional part of the Coriolis acceleration (see Eqs. 9 and 10) may play a role in its development, as tests, in which this part was switched off, indicate. This might result from the fact that our discretization does not satisfy $[\mathit{v}\cdot (\mathrm{2}\mathbf{\Omega}\times \mathit{v}){]}_{\text{discretized formulation}}=\mathrm{0}$; i.e., it violates momentum and energy conservation. A formulation of the Coriolis acceleration that satisfies these conservation principles is not possible within the discrete formulation of Eq. (1) used in ICON for reasons discussed by Gassmann (2011, 2018). In simulations of the standard Earth, the violation of the conservation principles is typically so small that it does not pose a problem. In contrast, the smallEarth approach seems to amplify the violation to such an extent that numerical instabilities may be triggered under certain conditions (e.g., a relatively high horizontal grid resolution in this case). However, a clear statement about this is difficult, since the stationary analytical atmospheric background state prescribed for this test case follows from balances that include the nontraditional part of the Coriolis acceleration. As a consequence, the model atmosphere undergoes an adjustment process right away from the beginning of the simulation if this part is switched off. This makes it difficult to evaluate if the instability is absent in the considered time period because of the disabled Coriolis acceleration or because of the atmospheric state being changed by the adjustment process. Nevertheless, given that the instability is absent in standard Earth simulations and measures against it could at most treat the symptoms, not the cause (leaving aside an extensive reformulation of the dynamical core that satisfies the aforementioned conservation principles), we decided to postpone further investigations. Apart from that, the comparison of Fig. 6 with the benchmark from Ullrich et al. (2014) makes us confident that our deepatmosphere implementation is satisfactory for our purposes.
3.2 Climatological test cases
3.2.1 Simulation setup
For the evaluation of the model climatology, a UAICON simulation with the upperatmosphere physics coupled to the ECHAM physics package has been performed. The deepatmosphere dynamics is also switched on. The model was integrated for 20 years with climatological boundary conditions: sea surface temperature and sea ice concentration are averaged for each calendar month from the Program for Climate Model Diagnosis and Intercomparison (PCMDI) Atmospheric Model Intercomparison Project (AMIP) dataset (Taylor et al., 1998) version 1.1.2 over 1979–2014; concentrations of radiatively active gases, namely O_{3}, CO_{2}, O_{2}, O, NO, CH_{4}, and N_{2}O, are averaged in the same manner from a 35year HAMMONIA simulation with fixed presentday boundary conditions; concentrations of CFC11 and CFC12 are fixed at 214.5 and 371.1 pptv, respectively; the 1865 condition of the tropospheric background aerosol from the MACv1 dataset (Kinne et al., 2013) is used; no volcanic or anthropogenic aerosols are used; landsurface parameters for the parameterization of the effects of subgridscale orography and for the embedded version of the Jena Scheme for Biosphere Atmosphere Coupling in Hamburg (JSBACH) landsurface model (v4; Giorgetta et al., 2018) are fixed as described by Giorgetta et al. (2018). The total solar irradiance is held constant at 1361.371 W m^{−2}, and the F10.7 index for the calculation of EUV heating rates is fixed at 150 sfu ($\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{sfu}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{22}}$ $\mathrm{W}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{Hz}}^{\mathrm{1}}$).
The simulation uses the R2B4 grid, which has a horizontal mesh size of about 160 km. In the vertical, the model uses 120 layers for the altitude range from the surface up to 150 km. Rayleigh damping (Klemp et al., 2008; Zängl et al., 2015) is applied above 120 km with a maximum damping coefficient of 10 s^{−1} at the top. Such strong damping is necessary to allow for a reasonable computational time step despite the occasionally very large vertical velocities in the thermosphere. The model was integrated with a (physical) time step of 4 min and five dynamical substeps each physical time step. Radiation parameterizations – i.e., the PSrad radiation scheme of ICON, the shortwave radiation in the SRBC and EUV, the nonLTE longwave radiation, the NO heating – are evaluated once every hour, whereas all other parameterizations are evaluated every time step. For all physics parameterizations, we apply the “allfast” treatment described by Giorgetta et al. (2018). For nonorographic gravity wave drag, a cutoff maximum vertical wavelength of 12 km is applied, thus prohibiting long gravity waves. Disabling these long gravity waves is physically sensible, as they are believed to strongly propagate horizontally and be subject to internal reflection before reaching the mesopause (Hines, 1997b).
Companion simulations have been performed with two ICON configurations using a standard model lid at 80 km, Rayleigh damping (maximum damping coefficient 1 s^{−1}) applied above 50 km, and 100 vertical levels exactly following the lower part of the vertical grid applied in the UAICON simulations. In the first configuration (referred to as ICON in the following), the deepatmosphere dynamics and the upperatmosphere physics are disabled. All other numerical and physical settings are identical to the UAICON run. The second configuration (referred to as ICON(UA)) additionally has the deepatmosphere dynamics and upperatmosphere physics enabled. With the help of these two configurations, we can estimate which of the differences between ICON and UAICON are due to the vertical extension and which are related to the application of extended physics and dynamics also below 80 km. The most important differences between these test cases are listed in Table 3.
Finally, a third configuration, denoted UAphysICON, has been simulated, which differs from the UAICON setup only in the deepatmosphere modification of the dynamics being switched off (see Table 3). We regard the comparison of UAICON and UAphysICON as a possibility to quantify the difference between the shallowatmosphere and the deepatmosphere dynamics.
When comparing the experiments ICON and ICON(UA), with their model top at 80 km, to UAICON, with its model top at 150 km, one has to take into account that the sponge layer affects the dynamics in ICON and ICON(UA) at an altitude range where there is no such impact in the UAICON configuration. As demonstrated by Shepherd et al. (1996), the distortion of the model dynamics by a sponge layer can extend even significantly below this layer by about two density scale heights. We have to accept this, since the use of a sponge layer is necessary to alleviate the adverse effects of a rigid model lid, such as wave reflection. Nevertheless, we assume that the model physics and dynamics of interest dominate to the extent that at least qualitative conclusions can be drawn from our comparison.
For the evaluation of the simulation results, a 15year (2002 to 2016) temperature climatology from observations of the Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) instrument on NASA's Thermosphere Ionosphere Mesosphere Energetics Dynamics (TIMED) satellite is used (v 2.0; Dawkins et al., 2018). A monthly zonalmean zonal wind climatology is taken from the Upper Atmosphere Research Satellite Reference Atmosphere Project (URAP; Swinbank and Ortland, 2003).
3.2.2 Comparison of simulated and observed climatologies
Figure 7 shows multiyear zonalmean temperatures for January and July from the UAICON and ICON simulations and from SABER. The observed temperature patterns are reasonably reproduced in the simulations for large parts of the stratosphere and mesosphere. UAICON simulates the low summer mesopause temperatures and its altitude well. The same is true for the weaker temperature minimum over the respective winter pole (the altitude range enclosed by the 200 K isoline near 100 km, at about 80^{∘} N in Fig. 7a and c, and at about 80^{∘} S in Fig. 7d and f). UAICON and ICON both simulate the stratopause and tropopause fairly accurately, although the winter stratosphere is slightly warmer than suggested by the observations. In addition, the tropical tropopause layer shows an increased thickness in the simulations as compared to the observations (compare, e.g., the equatorial region with temperatures between 180 and 200 K at an altitude of about 15 km between Fig. 7d and f). This may at least partly be related to the absence of stratospheric aerosol or to the prescribed climatological ozone.
Zonalmean zonal wind climatologies are presented in Fig. 8. Patterns simulated by UAICON and ICON agree qualitatively with the observationbased URAP climatology. The sign reversals of zonal wind in both hemispheres near the mesopause are simulated in UAICON but are in general too strong; i.e., the lowerthermospheric jets are too strong and peak at toolow altitudes. We are currently investigating if this can be adjusted by tuning the nonorographic gravity wave parameterization (without losing too much performance in lower parts of the atmosphere). Concerning stratospheric winds, the winter westerlies are too weak in ICON and UAICON, an issue already mentioned in the ICON evaluation by Crueger et al. (2018). While UAICON and ICON show very similar biases for the boreal winter jet, the problem is reduced in austral winter. It is no surprise that our UAICON simulation, performed with the same settings for the subgridscale orography (SSO) parameterization as used by Crueger et al. (2018), shows similar issues. A reduction of the orographic gravity wave sources would reduce this issue in particular in the Northern Hemisphere but has not been implemented by Crueger et al. (2018) as it would deteriorate nearsurface winds. Retuning of orographic and nonorographic gravity wave parameters is planned for future model versions.
In order to quantify how the vertical model extension, the upperatmosphere physics, and the deepatmosphere dynamics affect the state of the model atmosphere on climatological scales, we examine the zonalmean temperature difference as one possible measure for this purpose. Figure 9 shows the differences between UAICON and ICON, on the one hand, and between ICON(UA) and ICON on the other hand, again for the months of January and July. The differences below about 60 km are very similar in the left and the right panels, indicating that in this region they are mostly related to the extension of the dynamical and physical processes and not to the vertical extension. In most areas, the process extension leads to higher temperatures with the strongest signals of up to about 5 K in the summer middle stratosphere and even stronger in the winter lower mesosphere. Above about 60 km, the patterns of the differences are again similar but the magnitude is stronger for the difference between UAICON and ICON, meaning that here also the vertical extension adds a warming in comparison to the standard configuration of ICON. At the uppermost level of comparison, the temperature differences reach several tens of Kelvin. Despite the sponge layer gradually increasing in magnitude above 50 km in the two configurations, ICON and ICON(UA), the comparisons suggest that a vertical model extension beyond 80 km as implemented in UAICON even influences simulated climatological means down to at least 60 km.
The difference between UAphysICON and UAICON (i.e., the difference between shallowatmosphere and deepatmosphere dynamics) is shown in Fig. 10. The application of the deepatmosphere dynamics results in significantly lower temperatures above about 90 km as compared to the shallowatmosphere dynamics. This difference increases with height up to several tens of Kelvin at the model top and shows only relatively moderate meridional and temporal variation. A possible explanation for this might read as follows: in UAICON, the magnitude of the gravitational acceleration decreases with height (see Eq. 11), whereas it is constant in UAphysICON. So, at a given height, we would expect a larger air density on average in UAICON than in UAphysICON. In addition, the heating rates at a certain altitude that result from the parameterizations of radiative processes depend on the optical thickness of the air column above (for shortwave radiation), which in turn depends more or less directly on the air mass in the column above. That is, if we would use the air mass below the model top in a grid cell column as a vertical coordinate, and compare the mean temperatures from UAphysICON with UAICON in this coordinate, we would expect significantly smaller differences as compared to the comparison in the geometric height coordinate. To test this hypothesis, we computed the global horizontal mean of the temperature, denoted $\overline{T}\left(z\right)$, and the mean of the air mass in a grid cell column between the height z and the model top at height z_{top}, which we denote $\overline{m}\left(z\right)$. This was done for two simulation variants that are identical to the UAICON and UAphysICON simulations, except for using a shorter period. Each simulation variant was initialized with an IFS analysis for 1 November and integrated for 3 months, of which the data for January went into our analysis. We did this for the six winters from 2013/2014 to 2018/2019 and averaged over the January months. In addition, we combined the UA physics with the NWP physics instead of the ECHAM physics, given the relatively short integration period. For a more convenient comparison, we translated $\overline{m}$ into the “mass height” ${z}_{m}=H\mathrm{ln}[\overline{m}/{\overline{m}}_{\text{tot}}+(\mathrm{1}\overline{m}/{\overline{m}}_{\text{tot}}\left)\mathrm{exp}\right({z}_{\text{top}}/H\left)\right]$, which results from the hydrostatic balance (for the shallowatmosphere dynamics) $\partial \mathit{\rho}/\partial z=\mathit{\rho}/H$. Here, ${\overline{m}}_{\text{tot}}$ is the global mean total air mass in a grid cell column and we used a constant scale height of $H=\mathrm{R}T/g=\mathrm{7.6}$ km as a representative value for the observed range of $\overline{T}$. The results for $\overline{T}\left(z\right)$ and $\overline{T}\left({z}_{m}\right)$ are shown in Fig. 11. As expected, the difference $\left\mathit{\right\{}\overline{T}\left({z}_{m}\right){\mathit{\}}}_{\text{SA}}\mathit{\left\{}\overline{T}\right({z}_{m}\left){\mathit{\}}}_{\text{DA}}\right$ is significantly smaller than $\left\mathit{\right\{}\overline{T}\left(z\right){\mathit{\}}}_{\text{SA}}\mathit{\left\{}\overline{T}\right(z\left){\mathit{\}}}_{\text{DA}}\right$ above the mesopause region, in the geometric altitude range between about 95 km and about 130 km (between about 110 km and about 130 km in terms of the mass height). Here, “SA” and “DA” stand for the shallow atmosphere and deep atmosphere (i.e., UAphysICON and UAICON), respectively. Below and above this range, at least one of our assumptions seems to hold no longer. Below 95 km, in the mesosphere, the interplay of the radiative processes with the other physics may change. Above about 130 km, the proximity of the model top might be one reason why the hypothesis becomes invalid. As the air mass above the model top is neglected, the optical thickness in zenith direction approaches zero there. In addition to the comparison of $\overline{T}\left(z\right)$ and $\overline{T}\left({z}_{m}\right)$, we inspected the heating rates from the shortwave radiative processes within the geometric altitude range between 95 and 130 km at the beginning of the simulations, when the model atmosphere is still in its hydrostatically balanced initial state. They are systematically stronger in the UAphysICON setup than in the UAICON setup (not shown). As mentioned above, this might be attributed to different optical thicknesses of the overlying air. To conclude, a significant part of the temperature difference shown in Fig. 10 might be explained with the aforementioned hypothesis. However, more sensitivity simulations and further diagnostics would be necessary to specify the contributions from the dynamics and the individual physics parameterizations. We think that the differences visible in Fig. 10 are significant enough to provide an argument for the use of the deepatmosphere dynamics in a nonhydrostatic, geometricheightbased general circulation model that extends into the thermosphere.
The examination of zonal climatological means of the temperature and the zonal wind component provides the first part of a comprehensive evaluation of UAICON. It served us as one guideline for the general performance of UAICON and the identification of deficiencies during the main development phase. This development has now reached a stage which we regard as mature enough to form the first completed version of UAICON. Nevertheless, we are still facing significant biases in certain regions of the middle and upper atmosphere, which are especially apparent in the zonal wind component, as mentioned earlier. We assume that a retuning of some parameters of the parameterizations of orographically and nonorographically forced unresolved gravity waves is a promising, although not straightforward, route to reduce the biases.
An upperatmosphere extension of the ICON general circulation model has been presented. This includes the extension of the dynamical core from a shallowatmosphere to a deepatmosphere formulation in order to account for the spherical shape of the atmosphere and the gravitational field as well as to account for the nontraditional part of the Coriolis acceleration. In addition, the physics parameterizations have been complemented by processes which become relevant in the rarified air of the upper mesosphere and lower thermosphere. For instance, molecular diffusion takes over the lead role for mixing from turbulence, and various processes are linked to the relatively strong highfrequency solar irradiance in the upper atmosphere: its absorption is a source of heat; air compounds being ionized by the radiation align with the electromagnetic field of the Earth, forming currents that in turn interact with the neutral air flow via ion drag and Joule heating; atoms and radicals, as the products of photolysis, undergo a recombination reaction that is accompanied by chemical heating.
The new implementations, subsumed under the configuration name UAICON, have been validated by means of idealized test cases in terms of the dynamical core, and by climate simulations the results of which were compared to satellite observations. Two test cases were performed. The first one follows a method proposed by Läuter et al. (2005) and considers the propagation of a spherically symmetric sound wave in an atmosphere at rest in an absolute frame and so in a state of solidbody rotation relative to the rotating Earth. This requires accounting for the centrifugal acceleration explicitly in the dynamical core. In addition, gravity is switched off in order to simplify the derivation of an analytical solution, the comparison to which allows to quantify, e.g., how well the isotropic sound wave propagation is reproduced by the dynamical core on the anisotropic spherically curved grid, and how well it maintains the solidbody rotation of the air. The second test case is the Jablonowski–Williamson baroclinic instability test case in its formulation for the deep atmosphere by Ullrich et al. (2014). Its focus is on testing the representation of synopticscale flows through the simulation of the life cycle of a baroclinic wave on a hydrostatically and geostrophically balanced, zonally symmetric background state of the atmosphere. As with the first test case, the smallEarth approach of Wedi and Smolarkiewicz (2009) is applied to amplify the effects of the spherical curvature. No analytical solution is available for the second test case, so it relies on a model intercomparison. In either test case, UAICON showed satisfying performance, and no indication of a severe deficiency of the deepatmosphere modification of the dynamical core was found.
For the evaluation of the upperatmosphere physics, four AMIPtype simulations, each spanning 20 years, were performed and compared to temperature and zonal wind climatologies obtained from measurements of the SABER and TIMED satellite instruments and the URAP project, respectively. The first simulation uses a setup which we regard as typical for our envisaged first applications of UAICON. The model top is located at a height of 150 km, and the newly implemented upperatmosphere physics are enabled in addition to the standard physics of ICON. The deepatmosphere modifications of the dynamical core are applied as well. In the second simulation, the model top is lowered to 80 km, and the sponge layer is adjusted accordingly. Apart from that, the settings are identical to those of the first simulation. The third simulation equals the second simulation but is a standard ICON simulation with the upperatmosphere physics and the deepatmosphere modifications of the dynamics switched off. Finally, a fourth simulation is identical to the simulation with a model top at 150 km, except for the deepatmosphere modification of the dynamical core being switched off. By comparing the four simulation variants, we try to quantify the effects of the vertical extension of the model domain into the lower thermosphere and of the upperatmosphere physics on the middle and lower atmosphere. In addition, the difference between the deepatmosphere and the shallowatmosphere dynamics is quantified.
The temperature climatologies for January and July from the first type of simulation with UAICON are generally in good agreement with the observations. For instance, the atmospheric temperature minimum at the summer mesopause is relatively well captured. Some biases remain, however; a slightly toowarm winter stratosphere and toothick tropical tropopause layer might be due to the absence of stratospheric aerosol in the simulation or result from the employed ozone climatology. When it comes to the climatologies of the zonal wind for January and July, we find again that UAICON reproduces the qualitative structure of the wind in that part of the atmosphere observed by URAP. One remaining issue is that the westerly jets in the lower thermosphere are stronger in UAICON, roughly by 50 %. The comparison of UAICON simulations and two different model configurations with a lower model top at 80 km has shown that the addition of upperatmosphere physics and dynamics also affects stratospheric temperatures (increasing them by up to 5 K), and that the vertical extension has noticeable effects at least down to about 60 km. The difference between the shallowatmosphere and the deepatmosphere dynamics manifests itself in a significant decrease in temperature in the uppermesosphere–lowerthermosphere region if the latter is applied.
As an important application for UAICON, we have in view the investigation of the impact of gravity waves on the global atmospheric circulation with a special focus on the feedback of the middle and upper atmosphere on the tropospheric weather and climate. The high scalability of ICON on massively parallel computers will allow us to employ model configurations with much higher grid resolutions and eventually resolve large parts of the atmospheric gravity wave spectrum instead of parameterizing it.
The ICONSoftware is freely available to the scientific community for noncommercial research purposes under a license of DWD and MPIM. The license in its current form can be viewed on ICON developers (2018). UAICON existed in two configurations at the time of writing this article: one for the upperatmosphere developments within the development framework of the ICONSoftware for climate simulations and another one for the upperatmosphere developments within the NWP development framework of the ICONSoftware. If you would like to obtain UAICON, please contact icon@dwd.de. You will be provided an institutional license, which needs to be signed by the representative of your research institute and sent back to the DWD. The ICONSoftware is controlled by a Git version control system, and upon license agreement, a tar ball of the version of UAICON that was used to produce the results presented in this article is provided to you. The data shown in Sect. 3 as well as the scripts used for their production, postprocessing, and plotting are provided as a zipped folder in the Supplement of this article.
We assume the data used to define the initial state of the model atmosphere of UAICON to originate from hydrostatic models (this is true for the IFS data, as well as for the climatology employed here). The altitude assigned to such a data point is actually the geopotential height. This is no issue if the shallowatmosphere approximation is made and the geopotential height coincides with the geometric height employed in the dynamical core of ICON. However, in a deepatmosphere model, the geometric and geopotential heights differ. Therefore, the data initialization of UAICON takes place on geopotential heights, to and from which the geometric heights of the grid levels are transformed using
The geopotential height below which initial data, e.g., from the IFS model, are available, will be denoted 𝔷_{g} in the following (we use a value of 𝔷_{g}=70 km). Climatological temperatures are taken from Fleming et al. (1988), who offer tables with zonally averaged, monthly temperature values, denoted T_{F}, from mean sea level to a geopotential height of 120 in 5 km intervals. For our current simple approach, these datasets are averaged temporally and meridionally to obtain a single vertical temperature profile (see Fig. A1). Temperature values within the 5 km intervals are computed by a cubic spline interpolation (e.g., Bronstein et al., 2001). Above 120 km, an analytical temperature profile from Bates (1959) (see also Hedin, 1983) is used, which is formally identical to Eq. (7):
where ${T}_{\mathrm{120}\phantom{\rule{0.125em}{0ex}}\mathrm{km}}={T}_{F}({z}_{\mathrm{g}}=\mathrm{120}\phantom{\rule{0.125em}{0ex}}\mathrm{km})$, and T_{∞} is approximately the temperature for the limit z→∞. This limit corresponds to a geopotential height of z_{g}=a, which follows from Eq. (A1) by multiplying the righthand side of the first equation with $\mathrm{1}=(\mathrm{1}/z)/(\mathrm{1}/z)$ and applying the limit. The value of T_{∞} could be set, for instance, to the mean exospheric temperature of about 1035 K (in practice, we use a value of 400 K, since higher values could challenge the numerical stability in the initial phase of the simulation). For a steady transition between T_{F} and T_{B}, the scale height is set to ${H}_{B}=({T}_{\mathrm{\infty}}{T}_{\mathrm{120}\phantom{\rule{0.125em}{0ex}}\mathrm{km}})/(\mathrm{d}{T}_{F}/\mathrm{d}{z}_{\mathrm{g}}{)}_{{z}_{\mathrm{g}}=\mathrm{120}\phantom{\rule{0.125em}{0ex}}\mathrm{km}}$. In addition, the temperature blending requires extrapolating the temperature data from below 𝔷_{g} to above, which is done by a simple linear extrapolation:
with $\mathit{\gamma}=\mathrm{d}{T}_{\text{IFS}}/\mathrm{d}{z}_{\mathrm{g}}{}_{{z}_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}}}$. To obtain a statically stable stratification, γ is limited by the dry adiabatic lapse rate: $\mathit{\gamma}\ge {\mathrm{\Gamma}}_{\mathrm{d}}=g/{c}_{p}$ (e.g., Holton, 2004). The blending reads
where T_{clim} is T_{F} or T_{B}, respectively, and H_{blend} is a tunable blending scale height, which allows to control over what distance the transition between T_{IFS} and T_{clim} takes place. Currently, we use a value of H_{blend}=10 km to avoid negative absolute temperatures which could result from Eq. (A3), although somewhat larger values for H_{blend} do likely satisfy this as well. The blending factor α satisfies $\mathrm{d}\mathit{\alpha}/\mathrm{d}{z}_{\mathrm{g}}{}_{{z}_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}},{\mathfrak{z}}_{\mathrm{g}}+{H}_{\text{blend}}}=\mathrm{0}$ in order to guarantee a steady transition at z_{g}=𝔷_{g} and ${z}_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}}+{H}_{\text{blend}}$. Given the temperature field (A4), the other variables above 𝔷_{g} are determined by the hydrostatic and geostrophic balances. On the one hand, this provides a relatively simple method for their computation, and on the other hand it helps to reduce the magnitude of the dynamic tendencies and therefore the strength of the adjustment process during the first time steps of the numerical integration. The pressure is computed from a numerical integration of the discretized form of $\partial p/\partial {z}_{\mathrm{g}}+gp/\left(RT\right)=\mathrm{0}$, starting at 𝔷_{g} with $p(\mathit{\lambda},\mathit{\phi},{z}_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}})={p}_{\text{IFS}}(\mathit{\lambda},\mathit{\phi},{z}_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}})$, where the deepatmospherespecific terms, underlined in Eq. (10), are neglected to simplify matters. Once temperature and pressure are known, ρ, π, and θ can be diagnosed. The horizontal wind is determined from a blending formally identical to the temperature blending (Eq. A4) (using the same α). The IFS part for the blending above 𝔷_{g} is a simple linear extrapolation of the horizontal velocity, formally identical to Eq. (A3). The “climatological” part shall satisfy ${f}_{\mathrm{r}}{\mathit{e}}_{\mathrm{r}}\times {\mathit{v}}_{\mathrm{h},\text{clim}}=\mathit{\beta}\left(\mathit{\phi}\right)\left({\mathbf{\nabla}}_{\mathrm{h}}p\right)/\mathit{\rho}$, where p and ρ are the hydrostatic pressure and density, respectively. Associated with the thermal wind balance (Zdunkowski and Bott, 2003; Holton, 2004), relatively strong horizontal temperature gradients between 𝔷_{g} and 𝔷_{g}+H_{blend} can cause the magnitude of v_{h,clim} to increase with height and reach values which violate the CFL stability criterion (Zdunkowski and Bott, 2003; Holton, 2004). To avoid this, v_{h,clim} is multiplied by a factor $[\mathrm{1}+({z}_{\mathrm{g}}{\mathfrak{z}}_{\mathrm{g}})/{H}_{{v}_{\mathrm{h}}}]\mathrm{exp}[({z}_{\mathrm{g}}{\mathfrak{z}}_{\mathrm{g}})/{H}_{{v}_{\mathrm{h}}}]$, with a tuneable decay scale height ${H}_{{v}_{\mathrm{h}}}$ (we use a value of ${H}_{{v}_{\mathrm{h}}}=\mathrm{10}$ km). The value of this factor and its vertical derivative is 1 and 0 at z_{g}=𝔷_{g}, respectively, so as not to affect the continuity of the extrapolated wind at that altitude. Of course this factor causes v_{h,clim} to violate the geostrophic balance to some extent, but it turned out that this is less severe for the spinup phase than an initial wind field violating the CFL criterion locally. In addition, the horizontal pressure gradient is multiplied by the factor
in order to reduce its magnitude smoothly to zero towards the Equator, where the geostrophic balance does not apply. ${\mathit{\phi}}_{\text{trop}}>\mathrm{0}{}^{\circ}$ is a tuneable tropical latitude (we use ${\mathit{\phi}}_{\text{trop}}=\mathrm{10}{}^{\circ}$). Finally, the vertical wind above 𝔷_{g} is computed from a blending, again formally identical to Eq. (A4), of a linearly extrapolated w_{IFS} and a w_{clim}=0, in accordance with the hydrostatic balance and the boundary condition of a vanishing vertical wind at the model top.
B1 Formulation in the rotating frame
Here, we consider the transformation of a solution to the deepatmosphere equations, linearized about an atmosphere at rest in the absolute frame, to a frame rotating with angular velocity Ω. For the atmosphere at rest in the absolute frame, ${\mathit{v}}_{\mathrm{a}}={\mathit{v}}_{F}+\mathit{v}=\mathrm{0}$ holds, where ${\mathit{v}}_{F}=\mathbf{\Omega}\times \mathit{r}$, v being the velocity observed in the rotating frame, and $\mathit{r}=\mathit{X}\mathit{A}$ (with the center of the Earth A and an arbitrary point X not coincident with A). Therefore, in the rotating frame, the motion of the air follows
where Eqs. (B1) and (B2) are equivalent. In the last step, we have introduced the identity tensor 1 and the antisymmetric Coriolis tensor $\mathbf{W}=\mathbf{\Omega}\times \mathbf{1}$ to simplify the following considerations (compare Wilson, 1929; Zdunkowski and Bott, 2003). The solution to Eq. (B2) is
where ${\mathit{r}}_{\mathrm{0}}=\mathit{r}(t=\mathrm{0})$. Introducing $\widehat{\mathbf{\Omega}}=\mathbf{\Omega}/\mathrm{\Omega}$, and $\widehat{\mathbf{W}}=\mathbf{W}/\mathrm{\Omega}$, with $\mathrm{\Omega}=\left\mathbf{\Omega}\right$ and using
where ${\widehat{\mathbf{W}}}^{\mathrm{2}}=\widehat{\mathbf{W}}\cdot \widehat{\mathbf{W}}$, for instance, in combination with
where $\widehat{\mathbf{1}}=\mathbf{1}\widehat{\mathbf{\Omega}}\widehat{\mathbf{\Omega}}$ and $\widehat{\mathbf{\Omega}}\widehat{\mathbf{\Omega}}$ denote the dyadic product of $\widehat{\mathbf{\Omega}}$ with itself, the exponential factor in Eq. (B3) can be expanded as
So Eq. (B3) can be rewritten in a form which is easier to evaluate (compare Wilson, 1929; Bronstein et al., 2001; Zdunkowski and Bott, 2003):
Next, given, for instance, an analytical solution in the absolute frame, which is spherically symmetric with respect to the point B, therefore depending only on $\mathit{X}\mathit{B}=\left\right(\mathit{X}\mathit{A})(\mathit{B}\mathit{A}\left)\right=\mathit{\left\{}\right[(\mathit{X}\mathit{A})(\mathit{B}\mathit{A})]\cdot [(\mathit{X}\mathit{A})(\mathit{B}\mathit{A})]{\mathit{\}}}^{\mathrm{1}/\mathrm{2}}$, we would find the solution unaltered in the rotating frame relative to $\mathrm{exp}(\mathbf{W}t)\cdot (\mathit{B}\mathit{A})$ since, e.g.,
Here, (⋅)^{⊺} denotes the transpose, and we have used the identity $\mathrm{exp}(\mathbf{W}t{)}^{\mathit{\u22ba}}=\mathrm{exp}(\mathbf{W}t)$ which can be derived from Eq. (B6) and $\widehat{\mathbf{\Omega}}{\widehat{\mathbf{\Omega}}}^{\mathit{\u22ba}}=\widehat{\mathbf{\Omega}}\widehat{\mathbf{\Omega}}$, 1^{⊺}=1 and ${\widehat{\mathbf{W}}}^{\mathit{\u22ba}}=\widehat{\mathbf{W}}$. In addition, $\mathrm{exp}\left(\mathbf{W}t\right)\cdot \mathrm{exp}(\mathbf{W}t)=\mathbf{1}$ can be shown to hold, using Eqs. (B6), (B5), and the definition of $\widehat{\mathbf{W}}$. The above is not a mathematically rigorous proof for the general case, but we hope it can help to illustrate the basic idea behind the method. For detailed mathematical examinations, we refer the reader to Läuter et al. (2005); Staniforth and White (2008). In the following section, such a spherically symmetric solution to be rotated is derived.
B2 Derivation of the sound wave solution
Here, we derive the analytical solution for the expansion of sound waves for deepatmosphere Eqs. (1)–(4) reduced to the reversible processes (i.e., F=0 and Q=0), without gravity (g=0), linearized about an isothermal atmosphere of temperature T_{0}=const., at rest in the absolute frame (i.e., ${\mathit{v}}_{\mathbf{a}}={\mathit{v}}_{F}+{\mathit{v}}_{\mathrm{0}}=\mathbf{\Omega}\times \mathit{r}+{\mathit{v}}_{\mathrm{0}}=\mathrm{0}$). In particular, the neglect of gravity leads to a constant base state for pressure π_{0}=const. (from Eq. 1) and density ρ_{0}=const. (from Eq. 4), too. Otherwise, the derivation of an analytic solution even for the linearized equations would be very difficult.
Linearization of Eqs. (1) to (4) about this base state leads to
resulting in the sound wave equation:
with the squared speed of sound ${c}_{\mathrm{s}}^{\mathrm{2}}=({c}_{p}/{c}_{v})\mathrm{R}{T}_{\mathrm{0}}$.
Without gravity and disregarding the spherical boundaries of the atmosphere, there is no longer a distinct direction (or rather a distinct point). Therefore, we seek a solution which is spherically symmetric with respect to an arbitrary point B within the boundaries. Using the expansion of the Laplace operator in spherical coordinates (e.g., Bronstein et al., 2001), we can write
with the radial coordinate $x=\mathit{X}\mathit{B}$ (to avoid any confusion with the distance r from the center of the Earth). This equation is valid for x>0. Using the transform
one finds for x⋅ (Eq. B13) (Nolting, 2004)
In addition, we can derive from Eq. (B11)
where ${\mathit{v}}^{\prime}={v}^{\prime}{\mathit{e}}_{x}$, ${\mathit{e}}_{x}=(\mathit{X}\mathit{B})/\mathit{X}\mathit{B}$, and $\stackrel{\mathrm{\u0303}}{v}=x{v}^{\prime}$. Equation (B16) is required for the specification of the initial conditions. Given 1D wave Eq. (B15), a general solution is of the form (Nolting, 2004)
where f_{1} and f_{2} denote incoming and outgoing waves of arbitrary shape, respectively. Note that the transform (Eq. B14) induces a boundary condition at x=0, because finite values of π^{′} require $\stackrel{\mathrm{\u0303}}{\mathit{\pi}}(x=\mathrm{0},t)=\mathrm{0}$. Therefore, the outgoing wave for later times is determined by the incoming wave at the origin:
Since wave Eq. (B15) is of second order in time, initial conditions for $\stackrel{\mathrm{\u0303}}{\mathit{\pi}}(x,t=\mathrm{0})$ and for its derivative
are required. Instead, we require that the incoming wave vanishes (f_{1}=0) and prescribe the initial wind field:
where δv is a constant wind amplitude, n denotes the number of wave crests, and $x={b}_{\mathrm{1}}>\mathrm{0}$ and $x={b}_{\mathrm{2}}>{b}_{\mathrm{1}}$ are the radii, between which the sound wave has a nonvanishing amplitude at t=0. The Heaviside step function is given by (Bronstein et al., 2001)
The function in Eq. (B19) was chosen not least because it satisfies ${v}^{\prime}(x={b}_{\mathrm{1}},{b}_{\mathrm{2}})=\mathrm{0}$ and $\mathrm{d}{v}^{\prime}/\mathrm{d}x{}_{x={b}_{\mathrm{1}},{b}_{\mathrm{2}}}=\mathrm{0}$, and guarantees that the derived fields π^{′}, θ^{′}, and ρ^{′} do likewise, so that discontinuities pose no extra challenge to the numerical simulation. For the implementation of the test case, we preferred to specify the wave amplitude by means of a temperature amplitude δT, instead of the wind amplitude δv, and use the relation $\mathit{\delta}v=({c}_{v}/\mathrm{R})(\mathit{\delta}T/{T}_{\mathrm{0}}){c}_{\mathrm{s}}$ for this purpose. From the initial wind field (B19), we can calculate the derivative df_{2}(x)∕dx by Eqs. (B16) and (B18). Integration of df_{2}∕dx yields the Exner pressure perturbation (compare Bronstein et al., 2001):
or, alternatively, expressed as pressure perturbation:
The initialization of the remaining thermodynamic fields (here, e.g., ρ and θ) can be chosen arbitrarily as long as the linearized ideal gas law is fulfilled. From Eq. (4), it follows
In our implementation, we used ${\mathit{\rho}}^{\prime}(x,t=\mathrm{0})={c}_{v}{\mathit{\rho}}_{\mathrm{0}}$ (Eq. B21) ∕(Rπ_{0}) and determined ${\mathit{\theta}}^{\prime}(x,t=\mathrm{0})$ from Eq. (B23).
Of course, the solution outlined above holds only until the wave front, initially at x=b_{2}, impinges on the solid boundaries, either the model bottom or top, and is reflected. Therefore, the integration time for this test case is limited.
The above solution is what an observer would see in the absolute frame. To transform it into the rotating frame, every point X (including B) has to be rotated according to the formula
which follows from Eq. (B3). Alternatively, the formulation in Eq. (B7) could be used. Equation (B24) is the coordinateindependent formulation of the rotation. It has to be measured in the coordinate system of choice for the application in practice. As shown in Eq. (B8), the shape of the sound wave relative to B(t) is unaffected by the rotation.
Finally, we will give a short outline for the implementation of this test case into a model. First of all, the gravitational acceleration g has to be set to zero, or at least to a relatively small value. If rotation is switched on (Ω≠0), it is essential to take into account the centrifugal acceleration $\mathbf{\Omega}\times [\mathbf{\Omega}\times (\mathit{X}\mathit{A}\left)\right]$ explicitly in the dynamical core, as shown in Eq. (B1). If desired, the radius and angular velocity of the Earth are rescaled, a→η_{1}a and Ω→η_{2}Ω, and the topography is set to zero. If possible, the slip boundary condition should be applied at model bottom and top. All physics parameterizations should be switched off. After specifying the background state of the atmosphere (at rest in the absolute frame) by T_{0} and p_{0}, and the sound wave parameters B, b_{1}, b_{2}, δv or δT, and n, the dynamic fields of the (dry) atmosphere are initialized according to
In order to keep the nonlinear dynamics of the sound wave as simulated by the dynamical core negligibly small, the nondimensional wind amplitude $\mathit{\delta}v/{c}_{\mathrm{s}}=({c}_{v}/\mathrm{R})(\mathit{\delta}T/{T}_{\mathrm{0}})$ should be small enough. Apart from that, we adhered closely to the guidelines given by Baldauf et al. (2014) in our implementation of the test case into UAICON.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1235412019supplement.
The upperatmosphere extension of the ICON model was a close collaboration between the authors. The focus of GuZ and HS was on the adaption, implementation, and evaluation of the upperatmosphere physics parameterizations. The focus of SB, MB, GüZ, and DR was on the implementation and evaluation of the deepatmosphere modifications to the dynamical core.
The authors declare that they have no conflict of interest.
The authors (Sebastian Borchert, Günther Zängl, Michael Baldauf, Guidi Zhou, and Hauke Schmidt) thank the German Research Foundation (DFG) for partial support through the research unit Multiscale Dynamics of Gravity Waves (MSGWaves) and through grants ZA 268/101 and SCHM 2158/51. We are grateful to the entire ICON development team (MPIM, DWD, German Climate Computing Center (DKRZ), and Karlsruhe Institute of Technology (KIT)). We greatly appreciated many helpful discussions with Elisa Manzini (MPIM) on aspects of modeling the middle atmosphere. Special thanks go to Florian Prill (DWD), Marco A. Giorgetta (MPIM), and Sebastian Rast (MPIM) for many helpful and illuminating discussions on the dynamical core, the physics packages, and the code infrastructure of the ICON model. We thank the anonymous reviewers for very helpful comments, which led to significant improvements of the manuscript. For his help during the review process, we thank the editor, Juan Antonio Añel. Finally, we are grateful to the URAP and SABER instrument teams for creating and maintaining the satellite datasets.
This research has been supported by the Deutsche Forschungsgemeinschaft (grant nos. ZA 268/101 and SCHM 2158/51).
This paper was edited by Juan Antonio Añel and reviewed by three anonymous referees.
Baldauf, M. and Brdar, S.: An analytic solution for linear gravity waves in a channel as a test for numerical models using the nonhydrostatic, compressible Euler equations, Q. J. Roy. Meteor. Soc., 139, 1977–1989, https://doi.org/10.1002/qj.2105, 2014. a
Baldauf, M. and Brdar, S.: 3D diffusion in terrainfollowing coordinates: testing and stability of horizontally explicit, vertically implicit discretizations, Q. J. Roy. Meteor. Soc., 142, 2087–2101, https://doi.org/10.1002/qj.2805, 2016. a
Baldauf, M., Reinert, D., and Zängl, G.: An analytical solution for linear gravity and sound waves on the sphere as a test for compressible, nonhydrostatic numerical models, Q. J. Roy. Meteor. Soc., 140, 1974–1985, https://doi.org/10.1002/qj.2105, 2014. a, b, c, d, e, f, g, h
Banks, P. M. and Kockarts, G.: Aeronomy, Part B, Elsevier, 1973. a, b
Bates, D. R.: Some problems concerning the terrestrial atmosphere above the 100 km level, Philos. T. R. Soc. Lond., 253, 451–462, 1959. a, b, c
Becker, E.: Sensitivity of the upper mesosphere to the Lorenz energy cycle of the troposphere, J. Atmos. Sci., 66, 647–666, https://doi.org/10.1175/2008JAS2735.1, 2009. a
Becker, E. and Vadas, S. L.: Secondary gravity waves in the winter mesosphere: Results from a highresolution global circulation model, J. Geophys. Res.Atmos., 123, 2605–2627, https://doi.org/10.1002/2017JD027460, 2018. a
Bonaventura, L. and Ringler, T.: Analysis of discrete shallowwater models on geodesic Delaunay grids with Ctype staggering, Mon. Weather Rev., 133, 2351–2373, https://doi.org/10.1175/MWR2986.1, 2005. a
Bronstein, I. N., Semendjajew, K. A., Musiol, G., and Mühlig, H.: Taschenbuch der Mathematik, Verlag Harri Deutsch, Frankfurt am Main, 5 edn., 2001. a, b, c, d, e, f
CharltonPerez, A. J., Baldwin, M. P., Birner, T., Black, R. X., Butler, A. H., and Calvo, N.: On the lack of stratospheric dynamical variability in lowtop versions of the CMIP5 models, J. Geophys. Res.Atmos., 118, 2494–2505, https://doi.org/10.1002/jgrd.50125, 2013. a
Crueger, T., Giorgetta, M. A., Brokopf, R., Esch, M., Fiedler, S., Hohenegger, C., Kornblueh, L., Mauritsen, T., Nam, C., Naumann, A. K., Peters, K., Rast, S., Roeckner, E., Sakradzija, M., Schmidt, H., Vial, J., Vogel, R., and Stevens, B.: ICONA, the atmosphere component of the ICON Earth system model: II. Model evaluation, J. Adv. Model. Earth Sy., 10, 1638–1662, https://doi.org/10.1029/2017MS001233, 2018. a, b, c, d, e
Davies, T., Cullen, M., Malcolm, A., Mawson, M., Staniforth, A., White, A., and Wood, N.: A new dynamical core for the Met Office's global and regional modelling of the atmosphere, Q. J. Roy. Meteor. Soc., 131, 1759–1782, https://doi.org/10.1256/qj.04.101, 2005. a
Dawkins, E. C. M., Feofilov, A., Rezac, L., Kutepov, A. A., Janches, D., Höffner, J., Chu, X., Lu, X., Mlynczak, M. G., and Russell, J.: Validation of SABER v2.0 operational temperature data with groundbased lidars in the mesospherelower thermosphere region (75–105 km), J. Geophys. Res.Atmos., 123, 9916–9934, https://doi.org/10.1029/2018JD028742, 2018. a
Dipankar, A., Stevens, B., Heinze, R., Moseley, C., Zängl, G., Giorgetta, M. A., and Brdar, S.: Large eddy simulation using the general circulation model ICON, J. Adv. Model. Earth Sy., 7, 963–986, https://doi.org/10.1002/2015MS000431, 2015. a, b
Fleming, E. L., Chandra, S., Shoeberl, M. R., and Barnett, J. J.: Monthly mean global climatology of temperature, wind, geopotential height and pressure for 0–120 km, Tech. Rep. NASATM100697, NASA, 1988. a, b, c
Fomichev, V. I. and Blanchet, J.P.: Development of the new CCC/GCM longwave radiation model for extension into the middle atmosphere, Atmos. Ocean, 33, 513–529, https://doi.org/10.1080/07055900.1995.9649543, 1995. a, b
Fomichev, V. I. and Blanchet, J.P.: Matrix parameterization of the 15 pm CO_{2} band cooling in the middle and upper atmosphere for variable CO_{2} concentration, J. Geophys. Res., 103, 11505–11528, 1998. a, b
Fomichev, V. I., Ward, W. E., Beagley, S. R., McLandress, C., McConnell, J. C., McFarlane, N. A., and Shepherd, T. G.: Extended Canadian Middle Atmosphere Model: Zonalmean climatology and physical parameterizations, J. Geophys. Res., 107, 4087, https://doi.org/10.1029/2001JD000479, 2002. a
Fritts, D. C. and Alexander, M. J.: Gravity wave dynamics and effects in the middle atmosphere, Rev. Geophys., 41, 1003, https://doi.org/10.1029/2001RG000106, 2003. a
Gassmann, A.: Inspection of hexagonal and triangular Cgrid discretizations of the shallow water equations, J. Comput. Phys., 230, 2706–2721, https://doi.org/10.1016/j.jcp.2011.01.014, 2011. a
Gassmann, A.: Discretization of generalized Coriolis and friction terms on the deformed hexagonal Cgrid, Q. J. Roy. Meteor. Soc., 144, 2038–2053, https://doi.org/10.1002/qj.3294, 2018. a
Gill, A. E.: Atmosphereocean dynamics, Academic Press, London, 4th Edn., 1982. a, b, c, d, e
Giorgetta, M. A., Brokopf, R., Crueger, T., Esch, M., Fiedler, S., Helmert, J., Hohenegger, C., Kornblueh, L., Köhler, M., Manzini, E., Mauritsen, T., Nam, C., Raddatz, T., Rast, S., Reinert, D., Sakradzija, M., Schmidt, H., Schneck, R., Schnur, R., Silvers, L., Wan, H., Zängl, G., and Stevens, B.: ICONA, the atmosphere component of the ICON Earth system model: I. Model description, J. Adv. Model. Earth Sy., 10, 1613–1637, https://doi.org/10.1029/2017MS001242, 2018. a, b, c, d, e
Hedin, A. E.: A revised thermospheric model based on mass spectrometer and incoherent scatter data: MSIS83, J. Geophys. Res., 88, 10170–10188, 1983. a, b, c
Heinze, R., Dipankar, A., Henken, C. C., Moseley, C., Sourdeval, O., Trömel, S., Xie, X., Adamidis, P., Ament, F., Baars, H., Barthlott, C., Behrendt, A., Blahak, U., Bley, S., Brdar, S., Brueck, M., Crewell, S., Deneke, H., Di Girolamo, P., Evaristo, R., Fischer, J., Frank, C., Friederichs, P., Göcke, T., Gorges, K., Hande, L., Hanke, M., Hansen, A., Hege, H.C., Hoose, C., Jahns, T., Kalthoff, N., Klocke, D., Kneifel, S., Knippertz, P., Kuhn, A., van Laar, T., Macke, A., Maurer, V., Mayer, B., Meyer, C. I., Muppa, S. K., Neggers, R. A. J., Orlandi, E., Pantillon, F., Pospichal, B., Röber, N., Scheck, L., Seifert, A., Seifert, P., Senf, F., Siligam, P., Simmer, C., Steinke, S., Stevens, B., Wapler, K., Weniger, M., Wulfmeyer, V., Zängl, G., Zhang, D., and Quaas, J.: Largeeddy simulations over Germany using ICON: a comprehensive evaluation, Q. J. Roy. Meteor. Soc., 143, 69–100, https://doi.org/10.1002/qj.2947, 2017. a
Hines, C. O.: Dopplerspread parameterization of gravitywave momentum deposition in the middle atmosphere. Part 1: Basic formulation, J. Atmos. Sol.Terr. Phy., 59, 371–386, https://doi.org/10.1016/S13646826(96)00079X, 1997a. a, b
Hines, C. O.: Dopplerspread parameterization of gravitywave momentum deposition in the middle atmosphere. Part 2: Broad and quasi monochromatic spectra, and implementation, J. Atmos. Sol.Terr. Phy., 59, 387–400, https://doi.org/10.1016/S13646826(96)00079X, 1997b. a, b, c
Holton, J. R.: An introduction to dynamic meteorology, Elsevier Academic Press, Oxford, UK, 4th Edn., 2004. a, b, c, d
Hong, S.S. and Lindzen, R. S.: Solar Semidiurnal Tide in the Thermosphere, J. Atmos. Sci., 33, 135–153, https://doi.org/10.1175/15200469(1976)033<0135:SSTITT>2.0.CO;2, 1976. a, b
Huang, T., Walters, S., Brasseur, G., Hauglustaine, D., and Wu, W.: Description of SOCRATES: A chemical dynamical radiative twodimensional model, Tech. Rep. NCAR/TN440+EDD, NCAR, Boulder, CO, 1998. a, b
Iacono, M. J., Delamere, J. S., Mlawer, E. J., Shephard, M. W., Clough, S. A., and Collins, W. D.: Radiative forcing by longlived greenhouse gases: Calculations with the AER radiative transfer models, J. Geophys. Res.Atmos., 113, d13103, https://doi.org/10.1029/2008JD009944, 2008. a
ICON developers: ICON : Icosahedral Nonhydrostatic Weather and Climate Model, available at: https://code.mpimet.mpg.de/projects/iconpublic, last access: 14 November 2018. a
Jablonowski, C. and Williamson, D. L.: A baroclinic instability test case for atmospheric model dynamical cores, Q. J. Roy. Meteor. Soc., 132, 2943–2975, https://doi.org/10.1256/qj.06.12, 2006. a, b
Kasahara, A.: The roles of the horizontal component of the Earth's angular velocity in nonhydrostatic linear models, J. Atmos. Sci., 60, 1085–1095, https://doi.org/10.1175/15200469(2003)60<1085:TROTHC>2.0.CO;2, 2003. a, b
Kim, Y.J., Eckermann, S. D., and Chun, H.Y.: An overview of the past, present and future of gravitywave drag parametrization for numerical climate and weather prediciton models, Atmos. Ocean, 40, 65–98, https://doi.org/10.3137/ao.410105, 2003. a
Kinne, S., O'Donnel, D., Stier, P., Kloster, S., Zhang, K., Schmidt, H., Rast, S., Giorgetta, M., Eck, T. F., and Stevens, B.: MACv1: A new global aerosol climatology for climate studies, J. Adv. Model. Earth Sy., 5, 704–740, https://doi.org/10.1002/jame.20035, 2013. a
Kinnison, D. E., Brasseur, G. P., Walters, S., Garcia, R. R., Marsh, D. R., Sassi, F., Harvey, V. L., Randall, C. E., Emmons, L., Lamarque, J. F., Hess, P., Orlando, J. J., Tie, X. X., Randel, W., Pan, L. L., Gettelman, A., Granier, C., Diehl, T., Niemeier, U., and Simmons, A. J.: Sensitivity of chemical tracers to meteorological parameters in the MOZART3 chemical transport model, J. Geophys. Res.Atmos., 112, D20302, https://doi.org/10.1029/2006JD007879, 2007. a
Kirchhoff, G.: Vorlesung über mathematische Physik: Mechanik, B. G. Teubner, Leipzig, 1876. a
Klemp, J. B., Dudhia, J., and Hassiotis, A. D.: An upper gravitywave absorbing layer for NWP applications, Mon. Weather Rev., 136, 3987–4004, https://doi.org/10.1175/2008MWR2596.1, 2008. a
Kockarts, G.: Nitric oxide cooling in the terrestrial thermosphere, Geophys. Res. Lett., 7, 137–140, https://doi.org/10.1029/GL007i002p00137, 1980. a, b
Läuter, M., Handorf, D., and Dethloff, K.: Unsteady analytical solutions of the spherical shallow water equations, J. Comput. Phys., 210, 535–553, https://doi.org/10.1016/j.jcp.2005.04.022, 2005. a, b, c, d
Liu, H.L., Bardeen, C. G., Foster, B. T., Lauritzen, P., Liu, J., Lu, G., Marsh, D. R., Maute, A., McInerney, J. M., Pedatella, N. M., Qian, L., Richmond, A. D., Roble, R. G., Solomon, S. C., Vitt, F. M., and Wang, W.: Development and validation of the Whole Atmosphere Community Climate Model with thermosphere and ionosphere extension (WACCMX 2.0), J. Adv. Model. Earth Sy., 10, 381–402, https://doi.org/10.1002/2017MS001232, 2018. a
Manzini, E. and McFarlane, N. A.: The effect of varying the source spectrum of a gravity wave parameterization in a middle atmosphere general circulation model, J. Geophys. Res., 103, 31523–31539, https://doi.org/10.1029/98JD02274, 1998. a
Manzini, E., McFarlane, N. A., and McLandress, C.: Impact of the Doppler spread parameterization on the simulation of the middle atmosphere circulation using the MA/ECHAM4 general circulation model, J. Geophys. Res., 102, 25751–25762, https://doi.org/10.1029/97JD01096, 1997. a
Miyahara, S., Yoshida, Y., and Miyoshi, Y.: Dynamic coupling between the lower and upper atmosphere by tides and gravity waves, J. Atmos. Terr. Phys., 55, 1039–1053, https://doi.org/10.1016/00219169(93)90096H, 1993. a
Mlawer, E. J., Taubman, S. J., Brown, P. D., Iacono, M. J., and Clough, S. A.: Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlatedk model for the longwave, J. Geophys. Res.Atmos., 102, 16663–16682, https://doi.org/10.1029/97JD00237, 1997. a
Mlynczak, M. G. and Solomon, S.: A detailed evaluation of the heating efficiency in the middle atmosphere, J. Geophys. Res., 98, 10517–10541, https://doi.org/10.1029/93JD00315, 1993. a, b
Nolting, W.: Grundkurs theoretische Physik 3: Elektrodynamik, SpringerVerlag, Berlin, 7th Edn., 2004. a, b
Ogibalov, V. and Fomichev, V.: Parameterization of solar heating by the near IR CO_{2} bands in the mesosphere, Adv. Space Res., 32, 759–764, https://doi.org/10.1016/S02731177(03)800698, 2003. a, b
Phillips, N. A.: The equations of motion for a shallow rotating atmosphere and the “traditional approximation”, J. Atmos. Sci., 23, 626–628, 1966. a, b, c, d
Pincus, R. and Stevens, B.: Paths to accuracy for radiation parameterizations in atmospheric models, J. Adv. Model. Earth Sy., 5, 225–233, https://doi.org/10.1002/jame.20027, 2013. a
Richards, P. G., Torr, M. R., and Torr, D. G.: The seasonal effect of nitric oxide cooling on the thermospheric U.V. heat budget, Planet. Space Sci., 30, 515–518, https://doi.org/10.1016/00320633(82)900629, 1982. a
Richards, P. G., Fennelly, J. A., and Torr, D. G.: EUVAC: A solar EUV Flux Model for aeronomic calculations, J. Geophys. Res.Space., 99, 8981–8992, https://doi.org/10.1029/94JA00518, 1994. a, b
Richter, J. H., Sassi, F., Garcia, R. R., Matthes, K., and Fischer, C. A.: Dynamics of the middle atmosphere as simulated by the Whole Atmosphere Community Climate Model, version 3 (WACCM3), J. Geophys. Res.Atmos., 113, D08101, https://doi.org/10.1029/2007JD009269, 2008. a
Roble, R. G.: Energetics of the Mesosphere and Thermosphere, American Geophysical Union, 1–21, 1995. a
Roeckner, E., Brokopf, R., Esch, M., Giorgetta, M. A., Hagemann, S., Kornblueh, L., Manzini, E., Schlese, U., and Schulzweida, U.: Sensitivity of simulated climate to horizontal and vertical resolution in the ECHAM5 atmosphere model, J. Climate, 19, 3771–3791, https://doi.org/10.1175/JCLI3824.1, 2006. a
Sadourny, R.: Conservative finitedifference approximations of the primitve equations on quasiuniform spherical grids, Mon. Weather Rev., 100, 136–144, 1972. a
Scaife, A. A., Spangehl, T., Fereday, D. R., and Cubasch, U.: Climate change projections and stratospheretroposphere interaction, Clim. Dynam., 38, 2089–2097, https://doi.org/10.1007/s0038201110807, 2012. a
Schmidt, H., Brasseur, G. P., Charron, M., Manzini, E., Giorgetta, M. A., Diehl, T., Fomichev, V. I., Kinnison, D., Marsh, D., and Walters, S.: The HAMMONIA chemistry climate model: sensitivity of the mesopause region to the 11year solar cycle and CO_{2} doubling, J. Climate, 19, 3903–3931, https://doi.org/10.1175/JCLI3829.1, 2006. a, b, c
Shepherd, T. G., Semeniuk, K., and Koshyk, J. N.: Sponge layer feedbacks in middleatmosphere models, J. Geophys. Res., 101, 23447–23464, https://doi.org/10.1029/96JD01994, 1996. a
Skamarock, W. C. and Klemp, J. B.: A timesplit nonhydrostatic atmospheric model for weather research and forecasting applications, J. Comput. Phys., 227, 3465–3485, https://doi.org/10.1016/j.jcp.2007.01.037, 2008. a
Smolarkiewicz, P. K., Deconinck, W., Hamrud, M., Kühnlein, C., Mozdzynski, G., Szmelter, J., and Wedi, N. P.: A finitevolume module for simulating global allscale atmospheric flows, J. Comput. Phys., 314, 287–304, https://doi.org/10.1016/j.jcp.2016.03.015, 2016. a
Staniforth, A. and White, A. A.: Unsteady exact solutions of the flow equations for threedimensional spherical atmospheres, Q. J. Roy. Meteor. Soc., 134, 1615–1626, https://doi.org/10.1002/qj.300, 2008. a, b
Staniforth, A. and Wood, N.: The deepatmosphere Euler equations in a generalized vertical coordinate, Mon. Weather Rev., 131, 1931–1938, https://doi.org/10.1175//2564.1, 2003. a, b
Staniforth, A. and Wood, N.: Aspects of the dynamical core of a nonhydrostatic, deepatmosphere, unified weather and climateprediction model, J. Comput. Phys., 227, 3445–3464, https://doi.org/10.1016/j.jcp.2006.11.009, 2008. a, b, c
Stevens, B., Giorgetta, M., Esch, M., Mauritsen, T., Crueger, T., Rast, S., Salzmann, M., Schmidt, H., Bader, J., Block, K., Brokopf, R., Fast, I., Kinne, S., Kornblueh, L., Lohmann, U., Pincus, R., Reichler, T., and Roeckner, E.: Atmospheric component of the MPIM Earth System Model: ECHAM6, J. Adv. Model. Earth Sy., 5, 146–172, https://doi.org/10.1002/jame.20015, 2013. a
Strobel, D. F.: Parameterization of the atmospheric heating rate from 15 to 120 km due to O_{2} and O_{3} absorption of solar radiation, J. Geophys. Res., 83, 6225–6230, https://doi.org/10.1029/JC083iC12p06225, 1978. a, b
Swinbank, R. and Ortland, D. A.: Compilation of wind data for the Upper Atmosphere Research Satellite (UARS) reference atmosphere project, J. Geophys. Res.Atmos., 108, 4615, https://doi.org/10.1029/2002JD003135, 2003. a
Taylor, K. E., Williamson, D., and Zwiers, F.: The sea surface temperature and seaice concentration boundary conditions for AMIP II simulations, Tech. Rep. PCMDI Report No. 60, Program for Climate Model Diagnosis and Intercomparison, Lawrence Livermore National Laboratory, Livermore, Ca, 1998. a
Thompson, D. W. J., Baldwin, M. P., and Wallace, J. M.: Stratospheric connection to northern hemisphere wintertime weather: implications for prediction, J. Climate, 15, 1421–1428, https://doi.org/10.1175/15200442(2002)015<1421:SCTNHW>2.0.CO;2, 2002. a
Thuburn, J. and White, A. A.: A geometrical view of the shallowatmosphere approximation, with application to the semiLagrangian departure point calculation, Q. J. Roy. Meteor. Soc., 139, 261–268, https://doi.org/10.1002/qj.1962, 2013. a, b
Thuburn, J., Wood, N., and Staniforth, A.: Normal modes of deep atmospheres. I: Spherical geometry, Q. J. Roy. Meteor. Soc., 128, 1771–1792, https://doi.org/10.1256/003590002320603403, 2002a. a
Thuburn, J., Wood, N., and Staniforth, A.: Normal modes of deep atmospheres. II: f–Fplane geometry, Q. J. Roy. Meteor. Soc., 128, 1793–1806, https://doi.org/10.1256/003590002320603412, 2002b. a
Tomita, H. and Satoh, M.: A new dynamical framework of nonhydrostatic global model using the icosahedral grid, Fluid Dyn. Res., 34, 357–400, https://doi.org/10.1016/j.fluiddyn.2004.03.003, 2004. a
Tort, M. and Dubos, T.: Dynamically consistent shallowatmosphere equations with a complete Coriolis force, Q. J. Roy. Meteor. Soc., 140, 2388–2392, https://doi.org/10.1002/qj.2274, 2014. a
Ullrich, P. A. and Jablonowski, C.: MCore: A nonhydrostatic atmospheric dynamical core utilizing highorder finitevolume methods, J. Comput. Phys., 231, 5078–5108, https://doi.org/10.1016/j.jcp.2012.04.024, 2012. a
Ullrich, P. A., Melvin, T., Jablonowski, C., and Staniforth, A.: A proposed baroclinic wave test case for deep and shallowatmosphere dynamical cores, Q. J. Roy. Meteor. Soc., 140, 1590–1602, https://doi.org/10.1002/qj.2241, 2014. a, b, c, d, e, f, g, h, i
Ullrich, P. A., Jablonowski, C., Kent, J., Lauritzen, P. H., Nair, R., Reed, K. A., Zarzycki, C. M., Hall, D. M., Dazlich, D., Heikes, R., Konor, C., Randall, D., Dubos, T., Meurdesoif, Y., Chen, X., Harris, L., Kühnlein, C., Lee, V., Qaddouri, A., Girard, C., Giorgetta, M., Reinert, D., Klemp, J., Park, S.H., Skamarock, W., Miura, H., Ohno, T., Yoshida, R., Walko, R., Reinecke, A., and Viner, K.: DCMIP2016: a review of nonhydrostatic dynamical core design and intercomparison of participating models, Geosci. Model Dev., 10, 4477–4509, https://doi.org/10.5194/gmd1044772017, 2017. a, b
Vallis, G. K.: Atmospheric and oceanic fluid dynamics: fundamentals and largescale circulation, Cambridge University Press, New York, 2006. a, b, c
Walko, R. L. and Avissar, R.: The OceanLandAtmosphere Model (OLAM). Part II: Formulation and tests of the nonhydrostatic dynamic core, Mon. Weather Rev., 136, 4045–4062, https://doi.org/10.1175/2008MWR2523.1, 2008. a
Wan, H.: Developing and testing a hydrostatic atmospheric dynamical core on triangular grids, PhD thesis, Max Planck Institute for Meteorology, Hamburg, Germany, 2009. a, b
Wan, H., Giorgetta, M. A., Zängl, G., Restelli, M., Majewski, D., Bonaventura, L., Fröhlich, K., Reinert, D., Rípodas, P., Kornblueh, L., and Förstner, J.: The ICON1.2 hydrostatic atmospheric dynamical core on triangular grids – Part 1: Formulation and performance of the baseline version, Geosci. Model Dev., 6, 735–763, https://doi.org/10.5194/gmd67352013, 2013. a, b
Watanabe, S. and Miyahara, S.: Quantifiaction of the gravity wave forcing of the migrating diurnal tide in a gravity waveresolving general circulation model, J. Geophys. Res., 114, D07110, https://doi.org/10.1029/2008JD011218, 2009. a, b
Watanabe, S., Kawatani, Y., Tomikawa, Y., Miyazaki, K., Takahashi, M., and Sato, K.: General aspects of a T213L256 middle atmosphere general circulation model, J. Geophys. Res., 113, D12110, https://doi.org/10.1029/2008JD010026, 2008. a
Watanabe, S., Sato, K., Kawatani, Y., and Takahashi, M.: Vertical resolution dependence of gravity wave momentum flux simulated by an atmospheric general circulation model, Geosci. Model Dev., 8, 1637–1644, https://doi.org/10.5194/gmd816372015, 2015. a
Wedi, N. P. and Smolarkiewicz, P. K.: A framework for the testing global nonhydrostatic models, Q. J. Roy. Meteor. Soc., 135, 469–484, https://doi.org/10.1002/qj.377, 2009. a, b, c
White, A. A. and Bromley, R. A.: Dynamically consistent, quasihydrostatic equations for global models with a complete representation of the Coriolis force, Q. J. Roy. Meteor. Soc., 121, 399–418, 1995. a, b, c, d
White, A. A., Hoskins, B. J., Roulstone, I., and Staniforth, A.: Consistent approximate models of the global atmosphere: shallow, deep, hydrostatic, quasihydrostatic and nonhydrostatic, Q. J. Roy. Meteor. Soc., 131, 2081–2107, https://doi.org/10.1256/qj.04.49, 2005. a, b, c, d
Wilson, E. B.: Vector analysis, a textbook for the use of students of mathematics and physics, founded upon the lectures of J. Willard Gibbs, Yale University Press, New Haven, 1929. a, b
Wood, N. and Staniforth, A.: The deepatmosphere Euler equations with a massbased vertical coordinate, Q. J. Roy. Meteor. Soc., 129, 1289–1300, https://doi.org/10.1256/qj.02.153, 2003. a, b
Wood, N., Staniforth, A., White, A., Allen, T., Diamantakis, M., Gross, M., Melvin, T., Smith, C., Vosper, S., Zerroukat, M., and Thuburn, J.: An inherently massconserving semiimplicit semiLagrangian discretization of the deepatmosphere global nonhydrostatic equations, Q. J. Roy. Meteor. Soc., 140, 1505–1520, https://doi.org/10.1002/qj.2235, 2014. a, b, c
Zängl, G.: Extending the numerical stability limit of terrainfollowing coordinate models over steep slopes, Mon. Weather Rev., 140, 3722–3733, https://doi.org/10.1175/MWRD1200049.1, 2012. a
Zängl, G., Reinert, D., Rípodas, P., and Baldauf, M.: The ICON (ICOsahedral Nonhydrostatic) modelling framework of DWD and MPIM: description of the nonhydrostatic dynamical core, Q. J. Roy. Meteor. Soc., 141, 563–579, https://doi.org/10.1002/qj.2378, 2015. a, b, c, d, e, f, g, h, i
Zdunkowski, W. and Bott, A.: Dynamics of the atmosphere: a course in theoretical meteorology, Cambridge University Press, Cambridge, 2003. a, b, c, d, e, f
It would be more correct to refer to the horizontal unit vectors as e_{𝔱}(r_{j}) and e_{𝔫}(r_{j}), since they are defined only at the grid triangle edges j. In contrast to the zonal and meridional unit vectors of the geographical coordinate system, e_{λ}(r) and e_{φ}(r), e_{𝔱} and e_{𝔫} do not converge to a differentiable continuum for the number of triangles going to infinity, since the jumps in orientation from one edge to another adjacent edge remain, regardless of the horizontal grid resolution. Nevertheless, if computations require spatial derivatives of vectorial quantities, they can be performed in geographical coordinates and the result can be projected onto the e_{𝔱}(r_{j}) and e_{𝔫}(r_{j}). So we regard e_{𝔱} and e_{𝔫} as an “indirectly differentiable” continuum and write e_{𝔱}(r) and e_{𝔫}(r) in order to simplify matters.
After finishing the simulations presented in Sect. 3.2, we discovered that there is a bug in the reduction of heating rates calculated in the PSrad scheme necessary to account for chemical heating and energy loss due to airglow, which is fully applied at altitudes above 80 km. The reduction was about twice as large as intended, leading to too little solar heating in particular in a region close to 80 km. This will be fixed in future model versions.