Grid refinement in ICON v2.6.4
 Deutscher Wetterdienst, Offenbach am Main, Germany
 Deutscher Wetterdienst, Offenbach am Main, Germany
Correspondence: Günther Zängl (guenther.zaengl@dwd.de)
Hide author detailsCorrespondence: Günther Zängl (guenther.zaengl@dwd.de)
This article describes the implementation of grid refinement in the ICOsahedral Nonhydrostatic (ICON) modeling system. It basically follows the classical twoway nesting approach known from widely used mesoscale models like MM5 or WRF, but it differs in the way feedback from fine grids to coarser grids is applied. Moreover, the ICON implementation supports vertical nesting in the sense that the upper boundary of a nested domain may be lower than that of its parent domain. Compared to the wellestablished implementations on quadrilateral grids, new methods had to be developed for interpolating the lateral boundary conditions from the parent domain to the child domain(s). These are based on radial basis functions (RBFs) and partly apply direct reconstruction of the prognostic variables at the required grid points, whereas gradientbased extrapolation from parent to child grid points is used in other cases. The runtime flow control is written such that limitedarea domains can be processed identically to nested domains except for the lateral boundary data supply. To demonstrate the functionality and quality of the grid nesting in ICON, idealized tests based on the Jablonowski–Williamson test case (Jablonowski and Williamson, 2006) and the Schär mountain wave test case (Schär et al., 2002) are presented. The results show that the numerical disturbances induced at the nest boundaries are small enough to be negligible for real applications. This is confirmed by experiments closely following the configuration used for operational numerical weather prediction at DWD, which demonstrate that a regional refinement over Europe has a significant positive impact on the forecast quality in the Northern Hemisphere.
The ICON (ICOsahedral Nonhydrostatic) modeling framework is jointly developed by the German Weather Service (DWD), the Max Planck Institute for Meteorology (MPIM), the German Climate Computing Center (DKRZ) and the Karlsruhe Institute for Technology (KIT), targeting a unified global numerical weather prediction (NWP) and climate modeling system (GCM). The development work started in 2004 with basic research on the model grid and the numerical formulation of the dynamical core in a highly idealized shallowwater framework (Bonaventura and Ringler, 2005; Rípodas et al., 2009). The next step was the implementation of a hydrostatic dynamical core using the same model equations, time integration scheme and vertical discretization as the spectral transform ECHAM model (Wan et al., 2013). This work focused on investigating the effects of the different horizontal grid formulation (icosahedral vs. spherical harmonics) and provided the basis for a first version with full physics coupling. In parallel, nonhydrostatic dynamical cores were developed on hexagonal grids (Gassmann and Herzog, 2008; Gassmann, 2013) and triangular grids (Zängl et al., 2015). As discussed in Gassmann (2011) and Danilov (2012), Cgridtype triangular discretizations suffer from a spurious computational mode, which manifests itself in a rapidly oscillating checkerboard pattern in the horizontal divergence field, giving rise to numerical stability problems. To date, this problem is lacking a rigorous mathematical solution. However, Zängl et al. (2015) showed that the problem could largely be mitigated by a specific averaging of the velocity components entering the divergence operator. This pragmatic solution allowed retaining the triangular Cgrid discretization. Moreover, a triangular grid proved to be more suitable for implementing a twoway grid nesting capability than a hexagonal one. Triangular cells can be recursively partitioned into successively smaller triangles, which leads to a unique relationship between parent and child cells. For hexagons, however, this is not the case, as the majority of child cells are shared between two adjacent parent cells. To our knowledge, this is one aspect that makes twoway grid refinement approaches developed for hexagonal grids significantly more complex (e.g., Dubos and Kevlahan, 2013).
Despite impressive advances in computational power over the last decades, the application of global models with uniform, convectionpermitting resolution on weather or even climate timescales is still too costly to be performed on a regular basis. To date, highresolution limitedarea models (LAMs) serve as a costeffective alternative for exploring the mesoscale and microscale, and they will continue to serve as a working horse for both the numerical weather prediction (NWP) and climate community. Limitedarea models have proven successful, but they are known to have conceptional deficiencies, such as the potential illposedness of lateral boundary conditions (Davies, 2014), possible inconsistencies with the driving model in terms of the governing equations, numerical formulations or physical parameterizations, and the lack of regional to globalscale interactions (Warner et al., 1997). In order to mitigate these deficiencies, a number of methods have been considered to achieve locally enhanced resolution in global models.
The oldest approach dating back more than 40 years is the usage of stretched grids (Schmidt, 1977; Staniforth and Mitchell, 1978). The socalled Schmidt transformation allows for enhanced resolution in a particular region of interest by redistributing the grid points of an initially uniform model grid. This approach has proven successful in various NWP and climate applications. A review of GCM applications is provided, e.g., by FoxRabinovitz et al. (2008), and Goto et al. (2015) discuss a more recent NWP application using NICAM, which is based on a modified Schmidt transformation on an icosahedral grid (Tomita, 2008). Grid stretching obviates the need for lateral boundary conditions, and it allows for an immediate interaction between global and regional scales. Probably the largest drawback of this method is the fact that grid points can only be redistributed rather than created. Enhancing the resolution in one region of the globe implies a coarsening in another region. Besides the inevitable coarsening itself, the insufficient resolution of disturbances passing through the coarsened region may also negatively affect the simulation in the refined region.
More recently, global models using locally refined unstructured meshes have been developed. Unlike the grid stretching approach, they allow new grid cells to be added in regions of interest (static hrefinement). This approach is pursued, for example, by the Model for Prediction Across Scales (MPAS) (Skamarock et al., 2012) and the spectral element dynamical core of the Community Atmosphere Model (CAM) (Zarzycki et al., 2014). Albeit being much more flexible, this approach faces similar challenges as the grid stretching approach mentioned previously. In both approaches, the time step is restricted by the smallest cell in the domain unless specific measures like substepping are taken for individual cells or regions or unless horizontally implicit numerical methods are chosen. Moreover, care must be taken that the parameterizations are applicable on a wide range of scales and turn off gradually when the respective processes become resolved on the model grid (such as the convection parameterization in the gray zone). The development of scaleaware parameterizations poses major challenges and is an area of active research (Gross et al., 2018). However, notable progress in terms of scale awareness has been reported, e.g., for the CAM5 parameterization suite when compared to CAM4 (Gettelman et al., 2018).
The approach we are pursuing closely resembles traditional twoway nesting, as known from many regional mesoscale models such as MM5 (Grell et al., 1994) or WRF (Skamarock et al., 2019). Twoway nesting differs from the previous singlegrid approaches by the fact that multiple grids of different resolution are overlaid onto each other such that individual points on the globe are covered by more than one prognostic grid cell. Scale interaction can be ensured by feeding the nestedgrid solution back to the underlying parent grid at regular time intervals. Similar to LAMs, this method requires lateral boundary conditions to be specified for every nested domain, and it may suffer from spurious wave reflections at nest boundaries due to the abrupt jump in resolution. On the other hand, it allows many model settings to be chosen individually for each domain, which improves numerical efficiency and reduces the requirements concerning the parameterization suite. Distinct time steps can be chosen for each domain, allowing, e.g., a proportional reduction in refined domains in order to meet stability constraints. The parameterizations can be tuned individually for each domain and resolution or even switched off, relaxing the requirement of scale awareness to some degree. While twoway nesting has been successfully applied in limitedarea modeling for decades, it is still much less common in global modeling. Focusing on recent global models which are being used in research or operational forecasting, the authors are only aware of the twoway nesting approach implemented in the Geophysical Fluid Dynamics Laboratory (GFDL) FiniteVolume CubedSphere Dynamical Core (FV3) (Harris and Lin, 2013, 2014; Mouallem et al., 2022).
The purpose of this article is to describe the implementation of grid nesting in ICON. Compared to previous twoway nesting approaches, it is the first implementation on triangular grids and differs in the way feedback from fine grids to coarser grids is realized. In addition, the ICON implementation allows for vertical nesting in the sense that the upper boundary of a nested domain may be lower than that of the parent domain, which is not supported by most previous nesting implementations in other models. Without specific discussion, we note that a limitedarea mode is available in ICON as a byproduct of the grid nesting implementation, differing from nesting only in the way the lateral boundary conditions are provided. A detailed description the grid nesting implementation will be provided in Sect. 2, followed by application examples in Sect. 3. A brief summary will be given in Sect. 4.
2.1 Description of the basic design
The understanding of ICON's nesting implementation requires some basic knowledge of ICON's mathematical–physical design. We therefore start by summarizing key elements of the dynamical core, time stepping and physics–dynamics coupling scheme. ICON's dynamical core solves the fully compressible, nonhydrostatic Euler equations on the sphere using either the shallow or deep atmosphere formulation (Borchert et al., 2019). The spatial discretization is performed on an unstructured icosahedral–triangular ArakawaC grid in the horizontal and a terrainfollowing heightbased SLEVE coordinate (Leuenberger et al., 2010) with Lorenztype staggering in the vertical. As described in Zängl et al. (2015), the prognostic variables encompass the edgenormal horizontal wind speed v_{n}, vertical wind speed w, total air density ρ, virtual potential temperature θ_{v} and mass fractions q_{k} of various moisture quantities. See Fig. 1 in Wan et al. (2013) for the variable placement on the triangular grid. The twotimelevel predictor–corrector time integration scheme is fully explicit except for the terms describing the vertical propagation of sound waves, which are treated implicitly (Zängl et al., 2015). Hence, the permitted integration time step is comparatively small and constrained by the ratio of the speed of sound to the horizontal mesh size.
To optimize computational efficiency, different integration time steps are used for the dynamical core on the one hand and additional subgrid physical processes (and tracer transport) on the other hand. The time steps will be denoted by Δτ for the dynamical core and Δt for physics. The dynamical core is substepped with respect to physical processes, with a single physics time step consisting of five dynamics substeps by default ($\mathtt{nsubs}=\mathrm{\Delta}t/\mathrm{\Delta}\mathit{\tau}=\mathrm{5}$).
The physics–dynamics coupling scheme in ICON further distinguishes between fast physical processes, having a timescale shorter than or comparable to the time step Δt, and slow processes, which are considered to have a large timescale compared to Δt. Processes that fall into the category of fast are saturation adjustment, gridscale microphysics or turbulent diffusion, while examples of slow processes are convection and radiation. Fast processes are integrated with the previously defined time step Δt, whereas slow processes may be integrated with processspecific larger time steps that are an integer multiple of Δt. We therefore call the time step Δt the fast physics time step. As the tracer transport is performed at Δt as well, we apply the coupling of nested domains at Δt for both dynamics and tracer variables because this simplifies maintaining consistency with continuity (Gross et al., 2002). Essentially it is required to provide timeaveraged mass fluxes (averaged over the substeps) for tracer transport at nest lateral boundaries (see Sect. 2.2.1) and to ensure that the childtoparent feedback increments for partial densities ρq_{i} sum up to the feedback increments for total density ρ (see Sect. 2.2.2).
The static mesh refinement in ICON is accomplished using multiple individual grids, whereby one or more higherresolution (child) domains are overlaid on a coarser base (parent) domain. The base domain can be a regional or global domain. Model integration on the child domain is performed in addition to that on the underlying part of the parent domain; i.e., there is no “hole” in the parent domain where the child domain is located.
Each child domain has a defined parent domain providing lateral boundary conditions, but a parent domain can have several child domains. The child domains can be located in different geographical regions and can also serve as parent domains for further subdomains, but domains having the same parent are not allowed to share the same parent grid cells because this would lead to ambiguities in combination with twoway nesting. Nested domains may also be switched on or off during runtime. Conceptually, the number of nested domains is arbitrary and controlled by the grid files provided as input, but of course not all choices would make sense from a physical point of view. Each domain can be regarded as separate instances of the same model that are coupled to each other using the same numerical operators and filters, time integration scheme, and physics–dynamics coupling. If desired, however, different physical settings can be chosen individually for each domain. For example, radiation can be called more frequently on subdomains, or the convection scheme can be tuned differently or even switched off completely.
The refinement ratio between the parent domain and a child domain is fixed to a value of 2, as each parent triangle is split into four child triangles (see the right part of Fig. 1). While higher refinement ratios would technically be possible, we decided for this restriction because it reduces the risk of numerical artifacts (e.g., by partial wave reflection) along nest boundaries. Consistent with the refinement ratio of 2, the model integration time step Δt is multiplied by a factor of 0.5 for each additional nesting level. The coupling time step between successive nesting levels is the fast physics time step Δt.
The parent–child coupling can be either oneway or twoway. A mixture of oneway and twoway coupled domains is also possible. Twoway versus oneway coupling means that the prognostic variables on the child domain are transferred back to the coarser parent domain at regular time intervals using a dedicated feedback mechanism. As a result, the solution on the parent domain benefits from the higher resolution of the child domain. In the case of oneway nesting, the feedback is switched off.
To perform the coupling, we conceptually split any nested domain into three zones, which we call the boundary interpolation zone, the nudging zone and the feedback zone. In order to identify grid points belonging to these zones, cells, edges and vertices are indexed according to their distance from the boundary (see Appendix A1 for details). These zones along with the grid point indexing are depicted in Fig. 1, which shows part of the boundary region and inner region of a nested domain. Limitedarea domains are treated technically like oneway nested domains except for the fact that the boundary interpolation zone is filled with external input data rather than data interpolated from the parent domain.
The boundary interpolation zone is nonprognostic and is meant to store the boundary conditions that are necessary to solve the governing equations in the child domain. Boundary conditions are needed for the prognostic variables v_{n}, w, ρ, θ_{v} and q_{k}. By a dedicated boundary update mechanism (see Sect. 2.2.1 for details), both the prognostic variables and their time tendencies are interpolated from the parent to the child domain, and the boundary conditions are updated at every child time step. The boundary interpolation zone has a fixed width of four cell rows. This is motivated by the technical constraint that the boundary zone needs to match parent cell rows (i.e., an odd number of cell rows is not allowed), combined with the fact that two cell rows would not be sufficient to cover all stencil operations performed in the dynamical core. For example, the ∇^{4}diffusion operator (Zängl et al., 2015) requires information from three adjacent cell rows. We note that halving the size of the boundary interpolation zone would be possible by modifying some stencil operations near the lateral boundary, such as replacing ∇^{4} with ∇^{2}. However, this would have no impact on computational efficiency in practical applications with MPI domain decomposition (see Appendix A2).
The nudging zone, which is active in the case of oneway nesting only, serves to damp differences between the driving solution in the adjacent boundary interpolation zone and the prognostic solution in the child domain. Essentially, the prognostic model state of the child domain is relaxed (nudged) to the parent state following the traditional Newtonian relaxation approach described by Davies (1976). Details of the implementation are provided in Sect. 2.2.3. For twoway nesting, no nudging is applied, and the boundary interpolation zone borders the feedback zone.
In the feedback zone, the model state on the parent domain is relaxed towards the updated model state on the child domain at every fast physics time step Δt_{p} of the parent domain. We refer to this as relaxationtype feedback. As a result, the parent and child domains remain tightly coupled, and the solution on the parent domain benefits from the enhanced resolution of the child domain. Feedback is applied to the prognostic variables v_{n}, w, ρ and θ_{v} as well as to the mass fractions of water vapor q_{v}, cloud water q_{c} and cloud ice q_{i}. Childtoparent feedback has already been successfully applied in mesoscale models like MM5 (Grell et al., 1994) or WRF (Skamarock et al., 2019), as well as in global simulations based on a cubedsphere grid (Harris and Lin, 2013). However, our approach differs in the sense that the parent state is relaxed towards the child state with an adjustable timescale, rather than being overwritten by the child state. Compared to the conventional direct feedback, which is available as an option in ICON as well, the relaxation feedback has the advantage of generating fewer numerical disturbances near vertical nest interfaces, leading to slightly better forecast quality in NWP applications (without vertical nesting, the quality difference is small). In addition, the relaxation feedback requires no adjustment of the model orography at the parent grid level, which is more convenient if nested domains are turned on or off during runtime. See Sect. 2.2.2 for further details.
2.2 Parent–child coupling
2.2.1 Lateral boundary update: parent → child
The boundary update mechanism provides the child domain with uptodate lateral boundary conditions for the prognostic variables v_{n}, w, ρ, θ_{v} and q_{k}. In order to prevent parenttochild interpolated values of ρ from entering the solution of the mass continuity equation, the above set of variables is extended by the horizontal mass flux ρv_{n}. This will allow for parent–child mass flux consistency, as described below. For the subsequent description of the algorithm, let the model state on the parent and child domain be denoted by ${\mathcal{M}}_{\mathrm{p}}^{n}$ and ${\mathcal{M}}_{\mathrm{c}}^{n}$, respectively, where n specifies the time step index.
In general, the boundary update works as follows. Let ${\mathit{\psi}}_{\mathrm{p}}^{n}$ and ${\mathit{\psi}}_{\mathrm{p}}^{n+\mathrm{1}}$ represent a prognostic variable on the parent domain at time steps n and n+1, respectively. Once the variables on the parent domain ℳ_{p} have been updated from n to n+1, the time tendency
is diagnosed. Both the field ${\mathit{\psi}}_{\mathrm{p}}^{n}$ at time level n and the tendency $\frac{\partial {\mathit{\psi}}_{\mathrm{p}}}{\partial t}$ are then interpolated (downscaled) from the parent grid cells or edges to the corresponding cells or edges of the child domain's boundary zone (orange cells in Fig. 1). With ℐ_{p→c} denoting the interpolation operator, we get
The interpolated tendencies are needed in order to provide the lateral boundary conditions at the right time levels, since two integration steps are necessary on the child domain in order to reach the model state ${\mathcal{M}}_{\mathrm{c}}^{n+\mathrm{1}}$, with each step consisting of several (typically five) dynamics substeps. The temporal update is performed at each dynamics substep Δτ_{c} for the prognostic variables of the dynamical core and at each large step Δt_{c} for the tracer variables. As an example, the boundary conditions at the first dynamics substep of the first and second (fast physics) integration step on the child domain read ${\mathit{\psi}}_{\mathrm{c}}^{n}$ and ${\mathit{\psi}}_{\mathrm{c}}^{n}+\mathrm{0.5}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}{t}_{\mathrm{p}}\partial {\mathit{\psi}}_{\mathrm{c}}/\partial t$, respectively.
Concerning the interpolation operator ℐ_{p→c}, we distinguish between cellbased variables (i.e., scalars) and edgebased variables (v_{n} and ρv_{n}). For cellbased variables, a 2D horizontal gradient is reconstructed at the parent cell circumcenter by first computing edgenormal gradients at edge midpoints, followed by a ninepoint reconstruction of the 2D gradient at the cell center based on radial basis functions (Narcowich and Ward, 1994). The interpolated value at the jth child cell center is then calculated as
with ∇ψ_{p} denoting the horizontal gradient at the parent cell center and d(p,c_{j}) the distance vector between the parent and jth child cell center. The same operator is applied to cellbased tendencies.
To prevent excessive overshoots and undershoots of ${\mathit{\psi}}_{{c}_{j}}$ in the vicinity of strong gradients, a limiter for ∇ψ_{p} is implemented. It ensures that
on all four child points, where ψ_{p,𝚖𝚒𝚗} and ψ_{p,𝚖𝚊𝚡} denote the minimum and maximum of ψ_{p}, respectively, on the abovementioned reconstruction stencil including the local cell center, and β=1.05 is a tuning parameter. To minimize interpolation errors above steep orography, perturbations from the reference state (Zängl et al., 2015) rather than the full values are interpolated for the thermodynamic variables ρ and θ_{v}. In addition, the interpolation operator ℐ_{p→c} is also applied to the model orography in the boundary interpolation zone in the setup phase of the model (before calculating the vertical grid) in order to ensure consistency of the model grids even if different raw data sets have been used to generate the model orographies. To allow a smooth transition to the orography in the interior of the nested domain, a linear blending with a width of eight cell rows (as for the abovementioned nudging zone) is applied.
Regarding the interpolation of edgebased variables (i.e., the edgenormal vector components v_{n} and ρv_{n}), we distinguish between outer child edges coinciding with the edges of the parent cell and inner child edges (see Fig. 2a).
Edgenormal vector components at the inner child edges are reconstructed using a direct RBF reconstruction based upon the fivepoint stencil indicated by solid red dots in Fig. 2a. For a given inner child edge the stencil comprises the edges of the corresponding parent cell and the two edges of the neighboring parent cells that (approximately) share the orientation of the inner child edge.
For the outer child edges, a more sophisticated reconstruction is applied in order to ensure that the mass flux across a parent edge equals the sum of the mass fluxes across the corresponding child edges. We start with an RBF reconstruction of the 2D vector of the respective variable at the parent triangle vertices (blue triangles in Fig. 2b) using the six (five at the original vertices of the icosahedron, the pentagon points) edge points adjacent to a vertex (red dots).
The edgenormal vector component ϕ at the child edge is then computed as
with d(p,c_{e}) denoting the distance vector between the parent and child edge midpoints for a given parent edge and ∇_{t}ϕ_{p} denoting the gradient of the edgenormal vector component ϕ_{p} tangential to the parent edge. The latter is computed by projecting the reconstructed 2D vectors at the two vertices of an edge onto the edgenormal direction and taking the centered difference. Since by construction $\mathit{d}(p,{c}_{\mathrm{1}})=\mathit{d}(p,{c}_{\mathrm{2}})$ holds on the ICON grid, the abovementioned mass flux consistency is ensured. It is noted that attempts to use higherorder polynomial interpolation methods, which are the standard in mesoscale models with regular quadrilateral grids, were unsuccessful on the triangular ICON grid because the ensuing equation system leads to the inversion of nearly singular matrices.
Rather than interpolating v_{n} and its time tendency, only the time tendency is interpolated and then used to update v_{n} at child level at every dynamics time step. The wind field v_{n} itself is interpolated only once during the initialization of the child domain. This methodology has been chosen because the comparatively inaccurate interpolation to the interior child edges tends to induce smallscale noise in v_{n}. To suppress the remaining noise arising from the interpolation of the time tendency, a secondorder diffusion operator is applied in the inner half of the boundary interpolation zone on v_{n}, and the default fourthorder diffusion applied in the prognostic part of the model domain (Zängl et al., 2015) is enhanced in up to five grid rows adjacent to the interpolation zone by an amount exponentially decaying with distance from the boundary. For the secondorder diffusion, a coefficient of $\mathrm{0.005}\phantom{\rule{0.125em}{0ex}}{a}_{\mathrm{e}}/\mathrm{\Delta}t$ is applied, where a_{e} is the area represented by the current triangle edge, and the scaling with Δt means that the amount of diffusion is effectively independent from the time step. This proved to be sufficient to suppress the development of spurious disturbances even in the data assimilation cycle, which turned out to be the most critical application mode in this respect. For the other prognostic variables, no special filtering is applied near nest boundaries. In the case of oneway nesting, the secondorder velocity diffusion is extended into the nudging zone of the nested domain, replacing the enhanced fourthorder diffusion. More details on the nudging zone are given in Sect. 2.2.3.
For the horizontal mass flux ρv_{n}, the time average over the dynamic substeps is interpolated instead of time level n, as the timeaveraged mass flux is used by the tracer transport scheme in order to achieve consistency with continuity. Using the mass flux time tendency that is interpolated as well, the related time shift is corrected for when applying the boundary mass fluxes at the child level. In the nested domain, the interpolated mass fluxes valid for the current time step are then prescribed at the interface edges separating the boundary interpolation zone from the prognostic part of the nested domain (edges nr. 9 in Fig. 1). Due to the fluxform scheme used for solving the continuity equation (Zängl et al., 2015), this implies that the interpolated values of ρ do not enter into any prognostic computations in the dynamical core. They are needed, however, for flux limiter computations in the transport scheme. Moreover, no mass fluxes at interior child edges are used so that the nonconservative interpolation method used for those edges does not affect the model's conservation properties. For θ_{v} and the tracer variables q_{k}, the values at the edges are reconstructed in the usual manner following Eq. (20) in Zängl et al. (2015) and then multiplied with the interpolated mass fluxes before computing the flux divergences.
2.2.2 Feedback: child → parent
If twoway nesting is activated, the model state ${\mathcal{M}}_{\mathrm{p}}^{n+\mathrm{1}}$ on the parent domain is relaxed towards the updated model state ${\mathcal{M}}_{\mathrm{c}}^{n+\mathrm{1}}$ on the child domain at every parent fast physics time step Δt_{p}. This relaxationtype feedback is only applied to the prognostic variables v_{n}, w, θ_{v} and ρ plus specific humidity q_{v} and the specific contents of cloud water q_{c} and cloud ice q_{i}. Precipitating hydrometeors are excluded because recommended relaxation timescales (see below) are larger than their typical falling times. Surface variables are excluded as well because they can easily adjust during runtime, and the tile approach used in ICON's landsurface module would require a rather complicated algorithm to avoid inconsistencies.
Let ψ denote any of the above variables. Conceptually, the feedback mechanism is based on the following three steps.

