Grid Reﬁnement in ICON v2.6.4

. This article describes the implementation of grid reﬁnement in the ICOsahedral Nonhydrostatic (ICON) modeling system. It basically follows the classical two-way-nesting approach known from widely used mesoscale models like MM5 or WRF, but differs in the way how feedback from ﬁne 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 well-established implementations on quadrilateral grids, new methods had to be developed for interpolating 5 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 gradient-based extrapolation from parent to child grid points is used in other cases. The technical implementation on the unstructured ICON grid is based upon sorting the boundary interpolation points at the beginning of the grid-point index vector, so that computations on boundary interpolation points can be excluded by appropriate start indices without the need of IF masks. The run-time 10 ﬂow control is written such that limited-area 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 Schär mountain wave test case (Schär et al., 2002) and the Jablonowski-Williamson test case (Jablonowski and Williamson, 2006) 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 conﬁrmed by experiments closely following the conﬁguration used for operational 15 numerical weather prediction at DWD, which demonstrate that a regional reﬁnement over Europe has a signiﬁcant positive impact on the forecast quality in the northern hemisphere.


Introduction
The ICON (ICOsahedral Nonhydrostatic) modeling framework is jointly developed by the German Weather Service (DWD), the Max Planck Institute for Meteorology (MPI-M), 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 shallow-water framework (Bonaventura and Ringler, 2005;Ripodas 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 https://doi.org/10.5194/gmd-2022-120Preprint.Discussion started: 23 May 2022 c Author(s) 2022.CC BY 4.0 License.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), C-grid type 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, together with the fact that a triangular grid proved to be more suitable for implementing a grid nesting capability than a hexagonal one, led us to retain the triangular C-grid discretization.
Despite impressive advances in computational power over the last decades, the application of global models with uniform, convection-permitting resolution on weather or even climate timescales is still too costly to be performed on a regular basis.To date, high-resolution limited area models (LAM) serve as a cost-effective alternative for exploring the meso-and microscale, and will continue to serve as a working horse for both the NWP and climate community.Limited area models have proven successful, but they are known to have conceptional deficiencies, such as the potential ill-posedness 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-global scale 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 so-called 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 on GCM applications is provided, e.g. by Fox-Rabinovitz 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 region(s) of interest (h-refinement).This approach is pursued, for example, by the Model for Prediction Across Scales (MPAS) (Skamarock et al., 2012), or the spectral element dynamical core of the Community Atmosphere Model (CAM) (Zarzycki et al., 2014).Albeit being much more flexible, this approach has to cope with similar challenges as the grid stretching approach mentioned previously, regarding numerical efficiency and the suitability of the physical parameterization suite.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 horizontally implicit numerical methods are chosen.Depending on the zoom factor, 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 scale aware 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 two-way nesting, as known from many regional mesoscale models such as MM5 (Grell et al., 1994) or WRF (Skamarock et al., 2019).Two-way nesting differs from the previous singlegrid approaches by the fact that multiple grids of different resolution are overlaid into 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 nested grid 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 resolution jump.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 two-way nesting has been successfully applied in limited-area modeling since 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 two-way nesting approach implemented in the High Resolution Atmosphere Model (HiRAM) on a cubed sphere grid (Harris andLin, 2013, 2014).
The purpose of this article is to describe the implementation of grid nesting in ICON.Compared to previous two-way-nesting approaches, it is the first implementation on triangular grids and it differs in the way how 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 limited-area mode is available in ICON as a by-product of the gridnesting implementation, differing from nesting only in the way how 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.

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 https://doi.org/10.5194/gmd-2022-120Preprint.Discussion started: 23 May 2022 c Author(s) 2022.CC BY 4.0 License.icosahedral-triangular Arakawa-C grid in the horizontal, and a terrain following height-based SLEVE coordinate (Leuenberger et al., 2010) with Lorenz-type staggering in the vertical.As described in Zängl et al. (2015), the prognostic variables encompass the edge-normal 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 two-time-level 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 hand and additional sub-grid 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 sub-stepped with respect to physical processes, with a single physics time step consisting of 5 dynamics substeps by default (nsubs = ∆t/∆τ = 5).
The physics-dynamics coupling scheme in ICON further distinguishes between fast physical processes, having a time scale shorter than or comparable to the time step ∆t, and slow processes, which are considered to have a time scale large compared to ∆t.Processes that fall into the category fast are saturation adjustment, grid scale 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 process-specific larger time steps that are an integer multiple of ∆t.We therefore call the time step ∆t the fast physics time step.The fact that tracer transport is performed at ∆t as well, in combination with the wish for consistency with continuity (Gross et al., 2002), necessitates that the coupling of nested domains is performed at the fast physics timestep ∆t rather than ∆τ , as will be described in the following sections.
The static mesh refinement in ICON is accomplished using a multi-grid approach, where one or more higher resolution (child) domains are overlaid on a coarser base (parent) domain.The base domain can be a regional or a 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 centered.
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 overlap.They may also be switched on or off during runtime.Conceptually, the number of nested domains is arbitrary, but of course not all choices would make sense from a physical point of view.Each domain can be regarded as a separate instance 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 convection can be reduced by stronger entrainment or even be 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 4 child triangles (see Fig. 1).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 one-way or two-way.A mixture of one-way and two-way coupled domains is also possible.Two-way versus one-way 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 case of one-way nesting, the feedback is switched off.
To perform the coupling, we conceptionally 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.
Limited-area domains are treated technically like one-way 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 has a fixed width of 4 cell rows.This is motivated by the technical constraint that the boundary zone needs to match with parent cell rows, combined with the fact that 2 cell rows would not be sufficient to cover all stencil operations performed in the dynamical core (e.g.∇ 4 -diffusion; Zängl et al., 2015).It is non-prognostic, 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.The nudging zone, which is active in the case of one-way 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 two-way nesting, no nudging is applied and the boundary interpolation zone borders on 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 relaxation-type 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, ρ, θ v as well as to the mass fractions of water vapour q v , cloud water q c , and cloud ice q i .Child-to-parent 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 cubed-sphere 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 less numerical disturbances near vertical nest interfaces, leading to slightly better forecast quality in NWP applications.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.

Lateral Boundary Update: Parent → Child
The boundary update mechanism provides the child domain with up-to-date lateral boundary conditions for the prognostic variables v n , w, ρ, θ v , q k .In order to prevent parent-to-child 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 M n p and M n c , respectively, where n specifies the time step index.In general, the boundary update works as follows: Let ψ n p , ψ n+1 Concerning the interpolation operator I p→c , we distinguish between cell-based variables (i.e.scalars) and edge-based variables (v n and ρv n ).For cell-based variables a 2D horizontal gradient is reconstructed at the parent cell center by first computing edge-normal gradients at edge midpoints, followed by a 9-point 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 cell based tendencies.
To prevent excessive over-and undershoots of ψ cj in the vicinity of strong gradients, a limiter for ∇ψ p is implemented.It ensures that on all four child points, where ψ p,min and ψ p,max 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 I 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 above-mentioned nudging zone) is applied.
Regarding the interpolation of edge-based variables (i.e. the edge-normal 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).
Edge-normal vector components at the inner child edges are reconstructed using a direct RBF reconstruction based upon the five-point stencil indicated by solid red dots in Fig. 2a.For a given inner child edge the stencil comprises the edges of the 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 pentagon points) edge points adjacent to a vertex.
The edge-normal 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 edge-normal 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 edge-normal direction and taking the centered difference.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 small-scale noise in v n .To suppress the remaining noise arising from the interpolation of the time tendency, a second-order 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 the five grid rows adjacent to the interpolation zone.For the other prognostic variables, no special filtering is applied near nest boundaries.In the case of one-way nesting, the second-order velocity diffusion is extended into the nudging zone of the nested domain, replacing the enhanced fourth-order 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 sub-steps is interpolated instead of time level n, as the time averaged 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 flux-form 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 non-conservative 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.