Upscaling. The updated variable ${\mathit{\psi}}_{\mathrm{c}}^{n+\mathrm{1}}$ is interpolated (upscaled) from the child domain to the parent domain. The upscaling operators for cellbased and edgebased variables will be denoted by ℐ_{c→p} and ${\mathcal{I}}_{c\to p}^{e}$, respectively.

Difference computation. The difference between the parent domain variable ${\mathit{\psi}}_{\mathrm{p}}^{n+\mathrm{1}}$ and the upscaled child variable ${\mathcal{I}}_{c\to p}\left({\mathit{\psi}}_{\mathrm{c}}^{n+\mathrm{1}}\right)$ is computed.

Relaxation. The variable on the parent domain is relaxed towards the upscaled child domain variable by an increment that is proportional to the difference computed in step 2.
For edgebased normal velocity v_{n}, the arithmetic average of the two child edges lying on the parent edge is taken.
For cellbased variables the upscaling consists of a modified barycentric interpolation from the four child cells to the corresponding parent cell:
The weights α_{j} are derived from the following constraints in Eqs. (5)–(7). First of all, a necessary property for the interpolation operator is that it reproduces constant fields, i.e., the weights are normalized:
Moreover, the interpolation shall be linear: with the four child cell circumcenters x_{j} ($j=\mathrm{1},\mathrm{\dots},\mathrm{4}$) and x_{p} denoting the parent cell center, i.e., the interpolation target, we set
To motivate this constraint, consider the special case of equilateral triangles in which the center point of the inner child cell x_{1} coincides with the parent center such that the term (x_{1}−x_{p}) vanishes. Equation (6) then defines a barycentric interpolation within the triangle spanned by the mass points of the three outer child cells $\mathit{\{}{c}_{\mathrm{2}},{c}_{\mathrm{3}},{c}_{\mathrm{4}}\mathit{\}}$ (see Fig. 2a), where the weights $\mathit{\{}{\mathit{\alpha}}_{\mathrm{2}},{\mathit{\alpha}}_{\mathrm{3}},{\mathit{\alpha}}_{\mathrm{4}}\mathit{\}}$ represent the barycentric coordinates. Generally, the constraints in Eqs. (5) and (6) ensure reversibility in the sense that ℐ_{c→p} (ℐ_{p→c}) returns the original parent cell value irrespective of the reconstructed gradient.
Of course, the contribution of the point x_{1} closest to the interpolation target is of particular importance. Therefore, the underdetermined system of Eqs. (5) and (6) is closed with a final constraint which reads
where ${a}_{{\mathrm{c}}_{\mathrm{1}}}$ and a_{p} denote the inner child and parent cell areas, respectively. In other words, the inner child cell c_{1} containing the parent cell circumcenter is given a predefined weight corresponding to its fractional area coverage. This can be interpreted as a conservation constraint for the special case of a very localized signal at the mass point of the inner child cell.
In summary, this method can be regarded as a modified barycentric interpolation for the mass points $\mathit{\{}{\mathit{x}}_{\mathrm{2}},{\mathit{x}}_{\mathrm{3}},{\mathit{x}}_{\mathrm{4}}\mathit{\}}$, which accounts for x_{1} as an additional fourth source point. A more stringent barycentric interpolation would require an additional triangulation based on the child mass points.
We note that the cellbased operator ℐ_{c→p} is not strictly massconserving and that strict mass conservation would require some means of areaweighted aggregation from the child cells to the parent cells, which is available as an option. The difficulty with such methods on the ICON grid is related to the fact that the mass points lie in the circumcenter rather than the barycenter of the triangular cells. These points coincide on planar equilateral triangles, but they do not on general spherical triangles, the differences being largest in the vicinity of the pentagon points. Due to this fact, using an areaweighted aggregation from the child cells to the parent cells would map linear horizontal gradients on the child grid into a checkerboard noise pattern between upward and downwardoriented triangles on the parent grid.
Another difficulty that was encountered in the context of mass conservation is related to the fact that the density decreases roughly exponentially with height. In the presence of orography, the atmospheric mass resolved on the model grid therefore increases with decreasing mesh size, assuming the usual areaweighted aggregation of the orographic raw data to the model grid. Feeding back ρ is thus intrinsically nonconservative. To keep the related errors small and nonsystematic and to generally reduce the numerical errors over steep mountains, perturbations from the reference state are used for upscaling ρ and θ_{v} to the parent grid. A closer investigation of the related conservation errors revealed that the previously mentioned differences between bilinear and areaweighted averaging are negligible (with real orography) compared to the meshsizerelated conservation error.
When combining the abovementioned steps, the feedback mechanism for ρ can be cast into the following form:
Here ${\mathit{\rho}}_{\mathrm{p}}^{n+\mathrm{1}}$ denotes the parent cell density, which has already been updated by dynamics and physics. The superscript ∗ signifies the final solution including the feedback increment. Δt_{p} is the fast physics time step on the parent domain, and τ_{fb} is a userdefined relaxation timescale that has a default value of τ_{fb}=10 800 s and is independent of the relaxed field. The smaller the value of τ_{fb}, the faster the parent state is drawn towards the child state. The chosen default value is optimized for our typical NWP applications and aims at filtering smallscale transient features from the feedback, while fully capturing synopticscale features.
Finally note that the upscaled density includes the correction term Δρ_{corr}, which has been introduced to account for height differences between the child and parent cell circumcenters. At locations with mountainous orography, the heights at parent cell circumcenters can differ markedly from those at the corresponding child cells. Without an appropriate correction, the feedback process would introduce a noticeable bias in the parent domain's mass field. The correction term is given by
with the parent–child difference of the reference density field
and the potential temperature perturbation ${\mathit{\theta}}_{v,c}^{\prime \phantom{\rule{0.125em}{0ex}}n+\mathrm{1}}={\mathit{\theta}}_{v,c}^{n+\mathrm{1}}{\mathit{\theta}}_{v\phantom{\rule{0.125em}{0ex}}\mathrm{ref},c}$. The term Δρ_{ref,p} is a pure function of the parent–child height difference and can be viewed as a leadingorder correction term. As a further optimization, the empirical factor $(\mathrm{1.05}\mathrm{0.005}\phantom{\rule{0.125em}{0ex}}{\mathcal{I}}_{c\to p}({\mathit{\theta}}_{v}^{\prime \phantom{\rule{0.125em}{0ex}}n+\mathrm{1}}\left)\right)$ was added, which means that the correction is roughly proportional to the actual air density. We further note that a possibly more accurate and less ad hoc approach would require a conservative remapping step in the vertical, prior to the horizontal upscaling.
Care is required in order to achieve consistency with continuity. For this purpose, feedback is not implemented for the tracer mass fractions directly, but for partial densities. Building upon the implementation for ρ, we get
Mass fractions are rediagnosed thereafter:
By summing Eq. (9) over all partial densities, we recover Eq. (8) for the total density.
A very similar approach is used for θ_{v}. As for ρ, only the increment of θ_{v} is upscaled from the child domain to the parent domain and added to the parent reference profile θ_{v ref,p}.
The same approach is taken for w, but the full field is upscaled.
In the case of v_{n} some secondorder diffusion is added to the ensuing feedback increment in order to damp smallscale noise:
with the feedback increment
and the diffusion coefficient $K=\frac{\mathrm{1}}{\mathrm{12}}\frac{{a}_{\mathrm{p},\mathrm{e}}}{\mathrm{\Delta}{t}_{p}}$, where a_{p,e} is the area of the quadrilateral spanned by the vertices and cell centers adjacent to the parent's edge.
2.2.3 Lateral nudging
If the feedback is turned off, i.e., if oneway nesting is chosen, a nudging of the prognostic child grid variables towards the corresponding parent grid values is needed near the lateral nest boundaries in order to accommodate possible inconsistencies between the two grids, particularly near the outflow boundary. Nudging is performed every fast physics time step of the child domain Δt_{c} following Davies (1976), by adding a forcing term to the prognostic equations for v_{n}, ρ, θ_{v} and q_{v} of the form
with the nudging coefficient α_{nudge} and the term in brackets denoting the nudging increment. Because lateral boundaries are in general not straight lines on the unstructured ICON grid, attempts to make an explicit distinction between inflow and outflow boundaries (e.g., by prescribing v_{n} at inflow boundaries only) were not successful.
To compute the nudging increment, the child grid variables are first upscaled to the parent grid in the same way as for the feedback (Eqs. 3 and 4), followed by taking the differences between the parent grid variables and the upscaled child grid variables. The differences are then interpolated back to the child grid using the same methods as for the lateral boundary conditions (Eqs. 1 and 2). The nudging coefficient α_{nudge} decreases exponentially from the inner margin of the boundary interpolation zone towards the interior of the model domain and is defined as
with the maximum nudging coefficient A_{0}, the cell row index r, the nudging zone start index r_{0}=5 (see Fig. 1), the nudging zone width L and the efolding width μ; the latter two are defined in units of cell rows. The coefficients A_{0}, L and μ may be adjusted by the user, with the default values given by A_{0}=0.1, L=8 cell rows and μ=2 cell rows. A secondorder diffusion on v_{n} is used near the lateral nest boundaries in order to suppress smallscale noise. As opposed to Eq. (10) for feedback, the diffusion operator acts on the velocity field itself rather than the velocity increment.
2.3 Vertical nesting
The vertical nesting option allows setting model top heights individually for each domain, with the constraints that the child domain height is lower than or at most equal to the parent domain height and that the child domain extends into heights at which the coordinate surfaces are flat. This allows, for instance, a global domain extending into the mesosphere to be combined with a child domain that extends only up to the lower stratosphere, which can save significant computational resources. However, a vertical refinement in the sense that the vertical resolution in the child domain may differ from that in the parent domain is not implemented.
Vertical nesting requires appropriate boundary conditions for all prognostic variables to be specified at the vertical nest interface level, i.e., the uppermost half level of the nested domain. This is crucial in order to prevent vertically propagating sound and gravity waves from being spuriously reflected at the nest interface. In the following, boundary conditions are derived for v_{n}, w, θ_{v}, ρ and q_{k} as well as the vertical mass flux ρw. We note that the boundary condition for w is required for the vertically implicit sound wave solver in the dynamics, whereas ρw is needed to compute the vertical flux divergence terms in the prognostic equations for ρ, π and ρq_{k}.
Due to the constraints mentioned above, boundary conditions can be derived by horizontal parenttochild interpolation, without the need for any boundary interpolation zone extending vertically away from the upper nest boundary. For w, θ_{v}, ρ and ρw the full fields at the nest interface level are horizontally interpolated from the parent to the child grid using the same RBFbased interpolation method as for the lateral boundary conditions (Eq. 1). Rather than interpolating instantaneous values as for the lateral boundaries, w, θ_{v}, ρ and ρw are averaged over all dynamics substeps constituting a fast physics time step in order to filter oscillations related to vertically propagating sound waves. Hence, for $\mathit{\psi}\in \mathit{\{}w,{\mathit{\theta}}_{v},\mathit{\rho},\mathit{\rho}w\mathit{\}}$ the parenttochild interpolated field reads
with s denoting an individual dynamics substep and 𝚗𝚜𝚞𝚋𝚜 denoting the total number of substeps. We apply the same interpolation method to the corresponding time tendencies $\partial {\mathit{\psi}}_{\mathrm{p}}/\partial t$, which are estimated by taking the difference of the state variables at the substeps s=1 and s=𝚗𝚜𝚞𝚋𝚜. This enables us to perform a linear interpolation in time in order to provide the boundary conditions at approximately the right time levels for every dynamics substep on the child domain.
A slightly different approach is taken for v_{n}, which turned out to be beneficial in order to reduce the magnitude of the horizontal interpolation errors. The differences between the nest interface level and the next halflevel below (denoted as Δv_{n,p} in the following) are interpolated rather than the full field, again using the same methods as for the lateral boundary conditions. After interpolating Δv_{n,p} to the child domain, it is added to v_{n,c} at the second interface level ($k=\mathrm{3}/\mathrm{2}$) on the child domain in order to obtain the upper boundary condition, i.e.,
Since Δv_{n} is less strongly affected by sound waves, only an average between the first and last dynamics substep is taken prior to the interpolation. The temporal interpolation is neglected.
For the tracer variables we refrain from interpolating the partial mass fluxes (ρwq_{k})_{p} directly in order to ensure tracer and air mass consistency. Instead, we make use of the vertical mass flux boundary condition $(\stackrel{\mathrm{\u203e}}{\mathit{\rho}w}{)}_{\mathrm{c}}$ and multiply it with proper mass fractions. On the parent domain the required mass fractions are derived by taking the ratio of the vertical tracer mass flux at the nest interface level calculated in the vertical tracer transport scheme (ρwq_{k})_{p} and the available mass flux (ρw)_{p}. The mass fractions are then interpolated to the child domain using Eq. (1). Hence, the flux boundary condition for an arbitrary tracer field q_{k} reads
with $(\stackrel{\mathrm{\u203e}}{\mathit{\rho}w}{)}_{\mathrm{c}}$ computed via Eq. (11.) We note that due to the lack of q_{k} values above the nest upper boundary (lack of a boundary interpolation zone), the flux computation for scalars at the second interface level ($k=\mathrm{3}/\mathrm{2}$) is only stable for vertical Courant numbers ${C}_{k}=\left{w}_{k}\right\mathrm{\Delta}{t}_{c}/\mathrm{\Delta}{z}_{k}{}_{k=\mathrm{3}/\mathrm{2}}\le \mathrm{1}$, with Δz denoting the vertical layer thickness.
2.4 Recursive algorithm for multidomain setups
The previous sections have focused on the coupling of a single child domain. The nesting capability of ICON, however, is not limited to a single domain but supports multiple nests at the same level and multilevel nesting, as well as a combination of both. In the literature, multilevel nesting is also referred to as telescoping nesting (Mouallem et al., 2022). An example of multiple samelevel nesting will be provided in Sect. 3.1, while for multilevel nesting we refer to the ICON simulations in Weimer et al. (2021), wherein a threedomain threelevel setup has been used to investigate mountainwaveinduced polar stratospheric clouds.
The coupling of multiple samelevel nested domains with a parent domain is rather straightforward, as it only requires the singlenest coupling strategy (Sect. 2.2.1–2.2.2) to be applied sequentially for each nest. The coupling strategy for repeatedly nested domains is probably less obvious and will be described here for clarity and in order to complement already existing applications of this feature (e.g., Weimer et al., 2021).
Figure 3 displays a basic multilevel nesting example, wherein a global domain is combined with two successively (twoway) nested domains. The global domain is indicated at the bottom, and the nested domains are vertically staggered on top of it. The red and blue regions show the boundary interpolation zones and feedback zones of the individual domains, respectively. The integration time step on the global domain is denoted by Δt. It is automatically reduced by a factor of 2 when progressing to the next child grid level.
The whole processing sequence for integrating all domains from time step n to n+1 is shown in the flowchart at the lower left of Fig. 3. The domains are ordered topdown. Open and filled black dots show model states without and with feedback increments included, black arrows indicate time integration with the fast physics time step, and red and blue arrows indicate lateral boundary data interpolation and feedback, respectively.
The flow control of ICON's hierarchical nesting scheme is handled by a recursive subroutine that cascades from the global domain (or outer limitedarea domain) down to the deepest nesting level and calls the dynamical core and the physics parameterizations for each domain in basically the same way as for the global domain. The basic processing sequence is as follows.

A single integration step with Δt is executed on the global domain, resulting in an updated model state ${\mathcal{M}}_{\mathrm{p}}^{n+\mathrm{1}}$, as indicated by an open black circle in Fig. 3.

Boundary data are interpolated from the global domain to nest 1 (red arrow), followed by an integration step on nest 1 over the time interval $\mathrm{\Delta}t/\mathrm{2}$ and resulting in the model state ${\mathcal{M}}_{c\mathrm{1}}^{n+\mathrm{1}/\mathrm{2}}$.

As another nested domain exists within nest 1, boundary fields based on the model state ${\mathcal{M}}_{c\mathrm{1}}^{n+\mathrm{1}/\mathrm{2}}$ are interpolated to the second nested domain. Afterwards, the model is integrated on nest 2 over 2 times the time interval $\mathrm{\Delta}t/\mathrm{4}$, resulting in the model state ${\mathcal{M}}_{c\mathrm{2}}^{n+\mathrm{1}/\mathrm{2}}$.

Feedback is conducted from nest 2 back to nest 1 (blue arrow), resulting in an updated model state ${\mathcal{M}}_{c\mathrm{1}}^{n+\mathrm{1}/\mathrm{2}\ast}$ on nested domain 1 (black filled dot). Then, on the nested domain 1 the model is again integrated in time to reach model state ${\mathcal{M}}_{c\mathrm{1}}^{n+\mathrm{1}}$.