Feedback: Child → Parent
If two-way nesting is activated, the model state M n+1 p on the parent domain is relaxed towards the updated model state M n+1 c on the child domain at every parent fast physics time step ∆t p .This relaxation-type feedback is only applied to the prognostic variables v n , w, θ v , ρ 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 time scales (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 land-surface 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: 1. Upscaling: The updated variable ψ n+1 c is interpolated (upscaled) from the child domain to the parent domain.The upscaling operators for cell based and edge based variables will be denoted by I c→p and I e c→p , respectively.

Difference computation:
The difference between the parent-domain variable ψ n+1 p and the upscaled child variable ) 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 edge based normal velocity v n , a simple arithmetic average of the two child edges lying on the parent edge is taken.
For cell based 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 ( 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 = 1, . . ., 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 {c 2 , c 3 , c 4 } (see Fig. 2a), where the weights {α 2 , α 3 , α 4 } represent the barycentric coordinates.Generally, the constraints ( 6) and ( 5) ensure reversibility in the sense that I c→p (I 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 equations ( 5), ( 6) is closed with a final constraint which reads as where a c1 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 pre-defined 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 {x 2 , x 3 , x 4 }, 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 cell based operator I c→p is not strictly mass-conserving and that strict mass conservation would require some means of area-weighted aggregation from the child cells to the parent cells, which is available as an option.The problem 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.Using an area-weighted aggregation from the child cells to the parent cells thus would map linear horizontal gradients on the child grid into a checkerboard noise pattern between upward and downward oriented 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 area-weighted aggregation of the orographic raw data to the model grid.Feeding back ρ is thus intrinsically non-conservative.To keep the related errors small and non-systematic, 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 area-weighted averaging are unimportant (with real orography) compared to the mesh-size-related conservation error.
When combining the above mentioned steps, the feedback mechanism for ρ can be cast into the following form: Here ρ n+1 p 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 τ f b is a user-defined relaxation time scale that has a default value of τ f b = 10800 s and is independent of the relaxed field.The smaller τ f b , the more the parent state is drawn towards the child state.The chosen default value is optimized for our typical NWP applications and aims at filtering small-scale transient features from the feedback, while fully capturing synoptic-scale 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 )) 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 re-diagnosed 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-to the parent domain and added to the parent reference profile θ v ref,p .
The same approach is taken for w, however the full field is upscaled.
In the case of v n some second-order diffusion is added to the ensuing feedback increment in order to damp small-scale noise.
with the feedback increment and the diffusion coefficient ∆tp , where a p,e is the area of the quadrilateral spanned by the vertices and centers adjacent to the parent's edge.