This is followed by a second lateral boundary data interpolation from nest 1 to nest 2 based on ${\mathcal{M}}_{c\mathrm{1}}^{n+\mathrm{1}}$. Nest 2 is integrated in time again to reach its state ${\mathcal{M}}_{c\mathrm{2}}^{n+\mathrm{1}}$.

As a final step, feedback is performed from nest 2 to nest 1, followed by feedback from nest 1 to the global domain.
We note that the presented coupling strategy is very similar to that in WRF or FV3, but differences exist in several details. For example, the childtoparent feedback in FV3 covers only temperature and the wind components (Mouallem et al., 2022), which avoids any impact on mass conservation but apparently necessitates executing the feedback before the physics call at the parent level in order to maintain numerical stability.
3.1 Jablonowski–Williamson test
To demonstrate the functionality of the grid nesting and to investigate the numerical errors related to the mesh size discontinuity along the nest boundary, we start with the Jablonowski and Williamson (2006) test (called the JW test hereafter) already used in Zängl et al. (2015) and many other studies to evaluate basic aspects of numerical accuracy. This test considers the formation of baroclinic waves in a geostrophically and hydrostatically balanced zonally symmetric basic state with a very strong Equatortopole temperature contrast. In its standard version, the waves are triggered by a small perturbation in the initial wind field imposed at 20^{∘} N, 40^{∘} E. The perturbation grows very slowly during the first few days but evolves into an explosive cyclogenesis between days 7 and 10 of the test case. Dropping the initial perturbation, which implies that the model ideally should maintain the initial state, allows for testing the numerical errors related, for instance, to grid irregularities.
Our first series of experiments considers the baroclinic wave test for a variety of configurations based on the model grids R2B4 (mesh size 160 km) and R2B5 (80 km). Experiment E1 uses the global R2B4 grid only, experiments E2–E4 use a global R2B4 grid plus nested R2B5 grids at various locations along the track of the baroclinic wave (see Table 1 for details), and E5 uses a global R2B5 grid for reference. To account for the fact that these resolutions are much coarser than in typical NWP applications, the relaxation timescale τ_{fb} is increased to 12 h. The results are summarized in Fig. 4, showing the surface pressure and the relative vorticity at 850 hPa. For the nested experiments E2–E4, the results are displayed for the coarse (R2B4) grid, implying that the impact of the nested domain(s) appears indirectly via the feedback mechanism described above. Selected results for the nested domain are added in Fig. 5, while combined results are displayed in Fig. 6, where the highestresolution data are used in any region of the plot.
As discussed in Zängl et al. (2015), the general behavior of ICON (and most other grid point models) in the JW test is that the baroclinic wave train exhibits a phase lag at coarse resolution. Moreover, the leading cyclone tends to be too weak at coarse resolution, whereas the trailing cyclone tends to be slightly too intense because grid irregularities act as an additional trigger for wavelike disturbances. In Fig. 4, these features become evident from comparing E1 (Fig. 4a) with E5 (Fig. 4i). Comparing the nested runs E2 and E3 (Fig. 4c, e) with E1 reveals that the western nested domain, which is shared by E2 and E3 and is crossed by the baroclinic wave during the early stage of its evolution (days 2–4), leads to a slight reduction of the phase lag and the intensity bias of the trailing cyclone. On the other hand, there is very little difference between E2 and E3, implying that the eastern nested domain entered by the leading cyclone at the beginning of its rapid intensification around day 7 has only a weak impact on the cyclone intensity. More pronounced differences appear with the large nested domain of E4, covering the baroclinic wave train during most of its life cycle. In the nested domain (see Fig. 6), the position and intensity of the cyclones are now very similar to E5. The relaxationbased feedback to the coarse domain shown in Fig. 4g shifts the cyclones to the right positions, but their intensity is weaker than in the reference experiment E5. The intensity difference is primarily related to the abovementioned long relaxation timescale of 12 h, combined with the fact that the middle and left cyclones propagate at a somewhat slower speed in the coarse domain. Reducing τ_{fb} gradually decreases the difference, but for values less than about 5Δt_{p} (i.e., 2 h), the solution in the nested domain starts to deviate from the E5 reference because numerical disturbances generated along the nest boundaries become nonnegligible (not shown). Note in this context that in typical NWP applications, the resolution in the global domain is always high enough that cyclones and fronts move at the same speed as in the nested domain(s) and possess about the same intensity (except for tropical cyclones), so the intensity bias encountered here is not of practical relevance. On the other hand, the standard relaxation timescale of 3 h filters transient wavelike motions that are not deterministically predictable anyway, thus avoiding a transfer of smallscale noise into the global domain.
The relative vorticity fields (right panels of Fig. 4) confirm that the nested domain in E2 (and E3) reduces the phase lag of the trailing cyclone as well as the intensity of an extra disturbance on the upwind side that does not appear in E5. However, the intensity differences between E5 and the other experiments appear much more pronounced than in the surface pressure field. This is primarily because the accuracy at which derivativebased quantities like vorticity can be calculated depends directly on the mesh size of the model grid. To obtain better comparability, Fig. 5 displays vorticity fields for the leading cyclone computed at the same mesh size (R2B5) for E3, E4 and E5. Compared to the reference result from E5 (Fig. 5c), the E3 vorticity field exhibits significant distortions and a somewhat reduced amplitude (Fig. 5a), reflecting the fact that the cyclone was only poorly resolved during part of its development stage. For E4 (Fig. 5b), the differences are much smaller and concentrate in the vicinity of the nest boundary at 130^{∘} W. Note that there are no evident disturbances near the transition between the lateral boundary interpolation zone and the prognostic computational domain of the nested domain at about 135^{∘} W. Further, there is no evidence of significant discontinuities along the parent–child domain interface, as can be deduced from the combined solution plot in Fig. 6.
To further examine the flow disturbances generated by the resolution jump at the nest boundary, Fig. 7 compares the surface pressure fields for the steadystate JW test after 9 d of integration for E1 and E2. This variant of the JW test accentuates the socalled grid imprinting errors arising from any irregularity in the model grid, in our case due to both the nonuniformity of the icosahedral grid and the resolution jump along the nest boundaries. As a side note we mention that previous nonnested ICON results presented in Lauritzen et al. (2010) were obtained with an early version of the hydrostatic dynamical core (Wan et al., 2013), which did not perform as well as the nonhydrostatic dynamical core presented here. The result for E1 (Fig. 7a) shows the wellknown regular wavenumberfive disturbance pattern characteristic for icosahedral grids (see, e.g., Jablonowski and Williamson, 2006; Lauritzen et al., 2010). With the presence of a nested domain, the disturbances become irregular and reach a somewhat stronger peak amplitude (Fig. 7b). The largest disturbances occur in the air mass that was initially located in the nest region because the nestinduced perturbations there have the longest time to grow. For smaller values of τ_{fb}, the amplitude of these disturbances gradually increases (not shown). Regarding their practical relevance, we note that the baroclinicity of the real atmosphere is generally much weaker than in the initial state of the JW test because it is continuously depleted by synopticscale disturbances. Moreover, the nestinduced perturbations decrease about proportionally to the grid imprinting in the global grid (see Fig. 2 in Zängl et al., 2015) when refining the model resolution. In DWD's operational NWP applications, we have never encountered noticeable artificial disturbances along nest boundaries.
3.2 Schär mountain test with vertical nesting
In order to examine the behavior of the upper boundary condition for vertically nested domains, we conducted the Schär mountain test case (Schär et al., 2002). Therein, a uniform flow of constant wind speed U and Brunt–Väisälä frequency N over idealized hilly terrain gives rise to a combination of smallscale, vertically decaying nonhydrostatic gravity waves and largerscale vertically propagating quasihydrostatic gravity waves. We expect this test case to challenge the upper boundary condition implementation, as the vertically propagating waves should pass through the domain interface without spurious reflections or any accumulation of noise.
Following Zängl et al. (2015), we performed a quasithreedimensional analogue of this test on a limitedarea domain, where the original 1D mountain profile is changed to a 2D ridgelike profile that decays towards zero in the crossflow direction as the lateral domain boundary is approached. In our earlier study, the limitedarea setup was chosen in favor of a small planet configuration (Klemp et al., 2015) because it allows for a better control on the upstream flow conditions in the presence of high mountains, and it is retained here for convenience. The ridgelike profile is given by
with the originally proposed parameter settings h_{m}=250 m, a=5000 m and λ=4000 m, as well as the ridge length scale β=10^{5} m. The wind speed is set to $U=\mathrm{10}\phantom{\rule{0.125em}{0ex}}{\mathrm{ms}}^{\mathrm{1}}$, and the Brunt–Väisälä frequency is given by $N=\mathrm{0.01}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$ for z<20 km and gradually increases to $N=\mathrm{0.03}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$ above. Note that the latter setting differs from the constant N used by Schär et al. (2002) in order to allow for the high model top desired for our test configuration (see below). This does not allow for comparing against an analytic solution, but this is not needed here because model results serving this purpose have already been provided by Zängl et al. (2015).
We have performed two types of simulations: one with and one without a vertically nested domain. The reference simulation without a vertical nest was performed on a single $\mathrm{3}{}^{\circ}\times \mathrm{3}{}^{\circ}$ limitedarea domain centered at the Equator with a horizontal mesh size of Δx≈620 m (R2B12), a vertically stretched grid with 50 vertical layers and a 40 km model top. To avoid the reflection of gravity waves at the model top, a Rayleigh damping layer acting on vertical velocity w was applied above z=22 km (i.e., encompassing the uppermost 10 vertical layers). The simulation was run for 12 h until an approximately steady state was reached. The steadystate solution for w is depicted in Fig. 8a. While the wave structure at low levels closely matches previously published results for configurations with constant N (e.g., Schär et al., 2002; Skamarock et al., 2012; Zängl et al., 2015) partial wave reflections probably related to the stability change starting at z=20 km become apparent in the upper part of the domain.
For the nested simulation, a $\mathrm{2.4}{}^{\circ}\times \mathrm{2.4}{}^{\circ}$ R2B12 domain with 39 vertical layers and a vertical interface at 20 km is twoway nested into a $\mathrm{3}{}^{\circ}\times \mathrm{3}{}^{\circ}$ R2B11 parent domain (Δx≈1.2 km) with 50 vertical layers and a 40 km model top. There is no Rayleigh damping active in the nested domain. Note that the vertical layer distribution and the horizontal mesh size of the nested domain exactly match those of the reference simulation. The only significant difference lies in the much lower but open domain top. Steadystate solutions for the nested domain are depicted in Fig. 8b and c. They differ in terms of the relaxation timescale τ_{fb} for childtoparent feedback, with values of τ_{fb}=10 800 and τ_{fb}=900 s for panels (b) and (c), respectively.
In general, the results for the nested domain and the reference result are fairly close to each other. There is no indication of substantial wave energy reflection or noise accumulation along the nest interface level. Deviations from the reference result are largely confined to the uppermost quasihydrostatic wave crest and trough and to the leeward propagating wave signal. This holds not only for the steadystate solution, but also for the spinup phase (not shown). The reason for the deviations is twofold: firstly, the computation of the boundary conditions at the nest interface level inevitably goes along with spatiotemporal interpolation errors. Second, and more importantly, the solutions on the parent and child domains are slightly different, implying that the vertical interface condition derived from the parent domain cannot exactly match the solution on the child domain. These differences primarily originate from the differences in the mesh size, but they additionally depend on the feedback timescale. Reducing the feedback timescale strengthens the domain coupling and reduces the parent–child differences, which in turn improves the vertical nest interface conditions. This can be seen by comparing Fig. 8b (τ_{fb}=10 800 s) with Fig. 8c (τ_{fb}=900 s), which shows that a shorter feedback timescale further reduces the difference to the nonnested reference (Fig. 8a). The corresponding solution on the R2B11 parent domain (with τ_{fb}=900 s) is depicted in Fig. 9. Due to its 2 times coarser mesh size it is lacking some of the leewardpropagating wave signals and shows a slightly reduced amplitude of the quasihydrostatic wave. The attenuated, though visible, leewardpropagating wave signal in Fig. 9 results from the childtoparent feedback and does not exist in an R2B11 simulation without a nest (not shown).
From this test it can be concluded that, even without any interpolation error, the boundary conditions at the nest interface will never match perfectly due to the resolutioninduced differences between the parent and child model states. On the other hand, this test has shown a small but noticeable positive impact of the childtoparent feedback mechanism on the quality of the nest interface conditions, which improve with decreasing feedback timescale.
3.3 Operational NWP applications
In an operational context, one ideally expects the beneficial impact on forecast quality of the regionally refined resolution to be transferred to the global domain in the nest overlap area and subsequently propagate downstream, which is usually eastward in the extratropics. This implicitly assumes that the scores used to quantify forecast quality do improve with increasing model resolution, which, according to our experience, is the case for mesh sizes coarser than about 10 km but not necessarily in the convective gray zone. In the operational global forecasting system of DWD, the deterministic part has a global horizontal mesh size of 13 km and a twoway nested domain referred to as “ICONEU” with 6.5 km mesh size covering Europe and some adjacent regions (24.5^{∘} W–63.5^{∘} E, 29–71^{∘} N; i.e., domain 2 in Fig. 1). For this resolution range, the impact of the nested domain on the global forecast quality tends to be rather small. A much clearer impact is found for the corresponding ensemble prediction system (EPS), which (at the time of writing this paper) uses mesh sizes of 40 and 20 km. To demonstrate the benefit of twoway nesting on NWP quality, we therefore use the EPS configuration of ICON, but the subsequent experiments are executed as deterministic forecasts for simplicity and easier reproducibility outside the operational environment of DWD. Moreover, our experiments are initialized with interpolated operational IFS analyses available from the European Centre for MediumRange Weather Forecasts (ECMWF), which has the conceptual advantages that none of our experiments may benefit from being initialized from its “own” assimilation cycle and we can verify our forecast results against independent (i.e., IFS) analyses.
The current operational configuration at DWD uses 90 vertical levels with a model top at 75 km and a vertically nested EU domain with 60 levels and a vertical interface to the global domain at about 23 km. The forecast lead time is 180 h in the global domain and 120 h in the nested one. These settings are also used for our subsequent experiments unless stated otherwise. Specifically, we consider our nested EPS configuration, henceforth denoted as R2B6N7, a global 40 km mesh without a nest (R2B6) and a global 20 km mesh without a nest (R2B7). All experiment suites are conducted for January 2021, starting from the 00:00 UTC IFS analysis for each day of the month. The physics configuration equals the operational status of autumn 2021 as far as possible, notable exceptions being that no ensemble perturbations are used and that some tuning options relying upon coupling with the operational data assimilation scheme are turned off (see Reinert et al., 2021, for further details on the ensemble perturbations). Besides the abovementioned verification against IFS analyses, we use DWD's operational verification against SYNOP (surface) stations and radiosonde ascents (TEMP), following the WMO standard in all cases (WMO, 2019).
To exemplarily demonstrate the resolution dependence of the forecast quality and the related nest impact, Fig. 10 shows the root mean square error (RMSE) against IFS analyses for 500 hPa geopotential and vector wind in the Northern Hemisphere. It is clearly evident that the errors of R2B7 are smaller than for R2B6 during the whole forecast range. The nested setup R2B6N7 exhibits slightly smaller errors than R2B6 and is, as may be expected due to the size of the refined domain, closer to R2B6 than R2B7. A closer look at Europe is given in Fig. 11, showing relative differences to R2B6 for easier readability and including an additional test in which the nest stays active until 180 h (R2B6N7180h). It is seen that during the first 3 forecast days, the quality improvement related to the nest feedback in R2B6N7 is comparable to what is obtained with a global 20 km mesh (R2B7). Afterwards, the improvement in R2B6N7 starts to decrease as the advection from the coarse (global) domain into the nesting region becomes more and more relevant. After 5 d, the errors in R2B6N7 become temporarily a bit larger than in R2B6 over Europe, and we infer from the negligible differences between R2B6N7 and R2B6N7180h that this is unrelated to the usual termination of the nest at 120 h.
More detailed information on the nest impact is provided by the SYNOP and TEMP verification, which includes an assessment of statistical significance at the 95 % level. The SYNOP verification (Fig. 12) shows a pronounced beneficial impact in the nest overlap region (EUNEST) during the first 5 forecast days, which is largest for surface pressure (PS) with an improvement of almost 10 % on day 2 and still substantial for 10 m wind speed (FF), 2 m temperature (T2M) and 2 m humidity (RH2M). After the termination of the nested domain, the differences become small and insignificant. When the nest runtime is extended to 180 h, the improvements continue to be present for T2M and RH2M but are negligible for PS and FF (not shown). The results for Asia show that the improvements related to the EUNEST propagate downstream with a delay of about 3 d. They are a bit smaller than in the nest region but still partly reach the 95 % significance level, indicating an obvious benefit of our twoway nesting methodology. A more detailed view of the spatiotemporal evolution of the forecast skill improvement related to the nest feedback is provided in Fig. 13, showing a stationwise verification of surface pressure against all SYNOP stations used in the operational data assimilation cycle at DWD. On forecast days 1 and 2, the surface pressure RMSE is almost entirely smaller for R2B6N7 than R2B6, exceptions being restricted to a few individual stations. Later on, the spatial variability increases, with the first coherent regions with degradations appearing after about 3 d. Nevertheless, the impact of the nest stays predominantly positive, as seen in Fig. 12, and the eastward propagation of the average skill improvement is confirmed.
The TEMP verification (Fig. 14) additionally shows that the positive signal extends throughout the troposphere in both Europe and Asia. In the lower stratosphere, a slight degradation can be seen over Europe for wind speed and direction as well as temperature, which might be suspected to be related to the vicinity of the vertical nest interface at about 30 hPa. However, an additional sensitivity test with the nest extending up to 75 km (not shown) revealed that this is not the case. In this test, the degradation actually extends further up into the middle stratosphere, indicating that it is probably a doublepenalty effect related to resolving a larger part of the gravity wave spectrum, an issue we generally observe when increasing the horizontal model resolution (but the discussion of which is beyond the scope of this paper). On average over the Northern Hemisphere, the improvements related to the EUNEST are still significant for all variables, the magnitude in the TEMP verification being largest for wind speed and direction. We finally mention that the improvement of upperair geopotential is smaller than for PS, which appears to be a contradiction but is related to the fact that the density of radiosonde stations is more homogeneous than that of the surface stations, giving Europe an overly large weight in the SYNOP verification. Nevertheless, this does not affect our conclusion that the regional refinement over Europe has a clear and significant benefit for the forecast quality in the Northern Hemisphere, demonstrating that the twoway nesting capability in ICON fully serves its purpose.
Further analysis has been undertaken to investigate potential numerical issues related to the nesting, such as possible artifacts along the lateral boundaries of the nested domain and the loss of exact mass conservation mentioned in Sect. 2. Regarding boundary artifacts, precipitation is known to be one of the most sensitive variables because unphysical flow convergences or divergences tend to induce lineshaped structures along the boundaries that are readily recognized as unrealistic. To consider this aspect, we selected a forecast run in which substantial precipitation amounts occurred along the lateral boundaries (15 January 2021). The result after 5 forecast days is displayed in Fig. 15 for the nested domain and the immediate surroundings. The boundary interpolation zone, in which no prognostic precipitation is calculated in the nested domain, has been filled with interpolated values from the global domain. Although it is visible at a few spots that the interpolated precipitation along the boundaries is smoother than the prognostic one in the interior of the domain (as expected from the coarser model resolution at which it has been calculated), there is no evident discontinuity pointing to the presence of numerical problems.
To examine the impact of the twoway nesting on mass conservation, an additional set of forecast experiments has been conducted in which the lead time was extended to 30 d and the nested domain stays active until the end. Results for 18 arbitrarily selected initial dates are summarized in Fig. 16, showing the relative mass change in the global domain due to the presence of the twoway nest (the accuracy is close to machine precision otherwise). While there appears to be a slight systematic loss of mass during the first forecast days, which appears to be related to spinup effects, there is no systematic trend visible later on. The relative errors are on the order of 10^{−6}, which is small enough to be irrelevant for NWP purposes, and the absence of a trend indicates that the remaining conservation accuracy would even be sufficient on climate timescales (although this has not been further examined). Additional experiments indicated that the typical mass change during the first 2–3 d is smaller and not systematically negative when starting from operational analyses (not shown). Thus, no further efforts have been undertaken to investigate the apparent spinup effect with interpolated initial conditions.
This article provides a technical description of the grid nesting implementation in the ICOsahedral Nonhydrostatic (ICON) modeling system. The available options comprise oneway and twoway nesting, with one or more domains per nesting level, and vertical nesting in the sense that the upper boundary of a nested domain may be lower than that of its parent domain. In addition, a limitedarea mode is available as a byproduct of the grid nesting implementation, which differs from oneway nesting only in the way the lateral boundary conditions are provided. The modelinternal flow control basically follows the nesting approach known from mesoscale models like MM5 (Grell et al., 1994) or WRF (Skamarock et al., 2019), but the feedback applies a Newtonian relaxation rather than a direct replacement of the prognostic variables in the overlapping part of the parent domain. The relaxationtype feedback was found to produce slightly better forecast quality in NWP applications and obviates the need to adjust the model orography at the parent grid level. The main innovation that was needed compared to the wellestablished nesting implementations on quadrilateral grids was the development of appropriate interpolation algorithms for the lateral boundary conditions of the nested grids, as the usual higherorder polynomial interpolation leads to the inversion of a nearly singular matrix on the triangular ICON grid. It was found that radial basis functions constitute an adequate method. They are used either for reconstructing spatial gradients on the parent grid, followed by a linear extrapolation of the prognostic variables from the parent grid to the child grid, or for directly reconstructing the variable values on the child grid points.
To demonstrate the functionality and quality of the grid nesting in ICON, idealized tests restricted to the dry dynamical core are presented as are real forecast experiments using a model configuration very close to that used for operational NWP at DWD. The horizontal grid nesting is addressed by two variants of the Jablonowski–Williamson test case (Jablonowski and Williamson, 2006): one considering the evolution of a baroclinic wave train triggered by a disturbance in the initial condition and one ideally remaining in steady state in the absence of grid irregularities. Choosing the global mesh size such that the baroclinic wave is poorly resolved, the former test shows that a refined grid passed by the wave during the early stage of its development has some beneficial impact on the quality of the simulation, delivering a result lying between nonnested reference runs with the respective coarse and fine mesh sizes. Numerical disturbances related to the resolution jump along the nest boundaries are clearly smaller than the benefit of the regionally refined resolution. The steadystate variant of the test nevertheless shows that this kind of disturbance does exist and is somewhat larger than the disturbances induced by the irregularities of the uniform global icosahedral grid. The accuracy of the vertical nesting is examined with the Schär mountain test case (Schär et al., 2002), which considers a steadystate multiscale orographic gravity wave composed of vertically propagating and vertically decaying wave components. The results show that the wave reflections generated at the vertical nest interface are small enough to be irrelevant for practical applications. The findings from these idealized tests are corroborated by fullphysics NWP experiments using the operational DWD domain configuration with a vertically nested grid over Europe and adjacent regions. Provided that the mesh sizes are chosen to lie in a range in which the forecast quality has a pronounced resolution dependence, these tests demonstrate that the feedback from the regionally refined domain to the global one has a significant beneficial impact on a variety of standard NWP scores, which also becomes evident downstream (over Asia) with a delay of about 3 d. We are thus able to conclude that the ICON grid nesting, which has been used for operational weather forecasting at DWD since July 2015, has a practical and measurable benefit for our forecast quality.
A1 Grid point indexing
In order to identify grid points belonging to certain regions of a domain and to control the feedback between parent and child domains, dedicated integer fields named 𝚛𝚎𝚏𝚒𝚗_𝚌𝚝𝚛𝚕 exist for cells, edges and vertices (see Fig. 1). They are provided by ICON's grid generator. Along lateral nest boundaries, these fields contain positive numbers. For cells, they indicate the shortest distance to the boundary in units of cell rows. For example, a value of 1 indicates the outermost cells. A similar counting is applied for vertices. For edges, however, the counting proceeds twice as fast; i.e., the edgebased counter is the sum of two neighboring cell counters. This finegranuled counting enables the explicit specification of cell rows or edges up to which individual numerical operators (such as divergence, gradient or diffusion operators) are evaluated. By default, the 12 outermost cell rows (encompassing the boundary interpolation and feedback zones) are flagged by positive numbers. For the example in Fig. 1, ICON's grid generator had been configured to flag the 14 outermost cell rows, which then allows for an extended nudging zone of 10 cell rows at most.
In regions with an overlapping child domain the counters for cells, edges and vertices are filled with negative numbers (see lightblue area in Fig. 1). The counters start from −1 at the location of the child domain's outer boundary and descend towards the interior of the nest overlap region. The procedure is equivalent to the procedure along the lateral boundary, with the exception that it already stops after three cell rows. All remaining interior points are flagged with a value of −4 for cells and vertices and −8 for edges. This flagging of cells overlapping a child domain allows easy access to all the points involved in the childtoparent feedback process described in Sect. 2.2.2.
A2 Distributedmemory parallelization
Several measures are taken in order to optimize the computational efficiency of the nesting implementation.
Model grid points in ICON are stored in (blocked) 1D vectors. In order to achieve efficient runtime flow control without the need for masking operations within loops, the grid points are ordered by their distance from the lateral boundary by making use of the 𝚛𝚎𝚏𝚒𝚗_𝚌𝚝𝚛𝚕 fields described in Appendix A1. Hence, grid points lying at or near the lateral boundary of a nested domain are shifted to the beginning of the index vector. This allows excluding boundary points from prognostic computations accessing nonexisting neighbor points without masking operations. In the present implementation, the four outer cell rows constituting the boundary interpolation zone (Fig. 1) and the adjacent fifth one participate in the reordering. The remaining grid points follow in a nonordered fashion, while the 𝚛𝚎𝚏𝚒𝚗_𝚌𝚝𝚛𝚕 indexing proceeds with the remaining grid points in the nudging zone, followed by the prognostic grid points not overlapping a child domain (𝚛𝚎𝚏𝚒𝚗_𝚌𝚝𝚛𝚕=0), and by the overlapping prognostic points (𝚛𝚎𝚏𝚒𝚗_𝚌𝚝𝚛𝚕<0) (see Fig. 1). In the presence of an MPI domain decomposition, the related halo points are shifted to the very end of the index vector. Moreover, the two outermost cell rows are counted with zero weight when computing the domain decomposition, implying that the benefit of these grid points doing very little computation is not lost by load imbalance.
Regarding distributedmemory (MPI) parallelization, the default strategy adopted in ICON is to distribute all model domains among all compute processors. As this implies that child grid points are in general owned by a different processor than the corresponding parent grid point, an intermediate layer having the mesh size of the parent grid but the domain decomposition of the child grid is inserted in order to accommodate the data exchange required for boundary interpolation and feedback.
To reduce the amount of MPI communication for complex nested configurations, multiple nested domains at the same nesting level can be merged into one logical domain, which is then not geometrically contiguous. This needs to be specified by the user during the grid generation process by indicating a list of individual domains that are supposed to be merged. The lateral boundary points belonging to all components of the merged domain are then collected at the beginning of the index vector. For all prognostic calculations, the multiple domains are treated as a single logical entity, and just the output files may be split according to the geometrically contiguous basic domains. As oneway and twoway nesting cannot be mixed within one logical domain, there may still need to be two logical domains on a given nest level.
To further optimize the amount of MPI communication, a socalled processor splitting is available that allows for executing several nested domains concurrently on processor subsets whose size can be determined by the user in order to minimize the ensuing load imbalance. Unlike domain merging, this allows parallel execution of oneway and twoway nested domains. This option is currently restricted to the step from the global domain to the first nesting level in order to keep the technical complexity at a manageable level.
The ICON release version icon2.6.4 is freely available under a personal noncommercial research license. Information on the license and instructions for downloading the code can be found at https://code.mpimet.mpg.de/projects/iconpublic/wiki/Instructions_to_obtain_the_ICON_model_code_with_a_personal_noncommercial_research_license (MPIM, 2019). By downloading the code the user accepts the license agreement. All primary data and scripts which are necessary to validate the research findings of Sect. 3.1 and 3.2 are freely available for download from the Open Research Data Repository of the Max Planck Society (Edmond) under https://doi.org/10.17617/3.NOC2AE (Zängl et al., 2022). Data and scripts for the NWP applications of Sect. 3.3 are not included due to the huge size of the data set (about 3.3 TB) and the fact that installing and running the verification toolchain outside DWD's computer systems is beyond the authors' expertise. Access to these data will be granted upon request (contact: guenther.zaengl@dwd.de).
All authors contributed to the writing of this paper. GZ conducted the initial nesting implementation and ran the experiments of Sect. 3.1 and Sect. 3.3. DR worked on the vertical nesting and ran the experiments of Sect. 3.2. FP developed the grid generator and worked on parallelization aspects and performance optimizations. Some elements of the description in Sect. 2 have already been published in the ICON user tutorial (Prill et al., 2020).
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to thank the entire ICON team for their excellent work, without which this would not have been possible. Our additional thanks go to Felix Fundel (DWD), who developed the verification toolchain used in this work. The authors also appreciate the constructive comments made by the two anonymous reviewers, which led to a significant improvement of the paper.
This paper was edited by Ludovic Räss and reviewed by two anonymous referees.
Bonaventura, L. and Ringler, T.: Analysis of discrete shallow water models on geodesic Delaunay grids with Ctype staggering, Mon. Weather Rev., 133, 2351–2373, https://doi.org/10.1175/MWR2986.1, 2005. a
Borchert, S., Zhou, G., Baldauf, M., Schmidt, H., Zängl, G., and Reinert, D.: The upperatmosphere extension of the ICON general circulation model (version: uaicon1.0), Geosci. Model Dev., 12, 3541–3569, https://doi.org/10.5194/gmd1235412019, 2019. a
Danilov, S.: On utility of triangular Cgrid type discretization for numerical modeling of largescale ocean flows, Ocean Dynam., 60, 1361–1369, https://doi.org/10.1007/s1023601003396, 2012. a
Davies, H.: A lateral boundary formulation for multilevel prediction models, Q. J. Roy. Meteor. Soc., 102, 405–418, https://doi.org/10.1002/qj.49710243210, 1976. a, b
Davies, T.: Lateral boundary conditions for limited area models, Q. J. Roy. Meteor. Soc., 140, 185–196, https://doi.org/10.1002/qj.2127, 2014. a
Dubos, T. and Kevlahan, N. K.R.: A conservative adaptive wavelet method for the shallowwater equations on staggered grids, Q. J. Roy. Meteor. Soc., 139, 1997–2020, https://doi.org/10.1002/qj.2097, 2013. a
FoxRabinovitz, M., Cote, J., Dugas, B., Deque, M., McGregor, J. L., and Belochitski, A.: Stretchedgrid Model Intercomparison Project: decadal regional climate simulations with enhanced variable and uniformresolution GCMs, Meteorol. Atmos. Phys., 100, 159–178, https://doi.org/10.1007/s007030080301z, 2008. 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.: A global hexagonal Cgrid nonhydrostatic dynamical core (ICONIAP) designed for energetic consistency, Q. J. Roy. Meteor. Soc., 139, 152–175, https://doi.org/10.1002/qj.1960, 2013. a
Gassmann, A. and Herzog, H.J.: Towards a consistent numerical compressible nonhydrostatic model using generalized Hamiltonian tools, Q. J. Roy. Meteor. Soc., 134, 1597–1613, https://doi.org/10.1002/qj.297, 2008. a
Gettelman, A., Callaghan, P., Larson, V. E., Zarzycki, C. M., Bacmeister, J. T., Lauritzen, P. H., Bogenschutz, P. A., and Neale, R. B.: Regional Climate Simulations with the Community Earth System Model, J. Adv. Model Earth Sy., 10, 1245–1265, https://doi.org/10.1002/2017MS001227, 2018. a
Goto, D., Dai, T., Satoh, M., Tomita, H., Uchida, J., Misawa, S., Inoue, T., Tsuruta, H., Ueda, K., Ng, C. F. S., Takami, A., Sugimoto, N., Shimizu, A., Ohara, T., and Nakajima, T.: Application of a global nonhydrostatic model with a stretchedgrid system to regional aerosol simulations around Japan, Geosci. Model Dev., 8, 235–259, https://doi.org/10.5194/gmd82352015, 2015. a
Grell, G. A., Dudhia, J., and Stauffer, D.: A description of the fifthgeneration Penn State/NCAR Mesoscale Model (MM5), Tech. Rep. NCAR/TN398+STR, University Corporation for Atmospheric Research, https://doi.org/10.5065/D60Z716B, 1994. a, b, c
Gross, E. S., Bonaventura, L., and Rosatti, G.: Consistency with continuity in conservative advection schemes for freesurface models, Int. J. Numer. Meth. Fl., 38, 307–327, https://doi.org/10.1002/fld.222, 2002. a
Gross, M., Wan, H., Rasch, P. J., Caldwell, P. M., Williamson, D. L., Klocke, D., Jablonowski, C., Thatcher, D. R., Wood, N., Cullen, M., Beare, B., Willett, M., Lemarié, F., Blayo, E., Malardel, S., Termonia, P., Gassmann, A., Lauritzen, P. H., Johansen, H., Zarzycki, C. M., Sakaguchi, K., and Leung, R.: Physics–Dynamics Coupling in Weather, Climate, and Earth System Models: Challenges and Recent Progress, Mon. Weather Rev., 146, 3505–3544, https://doi.org/10.1175/MWRD170345.1, 2018. a
Harris, L. M. and Lin, S.J.: A TwoWay Nested GlobalRegional Dynamical Core on the CubedSphere Grid, Mon. Weather Rev., 141, 283–306, https://doi.org/10.1175/MWRD1100201.1, 2013. a, b
Harris, L. M. and Lin, S.J.: GlobaltoRegional Nested Grid Climate Simulations in the GFDL High Resolution Atmospheric Model, J. Climate, 27, 4890–4910, https://doi.org/10.1175/JCLID1300596.1, 2014. a
Jablonowski, C. and Williamson, D.: 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, c, d
Klemp, J. B., Skamarock, W. C., and Park, S.H.: Idealized global nonhydrostatic atmospheric test cases on a reducedradius sphere, J. Adv. Model. Earth Sy., 7, 1155–1177, https://doi.org/10.1002/2015MS000435, 2015. a
Lauritzen, P. H., Jablonowski, C., Taylor, M. A., and Nair, R. D.: Rotated Versions of the Jablonowski SteadyState and Baroclinic Wave Test Cases: A Dynamical Core Intercomparison, J. Adv. Model Earth Sy., 2, 15, https://doi.org/10.3894/JAMES.2010.2.15, 2010. a, b
Leuenberger, D., Koller, M., Fuhrer, O., and Schär, C.: A generalization of the SLEVE vertical coordinate, Mon. Weather Rev., 138, 3683–3689, https://doi.org/10.1175/2010MWR3307.1, 2010. a
Mouallem, J., Harris, L., and Benson, R.: Multiple samelevel and telescoping nesting in GFDL's dynamical core, Geosci. Model Dev., 15, 4355–4371, https://doi.org/10.5194/gmd1543552022, 2022. a, b, c
MPIM: Instructions for obtaining the ICON Code, https://code.mpimet.mpg.de/projects/iconpublic/wiki/Instructions_to_obtain_the_ICON_model_code_with_a_personal_noncommercial_research_license, last access: 18 November 2019. a
Narcowich, F. and Ward, J.: Generalized Hermite interpolation via matrixvalued conditionally positive definite functions, Math. Comput., 63, 661–687, https://doi.org/10.2307/2153288, 1994. a
Prill, F., Reinert, D., Rieger, D., and Zängl, G.: ICON Tutorial: Working with the ICON model, Deutscher Wetterdienst (DWD), https://doi.org/10.5676/dwd_pub/nwv/icon_tutorial2020, 2020. a
Reinert, D., Prill, F., Frank, H., Denhard, M., Baldauf, M., Schraff, C., Gebhardt, C., Marsigli, C., and Zängl, G.: DWD Database Reference for the Global and Regional ICON and ICONEPS Forecasting System, Deutscher Wetterdienst (DWD), https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/icon/icon_dbbeschr_aktuell.html (last access: 15 August 2022), 2021. a
Rípodas, P., Gassmann, A., Förstner, J., Majewski, D., Giorgetta, M., Korn, P., Kornblueh, L., Wan, H., Zängl, G., Bonaventura, L., and Heinze, T.: Icosahedral Shallow Water Model (ICOSWM): results of shallow water test cases and sensitivity to model parameters, Geosci. Model Dev., 2, 231–251, https://doi.org/10.5194/gmd22312009, 2009. a
Schär, C., Leuenberger, D., Fuhrer, O., Lüthi, D., and Girard, C.: A new terrainfollowing vertical coordinate formulation for atmospheric prediction models, Mon. Weather Rev., 130, 2459–2480, https://doi.org/10.1175/15200493(2002)130<2459:ANTFVC>2.0.CO;2, 2002. a, b, c, d, e, f
Schmidt, F.: Variable fine mesh in spectral global models, Beitr. Phys. Atmos., 50, 211–217, 1977. a
Skamarock, W. C., Klemp, J. B., Duda, M. G., Fowler, L. D., Park, S.H., and Ringler, T. D.: A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tesselations and CGrid Staggering, Mon. Weather Rev., 140, 3090–3105, https://doi.org/10.1175/MWRD1100215.1, 2012. a, b
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Zhiquan, L., Berner, J., Wang, W., Powers, J. G., Duda, M. G., , Barker, D. M., and Huang, X.: A Description of the Advanced Research WRF Model Version 4, No. ncar/tn556+str, National Center For Atmospheric Research, Boulder, CO, https://doi.org/10.5065/1dfh6p97, 2019. a, b, c
Staniforth, A. N. and Mitchell, H. L.: A VariableResolution FiniteElement Technique for Regional Forecasting with the Primitive Equations, Mon. Weather Rev., 106, 439–447, https://doi.org/10.1175/15200493(1978)106<0439:AVRFET>2.0.CO;2, 1978. a
Tomita, H.: A Stretched Icosahedral Grid by a New Grid Transformation, J. Meteorol. Soc. Jpn. Ser. II, 86A, 107–119, https://doi.org/10.2151/jmsj.86A.107, 2008. a
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, c
Warner, T. T., Peterson, R. A., and Treadon, R. E.: A Tutorial on Lateral Boundary Conditions as a Basic and Potentially Serious Limitation to Regional Numerical Weather Prediction, B. Am. Meteorol. Soc., 78, 2599–2618, https://doi.org/10.1175/15200477(1997)078<2599:ATOLBC>2.0.CO;2, 1997. a
Weimer, M., Buchmüller, J., Hoffmann, L., Kirner, O., Luo, B., Ruhnke, R., Steiner, M., Tritscher, I., and Braesicke, P.: Mountainwaveinduced polar stratospheric clouds and their representation in the global chemistry model ICONART, Atmos. Chem. Phys., 21, 9515–9543, https://doi.org/10.5194/acp2195152021, 2021. a, b
WMO: Manual on the Global Dataprocessing and Forecasting System: Annex IV to the WMO Technical Regulations, WMOno. 485, World Meteorological Organization, Geneva, Switzerland, ISBN 9789263104854, 2019. 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, j, k, l, m, n, o
Zängl, G., Reinert, D., and Prill, F.: Grid Refinement in ICON v2.6.4 (research data), Edmond [data set], https://doi.org/10.17617/3.NOC2AE, 2022. a
Zarzycki, C. M., Jablonowski, C., and Taylor, M. A.: Using VariableResolution Meshes to Model Tropical Cyclones in the Community Atmosphere Model, Mon. Weather Rev., 142, 1221–1239, https://doi.org/10.1175/MWRD1300179.1, 2014. a