Lateral Nudging
If the feedback is turned off, i.e. if one-way 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 (Eq.( 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 (Eq. ( 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 e-folding width µ, the latter two 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 second-order diffusion on v n is used near the lateral nest boundaries in order to suppress small-scale noise.As opposed to Eq. ( 10) for feedback, the diffusion operator acts on the velocity field itself, rather than the velocity increment.

Vertical Nesting
The vertical nesting option allows to set model top heights individually for each domain, with the constraints that the child domain height is lower or at most equal to the parent domain height, and that the child domain extends into heights where 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 a significant amount of 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 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.Boundary conditions for v n , w, θ v , ρ, q k as well as the vertical mass flux ρw are derived as follows: For w, θ v , ρ and ρw the full fields at the nest interface level are interpolated from the parent to the child grid, using the same RBF based 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 ψ ∈ {w, θ v , ρ, ρw} the parent to child interpolated field reads with s denoting an individual dynamics substep, and nsubs denoting the total number of substeps.We apply the same interpolation method to the corresponding time tendencies ∂ψ p /∂t, which are estimated by taking the difference of the state variables 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 half level below (denoted as ∆v n,p in the following) are interpolated rather than the full field, using again 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 = 3/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 the 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 tracerand air mass consistency.Instead, we make use of the already interpolated mass flux (ρw) 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 (ρw) c computed via Eq.(11).

Recursive Algorithm for Multi-domain Setups
So far, we have focused on the coupling of an individual parent and child domain.The coupling of multiple and possibly repeatedly nested domains requires a well conceived processing sequence, whose basics will be described in the following.rameterizations for each domain in basically the same way as for the global domain.The basic processing sequence is as follows: 1.A single integration step with ∆t is executed on the global domain, resulting in an updated model state M n+1 p , as indicated by an open black circle in Fig. 3.
2. Boundary data are interpolated from the global domain to nest 1 (red arrow), followed by an integration step on nest 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) (JW) test 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 equator-to-pole 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 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 time scale τ f b 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.
As discussed in Zängl et al. (2015), the general behaviour 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 wave-like 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  , and E5 (i,j) after 9 days of integration.Dashed boxes show the nest locations.See Table 1 for the experiment configurations.
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 (not shown), the position and intensity of the cyclones is now very similar to E5.The relaxation-based 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 445 above-mentioned long relaxation time scale of 12 hours, combined with the fact that the middle and left cyclones propagate at a somewhat slower speed in the coarse domain.Reducing τ f b gradually decreases the differences, but for values less than about 5∆t p (i.e. two hours), the solution in the nested domain starts to deviate from the E5 reference because numerical disturbances generated along the nest boundaries become non-negligible (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 that the intensity bias encountered here is not of practical relevance.On the other hand, the standard relaxation time scale of three hours filters transient wave-like motions that are not deterministically predictable anyway, thus avoiding a transfer of small-scale 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 derivative-based quantities like vorticity can be calculated depends directly on the resolution of the model grid.To obtain better comparability, Fig. 5 displays vorticity fields for the leading cyclone computed at the same grid resolution (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 to 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.
To further examine the flow disturbances generated by the resolution jump at the nest boundary, Fig. 6 compares the surface pressure fields for the steady-state JW test after 10 days of integration for E1 and E2.This variant of the JW test accentuates the so-called grid imprinting errors arising from any irregularity in the model grid, in our case both due to the non-uniformity of the icosahedral grid and due to the resolution jump along the nest boundaries.The result for E1 (Fig. 6a) shows the wellknown regular wavenumber-five disturbance pattern characteristic for icosahedral grids.With the presence of a nested domain, the disturbances become irregular and reach a somewhat stronger peak amplitude (Fig. 6b).The largest disturbances occur in the airmass that was initially located in the nest region because the nest-induced perturbations there have the longest time to grow.For smaller values of τ f b , the amplitude of these disturbances gradually increases (not shown).However, at finer model resolution, the nest-induced perturbations decrease about proportionally to the grid imprinting in the global grid (see Fig. 2 in Zängl et al., 2015), and they have never become noticeable in our operational NWP applications.

Schär Mountain Test with Vertical Nesting
In order to examine the behaviour of the upper boundary condition for vertically nested domains, we conducted the Schär 475 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 small-scale, vertically decaying non-hydrostatic gravity waves and largerscale vertically propagating quasi-hydrostatic 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.

480
Following Zängl et al. (2015), we performed a quasi-three-dimensional analogue of this test, where the original 1D mountain profile is changed to a 2D ridge-like profile that decays towards zero in cross-flow direction, as the lateral domain boundary is approached.It is given by The wind speed is set to U = 10 ms −1 , and the Brunt-Väisälä-frequency is given by N = 0.01 s −1 for z < 20 km, and gradually increases to N = 0.03 s −1 above.
We have performed two types of simulations, one with and one without a vertically nested domain.The reference simulation without vertical nest was performed on a single 3 • × 3 • limited area 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 steady state solution for w is depicted in Fig. 7a.It may be compared against earlier ICON results in Zängl et al. (2015), or results from other numerical models at comparable resolutions (e.g.Schär et al., 2002;Skamarock et al., 2012).
For the nested simulation, a 2.4 • × 2.4 • R2B12 domain with 39 vertical layers and a vertical interface at 20 km is two-way 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.Spurious disturbances are largely confined to the uppermost quasi-hydrostatic wave crest and trough, and to the leeward propagating wave signal.This does not only hold for the steady state solution, but also for the spinup phase (not shown).The reason for the remaining disturbances is twofold: Firstly, the computation of the boundary conditions for 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 grid cannot exactly match the solution on the child grid.
This difference depends on the feedback time scale, and it can be seen by comparing Fig. 7b (τ f b = 10800 s) with Fig. 7c (τ f b = 900 s) that a shorter time scale further reduces the difference to the non-nested reference.The solution on the R2B11 parent domain (with τ f b = 900 s) is depicted in Fig. 8. Due to its two times coarser mesh size it is lacking some of the leeward propagating wave signals and shows a slightly reduced amplitude of the quasi-hydrostatic wave.We note in addition that the attenuated, though visible, leeward propagating wave signal in Fig. 8 results from the child to parent feedback and does not exist in a R2B11 simulation without 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 resolution-induced differences between the parent and child model states.On the other hand, this test has shown the ability of the child to parent feedback mechanism to have a small but noticeable positive impact on the quality of the nest interface conditions.

Operational NWP Applications
In an operational context, one ideally expects that the beneficial impact on forecast quality of the regionally refined resolution is transferred to the global domain in the nest overlap area and subsequently propagates downstream, which is usually eastward 520 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 two-way nested domain referred to as "ICON-EU" with 6.5 km mesh size covering Europe and some adjacent regions (24.5 • W -63.5 • E; 29 • N -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 manuscript) uses mesh sizes of 40 km and 20 km, respectively.
To demonstrate the benefit of two-way 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 Medium-Range 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 nest (R2B6), and a global 20-km mesh without nest (R2B7).All experiment suites are conducted for January 2021, starting from the 00UTC 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 above-mentioned verification against IFS analyses, we use our 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. 9 shows the 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. 10, showing relative differences to R2B6 for easier readability and including an additional test in which the nest stays active until 180 h (R2B6N7-180h).It is seen that during the first three 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 five days, the errors in R2B6N7 become temporarily a bit larger than in R2B6 over Europe, and we infer from the negligible differences between R2B6N7 and R2B6N7-180h 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 (Figs. 11 and 12), which include an assessment of statistical significance at the 95% level.The SYNOP verification shows a pronounced beneficial impact in the nest overlap region (EU-NEST) during the first five 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 throughout the troposphere in both Europe and Asia.In the lower stratosphere, a slight degradation can be seen over Europe for wind speed / direction and 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 double-penalty 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 EU-nest 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 upper-air 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 that that of the surface stations, giving Europe an overproportional 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 two-way nesting capability in ICON fully serves its purpose.

Conclusions
This article provided a technical description of the grid nesting implementation in the ICOsahedral Nonhydrostatic (ICON) modelling system.The available options comprise one-way and two-way 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 limited-area mode is available as a by-product of the grid-nesting implementation, which differs from one-way nesting only in the way the lateral boundary conditions are provided.The model-internal 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 relaxation-type 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 well-established nesting implementations on quadrilateral grids was the development of appropriate interpolation algorithms for the lateral boundary conditions of the nested grids, as the usual higher-order polynomial interpolation leads to the inversion of a nearly singular matrix on the triangular ICON grid.It was found that radial basis functions (RBFs) 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 well as 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 condi- tion, 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 non-nested 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 steady-state variant of the test nevertheless shows that this kind of disturbances 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 steady-state multi-scale orographic gravity wave composed of vertically propagating and vertically decaying wave components.The results show that the wave reflections generated at the vertical nest 605 interface are small enough to be irrelevant for practical applications.The findings from these idealized tests are corroborated by full-physics 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 where the forecast quality has a pronounced In regions with an overlapping child domain the counters for cells, edges and vertices are filled with negative numbers (see light-blue 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 3 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 with a child domain allows to easily access all the points involved into the child to parent feedback process described in Sect.2.2.2.

A2 Distributed-Memory 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 of masking operations within DO-loops, the grid points are ordered by their distance from the lateral boundary by making use of the refin_ctrl 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 non-existing 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 come behind in a non-ordered fashion, the refin_ctrl indexing proceeding with the remaining gridpoints in the nudging zone, the prognostic grid points not overlapping with a child domain (refin_ctrl = 0), and the overlapping prognostic points (refin_ctrl < 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.
Regarding distributed-memory (MPI) parallelization, the general 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 resolution 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 done during the grid generation process by indicating a list of domains.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 one-way and two-way nesting cannot be mixed within one logical domain, there may still be two logical domains on a given nest level.
To further optimize the amount of MPI communication, a so-called 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.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.

Figure 1 .
Figure 1.Basic structure of a nested or limited-area domain.Orange: boundary interpolation zone, having a fixed width of 4 cell rows.Ocher: nudging zone with adjustable width, only active for one-way nesting and in limited area mode.Blue and light-blue: Child to parent feedback zone.Light-blue: Another child domain overlapping with the depicted domain.Prognostic computations are restricted to the feedback and nudging zone.Black and white integers indicate the internal indexing which is used to assign cells and edges to individual zones.More details on the indicated sorting of cell rows are given in Appendix A2.
2.2.1 for https://doi.org/10.5194/gmd-2022-120Preprint.Discussion started: 23 May 2022 c Author(s) 2022.CC BY 4.0 License.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 fast physics timestep ∆t c .

p
represent a prognostic variable on the parent domain at time steps n and n + 1, respectively.Once the variables on the parent domain M p have been updated from n to n + 1, diagnosed.Both the field ψ n p at time level n and the tendency ∂ψp ∂t are then interpolated (downscaled) from the parent grid cells/edges to the corresponding cells/edges of the child domain's boundary zone (orange cells in Fig. 1).With I p→c denoting https://doi.org/10.5194/gmd-2022-120Preprint.Discussion started: 23 May 2022 c Author(s) 2022.CC BY 4.0 License. the interpolation operator, we get ψ n c = I p→c ψ n p ∂ψ c ∂t = I p→c ∂ψ p ∂t .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 M n+1 c , with each step consisting of several (typically 5) dynamics sub-steps.The temporal update is performed at each dynamics sub-step ∆τ 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 sub-step of the first and second (fast physics) integration step on the child domain read ψ n c and ψ n c + 0.5 ∆t p ∂ψ c /∂t, respectively.

Figure 2 .
Figure 2. Horizontal reconstruction stencil for edge-normal vector components at (a) inner child edges and (b) outer child edges.The child edge under consideration is highlighted in red.Black open dots indicate child edge midpoints, while black solid dots indicate cell circumcenters.Solid red dots represent the reconstruction stencil, i.e. the location of the parent edge-normal vector components entering the reconstruction, and blue triangles in (b) indicate the location of the reconstructed 2D vectors.See the text for details.
Since by construction d(p, c 1 ) = −d(p, c 2 ) holds on the ICON grid, the above mentioned mass flux consistency is ensured.It is noted that attempts to use higher-order 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.https://doi.org/10.5194/gmd-2022-120Preprint.Discussion started: 23 May 2022 c Author(s) 2022.CC BY 4.0 License.
p , with the parent-child difference in the reference density field ∆ρ ref,p = I c→p (ρ ref,c ) − ρ ref,p , and the potential temperature perturbation θ n+1 v,c = θ n+1 v,c − θ v ref,c .The term ∆ρ ref,p is a pure function of the parent-child height difference and can be viewed as a leading-order correction term.As a further optimization, the empirical factor (1.05 − 0.005 I c→p (θ n+1 v https://doi.org/10.5194/gmd-2022-120Preprint.Discussion started: 23 May 2022 c Author(s) 2022.CC BY 4.0 License.
https://doi.org/10.5194/gmd-2022-120Preprint.Discussion started: 23 May 2022 c Author(s) 2022.CC BY 4.0 License.at the substeps s = 1 and s = nsubs.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.

Figure 3
Figure 3 displays a schematic example where a global domain is combined with two successively (two-way) 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 top down.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 limited-area domain) down to the deepest nesting level and calls the dynamical core and the physics pa-

Figure 3 .
Figure 3. Schematic of a global domain with two successively nested domains (two-way).The processing sequence for the time integration of all domains from time step n to n + 1 is shown in the flowchart at the lower left.See the text for details.

Figure 6 .
Figure 6.Surface pressure (hPa) after 10 days of integration for the steady-state JW test and experiments E1 (a) and E2 (b).Dashed box shows the nest location.
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, as well as the horizontal resolution of the nested domain, exactly match those of the reference simulation.The only significant difference lies in the much lower, but open domain top.Steady state solutions for the nested domain are depicted in Fig. 7b,c.They differ in terms of the relaxation time scale τ f b for child to parent feedback, with values of τ f b = 10800 s and τ f b = 900 s for panels (b) and (c), respectively.

Figure 7 .
Figure 7. Vertical velocity w after 12 h for the Schär et al. (2002) test case (contour interval 0.05 ms −1 ).(a) Reference results for a single, non-nested R2B12 domain with 40 km model top.(b),(c) Results for a vertically nested R2B12 domain with 20 km model top, and two distinct relaxation time scales for child to parent feedback.The entire vertical extend of the nested domain is depicted.

Figure 8 .
Figure 8. Vertical velocity w after 12 h for the R2B11 parent domain, with τ f b = 900 s.

Figure 9 .
Figure 9. Verification against IFS analyses for the northern hemisphere (latitude > 20 • N) at 500 hPa: RMS errors of geopotential (m) and vector wind (ms −1 ) as a function of lead time (h).See text for experiment acronyms.

Figure 10 .
Figure 10.Verification against IFS analyses for Europe (35-72 • N, 12 • W-45 • E) at 500 hPa: Relative RMS difference (%) of geopotential and vector wind w.r.t.R2B6.Note that R2B6N7 and R2B6N7-180h are identical up to a lead time of 120 h.See text for experiment acronyms.
560shown).The results for Asia show that the improvements related to the EU-nest propagate downstream with a delay of about three days.They are a bit smaller than in the nest region but still partly reach the 95% significance level, indicating an obvious benefit of our 2-way nesting methodology.The TEMP verification (Fig.12) additionally shows that the positive signal extends https://doi.org/10.5194/gmd-2022-120Preprint.Discussion started: 23 May 2022 c Author(s) 2022.CC BY 4.0 License.

Figure 12 .
Figure 12.Scorecard for verification against radiosonde ascents for wind direction (DD), wind speed (FF), relative humidity (RH), temperature (T) and geopotential (Z) in the same regions as in Fig. 11.Colours indicate relative RMSE differences between R2B6N7 and R2B6, and wide boxes indicates statistical significance at the 95%-level.
the time interval ∆t/2, and resulting in the model state M 6.As a final step, feedback is performed from nest 2 to nest 1, followed by feedback from nest 1 to the global domain.https://doi.org/10.5194/gmd-2022-120Preprint.Discussion started: 23 May 2022 c Author(s) 2022.CC BY 4.0 License.

Table 1 .
List of model configurations for the Jablonowski-Williamson test.