Interactive comment on “ The ICON-1 . 2 hydrostatic atmospheric dynamical core on triangular grids – Part 1 : Formulation and performance of the baseline version ” by H

Abstract. As part of a broader effort to develop next-generation models for numerical weather prediction and climate applications, a hydrostatic atmospheric dynamical core is developed as an intermediate step to evaluate a finite-difference discretization of the primitive equations on spherical icosahedral grids. Based on the need for mass-conserving discretizations for multi-resolution modelling as well as scalability and efficiency on massively parallel computing architectures, the dynamical core is built on triangular C-grids using relatively small discretization stencils. This paper presents the formulation and performance of the baseline version of the new dynamical core, focusing on properties of the numerical solutions in the setting of globally uniform resolution. Theoretical analysis reveals that the discrete divergence operator defined on a single triangular cell using the Gauss theorem is only first-order accurate, and introduces grid-scale noise to the discrete model. The noise can be suppressed by fourth-order hyper-diffusion of the horizontal wind field using a time-step and grid-size-dependent diffusion coefficient, at the expense of stronger damping than in the reference spectral model. A series of idealized tests of different complexity are performed. In the deterministic baroclinic wave test, solutions from the new dynamical core show the expected sensitivity to horizontal resolution, and converge to the reference solution at R2B6 (35 km grid spacing). In a dry climate test, the dynamical core correctly reproduces key features of the meridional heat and momentum transport by baroclinic eddies. In the aqua-planet simulations at 140 km resolution, the new model is able to reproduce the same equatorial wave propagation characteristics as in the reference spectral model, including the sensitivity of such characteristics to the meridional sea surface temperature profile. These results suggest that the triangular-C discretization provides a reasonable basis for further development. The main issues that need to be addressed are the grid-scale noise from the divergence operator which requires strong damping, and a phase error of the baroclinic wave at medium and low resolutions.


Introduction
In the development of general circulation models (GCMs) for the purpose of numerical weather prediction (NWP) and climate research, one of the central tasks is to design numerically accurate, robust, and computationally efficient algorithms to solve the resolved-scale fluid dynamics equations that govern the atmospheric motions.This component of a GCM, commonly referred to as the dynamical core, provides meteorological background for and interacts with tracer transport as well as diabatic processes like radiative heating/cooling and cumulus convection.Numerical properties of the dynamical core play an important role in determining the behavior of an atmospheric model.
Conventionally, operational NWP and climate models are developed on latitude-longitude grids.The convergence of meridians and the resulting clustering of grid points near the North and South Poles not only causes numerical issues (Williamson, 2007), but also is expected to become a bottleneck in massively parallel computing (Staniforth and Thuburn, 2012).In the past two to three decades, the idea of using geodesic grids, first proposed in the late 1960s by Williamson (1968) and Sadourny et al. (1968), has been revisited by many researchers.The efforts of Heikes and Randall (1995a,b), Ringler and Randall (2002), Majewski et al. (2002), Tomita and Satoh (2004) and Giraldo and Rosmond (2004) are examples that developed models on icosahedral grids.Other grid configurations such as the cubed sphere (e.g., McGregor, 1996;Giraldo et al., 2003;Putman and Lin, 2007), Yin-Yang grid (e.g., Kageyama and Sato, 2004;Qaddouri et al., 2008;Li et al., 2007), and Fibonacci grid (Swinbank and Purser, 2006), have also been investigated.The numerical discretizations applied in these models include not only finite-difference and finite-volume schemes (e.g., Heikes and Randall, 1995a;Putman and Lin, 2007), but also modern techniques such as high-order continuous and discontinuous Galerkin methods (e.g.Taylor et al., 1997;Giraldo et al., 2003;Läuter et al., 2008).General overviews of the historical development and recent advances in this research area can be found in Staniforth and Thuburn (2012) and Bonaventura et al. (2012), along with more complete lists of references.
In the year 2001, a collaboration was initiated between the Max Planck Institute for Meteorology (MPI-M) and the German Weather Service (DWD) to develop a unified model system for weather prediction and climate applications.Based on previously identified model limitations in the basic assumptions, the conservation properties, as well as the computational performance, the following set of goals were set for the new development: 1.A nonhydrostatic dynamical core for high-resolution global weather forecasting and cloud-resolving-scale climate modeling.
2. Conservation of air and tracer mass, as well as consistency between the discrete continuity equation and the tracer transport algorithm.
3. The capability of two-way nesting with multiple refined regions per nesting level, to replace the current operational set-up at DWD that contains two different models (a global hydrostatic model and a regional nonhydrostatic model).
4. Capabilities of one-way horizontal nesting, vertical nesting, and stand-alone limited-area mode.
5. Scalability and efficiency on parallel computing architectures with O(10 4 ) or more cores.
Following the GME model (Majewski et al., 2002) of DWD which was already using geodesic grids to deliver operational forecasts at the time, the icosahedral grids were chosen for their quasi-uniformity, giving the name "ICON" (ICOsahedral Nonhydrostatic models) to the joint project.
To achieve the desired conservation properties (goals 2 and 3), and meanwhile keep the discretization stencils relatively small to reduce the amount of communication in parallel computing (goal 5) it was decided to employ the finitevolume/finite-difference discretization approach with C-type staggering.The remaining question then was what control volume to use.
A spherical icosahedral grid obtained by projecting and recursively subdividing an icosahedron (Sect.3) can be viewed as consisting of either triangular or hexagonal/pentagonal cells (cf.Fig. 12 in Staniforth and Thuburn, 2012).While a triangular cell can be easily divided into sub-triangles, a hexagonal cell can not exactly overlap a set of hexagons at higher resolutions.The hierarchical structure of the triangular mesh is thus very attractive for the implementation of mass-conserving discretizations for multi-resolution capabilities such as one-way and two-way nesting.Orographic features could be described with greater precision without directly affecting the global mesh, and physically important boundary fluxes and soil properties could be computed on a finer mesh and then exactly summed up onto a coarser one used by the other model components.Moreover, in NWP and climate applications, it is common for some of the parameterized diabatic processes to employ temporal and/or spatial resolutions lower than those used by the dynamical core, in order to reduce the computational cost of the whole model.Along this line, a triangular mesh allows for calculating, for example, radiative transfer with reduced spatial resolution in a straightforward way.For a hexagonal grid, it is far less obvious how to use a similar concept for particularly expensive physics parameterizations.These considerations have led to the preference of triangles as control volumes for the ICON models.
Regarding the numerical properties of the discretization methods, Ničković et al. (2002) showed that a C-grid staggering on the hexagonal grid produced an undesirable geostrophic mode with non-zero phase speed.This was considered, at the time, as a severe drawback and the major obstacle to applying such a grid to fluid dynamics on a rotating sphere.In the early phase of the ICON development, both triangular and hexagonal C-discretizations were implemented and tested for the shallow water equations, during which a computational mode similar to that described by Ničković et al. (2002) was observed.The absence of such a mode on the triangular grids, together with the abovementioned considerations on grid refinement, led to the decision to use the triangular version of the icosahedral grids.A complete horizontal discretization concept was presented by Bonaventura and Ringler (2005) in the context of solving the shallow water equations on the sphere.Further testing and a comparison with the GME model was carried out by Rípodas et al. (2009), in which it was concluded that the discretizations implemented in the ICON shallow water model were more accurate than those in GME, and meanwhile had better conservation properties than GME.
When the ICON shallow water model was further developed into the hydrostatic dynamical core described in this paper, grid-scale noise was noticed in numerical experiments, and later recognized as part of the inherent feature of the triangular C-grid (Sect.4, see also Wan, 2009).Around the same time, Thuburn (2008) found a way to overcome the geostrophic mode problem on regular planar hexagonal grid, who later extended the solution to arbitrarily structured Cgrids (Thuburn et al., 2009).Their discretizations then became the basis of the MPAS model (Model for Prediction Across Scales, Skamarock et al., 2012).
As is shown in later sections of this paper, our investigations revealed that on the one hand it was possible to suppress the grid-scale noise on the triangular C-grid through carefully chosen numerical diffusion, thus obtaining reasonable results in dry dynamical core tests.On the other hand, concerns remained that when parameterized diabatic processes were introduced into the model, their interactions with the dynamical core might amplify the grid-scale noise, and the fine-tuned numerical diffusion might have harmful effects in long-term climate simulations.Because of these concerns, a hexagonal model was developed by Gassmann (2013).Theoretical analysis revealed some clear advantages over triangles (see, e.g., Gassmann, 2011), although results from idealized dry dynamical core tests did not show dramatic differences as originally expected.
In this paper we describe the triangular version of the ICON hydrostatic dynamical core (hereafter referred to as the ICOHDC) developed based on the work of Bonaventura and Ringler (2005), analyze its numerical properties, and carry out a series of idealized tests to assess the ability of the dynamical core to reproduce important features of the large scale adiabatic and diabatically forced circulation.The purpose of the paper is not, however, to compare the triangular dynamical core with its hexagonal variant in terms of accuracy and efficiency, or to propose a choice of one discretization approach over another, but rather to provide the basis for a discussion on these issues.
The hydrostatic dynamical core is considered as an intermediate step towards the ultimate goal (i.e. a nonhydrostatic core) of the ICON project.At MPI-M rich experience has been accumulated in the past decades with the ECHAM model series which employ the hydrostatic assumption (Roeckner et al., 1992(Roeckner et al., , 1996(Roeckner et al., , 2006;;Giorgetta et al., 2012).In order to assess potential benefits or drawbacks of the triangular C-grids in an isolated manner, i.e. to separate them from other aspects such as choice of governing equa-tions and dynamics-physics interaction, the hydrostatic dynamical core we describe here has been developed and tested using the same governing equations as in ECHAM.Again to facilitate testing and understanding, the first version of the ICOHDC uses the same vertical discretization and time stepping method as in the spectral transform dynamical core of ECHAM.
The paper is organized as follows: Sects. 2 and 3 present the governing equations and the computational mesh used in the ICOHDC; Sects. 4 and 5 discuss the discrete formulation; Sects.6 and 7 evaluate the performance of the dynamical core in idealized test cases.The conclusions and an outlook are given in Sect.8.

Governing equations
The governing equations solved by the hydrostatic dynamical core are the primitive equations in velocity-temperature form.A generic pressure-based terrain-following vertical coordinate η is used, the value of which ranges from 0 at the top of the model to 1 at the earth's surface.The pressure value at the model top is set to 0 hPa to include the total mass of the atmosphere.The prognostic equations of our model read They are complemented by two diagnostic equations, Here, v stands for the horizontal wind vector, ∇ denotes the horizontal gradient, and k is the local unit vector pointing to the upward direction.ζ = (∇ × v) • k denotes the vertical component of the relative vorticity.K is the kinetic energy per unit mass, K = |v| 2 /2.φ is the geopotential.η = dη/dt stands for the material derivative of η.All the other symbols have their conventional meanings.
The horizontal grid generation algorithm used here starts from projecting a regular icosahedron onto the sphere, with two of the twelve vertices coinciding with the North and South Poles.The other five vertices in each hemisphere are located along the latitude circle of 26.6 • N/S with equal longitude intervals of 72 • .The second step of grid generation is to divide each great circle arc of the projected icosahedron into n r arcs of equal length, and each icosahedron face into n 2 r small triangles, as described by Sadourny et al. (1968).This is referred to as the root division.The resulting mesh is designated as the grid level 0. Further mesh refinement is achieved by bisecting each spherical triangle edge and connecting the midpoints by great circle arcs, yielding four small cells for each parent triangle.In our terminology, a root division of the original spherical icosahedron edge into n r arcs followed by n b recursive edge bisections leads to the "R n r B n b grid".
From grid level 0 onwards, only the twelve icosahedron vertices are surrounded by five triangles (they are thus also referred to as the pentagon points), while the other vertices are surrounded by six cells.This introduces irregularity into the grid, resulting in inequality in cell areas and edge lengths.Triangles closest to the pentagon points feature the most severe deformation.Like in Bonaventura and Ringler (2005, hereafter referred to as BR05), C-staggering is applied to the triangular cells by placing mass and temperature at triangle circumcenters.This particular choice of cell center (as opposed to, e.g.barycenter) results in the property that the arc connecting two neighboring mass points (i.e., the dual edge) is orthogonal to and bisects the shared triangle edge.These bisection points are used as velocity points, at which the component of horizontal wind perpendicular to the edge (denoted by v n in this paper, cf.Fig. 1) is predicted using Eq.(1).The velocity points, however, do not bisect the dual edges due to the irregularity of the spherical grids, causing a component of first-order discretization error in the directional gradient calculated at these locations (cf. the normal gradient operator defined in Sect.4.1 and Fig. 2b).In the shallow water model the grid optimization algorithm of Heikes and Randall (1995b) was employed to reduce the off-centering (Rípodas et al., 2009).The Heikes-Randall optimization produces more severely deformed triangular cells than on the unoptimized grid, which affects the performance of the hydrostatic core at medium and low resolutions.Here we use the the grid optimization method of Tomita et al. (2001) for the ICOHDC which connects the triangle vertices by identical linear springs of a tunable spring coefficient β.The iterative optimization algorithm relocates the vertices until the spring system reaches the lowest potential energy.The resulting icosahedral grids feature small off-centering that decreases with increasing resolution, moderate deformation of the triangles, and smooth transition of geometric properties throughout the horizontal domain.The parameter β is set to 0.9 in this work based on inspection of the grid properties and results from dynamical core tests.Although Table 1 reveals that the area ratio of the largest to smallest cells increases with resolution, the slight loss of uniformity does not show clear impact on the testing results.
In the vertical, the widely used hybrid p-σ coordinate of Simmons and Strüfing (1981) ("coordinate 4" in their report) is employed.The staggering follows Lorenz (1960), meaning that the horizontal wind and temperature are carried at "full levels" representing layer-mean values, while the vertical velocity is diagnosed at "half levels" (i.e.layer interfaces, cf.Fig. 1).The vertical grid is identical to that used in the ECHAM models (e.g., Roeckner et al., 2006).

C-grid discretization
In this section we briefly describe the C-grid discretization inherited from the ICON shallow water model of BR05, then present an analysis of its properties.

Basic operators
BR05 established a spatial discretization method for solving the shallow water equations on the spherical triangular grid described in the previous section.This method forms the basis for the hydrostatic model discussed here.Their discretization concept is a mimetic finite difference scheme consisting of the following elements: -The discrete model predicts the normal component of the horizontal wind v n with respect to triangle edges.The tangential component v τ , needed for the vorticity flux term in the momentum equation, is reconstructed from the normal components using vector radial basis functions (cf.references in Sect.5.2).
-Horizontal derivatives are represented by four discrete operators.The divergence operator div(v) applies the Gauss theorem on each triangular control volume to approximate the spatially averaged divergence over that cell (Fig. 2a).The curl operator curl(v) uses the Stokes' theorem to approximate the vertical component of the relative vorticity averaged over a dual (hexagonal or pentagonal) cell centered at a triangle vertex and bounded by arcs connecting the centers of all triangles sharing the vertex (Fig. 2a).The directional derivative of a scalar field at the midpoint of a triangle edge in the normal direction, grad n (•), is approximated by a straightforward finite-difference discretization involving two cell centers, and referred to as the (normal) gradient operator (Fig. 2b).The horizontal derivative tangential to the edge, grad τ (•), is also defined at the edge midpoint, approximated by a central difference using values at the two ends of the edge (Fig. 2b).The mathematical formulations of the four operators are given by Eqs. ( 4), ( 5), (7), and (8) in BR05.
-Higher order spatial derivatives (Laplacian and hyper-Laplacian) are constructed from the four basic operators outlined above (see, e.g.Eq. 14 in Sect.4.3).These derivatives are needed, for example, in semi-implicit time stepping schemes and for horizontal diffusion.
This discretization scheme is conceptually the same as the widely used C-type discretization on quadrilateral grids.The basic operators are simple, but nevertheless have nice properties.The divergence operator, essentially a finite-volume discretization, makes it straightforward to achieve mass conservation.The curl operator guarantees that the global integral of the relative vorticity vanishes.The divergence and gradient operators are mimetic in the sense that the rule of integration by parts has a counterpart in the discrete model (cf.Eqs. 9 and 10 in BR05), a convenient property for achieving conservation properties.The basic operators are also highly localized (i.e.defined on small stencils), which is beneficial in massively parallel computing.
However, a question remains whether the good properties of the quadrilateral C-grids in terms of the faithful representation of inertia-gravity waves are inherited by the triangular C-grid without limitation.For the hexagonal/pentagonal grids (which can been viewed as the dual meshes of the triangular grids), the wave dispersion analysis in Ničković et al. (2002) revealed that discretization approaches using C-staggering could produce spurious geostrophic modes.A technique to avoid such modes on the hexagonal/pentagonal meshes was proposed in Thuburn (2008) and further developed in (Thuburn et al., 2009).On the triangular C-grids, spurious modes have also been noticed (e.g. Le Roux et al., 2007;Danilov, 2010;Weller et al., 2012).Some recent articles discussed this issue by analyzing the linearized shallow water equations and the representation of vector fields in a trivariate coordinate system (Thuburn, 2008;Danilov, 2010;Gassmann, 2011;Weller et al., 2012).Here, we take a different perspective and use truncation error analysis to show that the divergence operator on the triangular C-grid described above inherently produces gridscale checkerboard error patterns.The same analysis leads to a proposal for estimating the specific amount of numerical hyper-diffusion necessary to reduce the impact of these systematic errors on numerical simulations.

Truncation error analysis
To focus on the triangular geometry and the C-staggering, we carry out the analysis on a planar grid consisting of equilateral triangles of edge length l.We associate a local Cartesian coordinate to each cell with the origin located at the triangle center and the x-axis parallel to one edge.We  then denote the normal outward unit vector at edge midpoints by n j where j ∈ [1, 3] is the edge index.A label δ is assigned to each cell to denote its orientation, with values of 0 and 1 for upward-and downward-pointing triangles, respectively.Thus, the three neighboring cells sharing edges with an upward-pointing triangle are downwardpointing, and vice versa.In the truncation error calculation, a downward-pointing triangle is understood as the image of the corresponding upward-pointing triangle mirrored at the x-axis (Fig. 3).For a generic vector field v differentiable to a sufficiently high order, we denote its components in the x and y directions by u and v, respectively.Applying the discrete divergence operator pointwise values of the vector field known at edge midpoints, denoted as v j , the 2-D Taylor expansion yields Here the subscript o denotes the function evaluation at the triangle center.The functions H and F read Equation ( 7) indicates that the discrete divergence operator applied to pointwise values of v is a first-order approximation of the divergence at the triangle center.More importantly, the first-order error term changes sign from an upward-pointing triangle to a downward-pointing one, which results in a checkerboard error pattern.11), calculated using Eq. ( 6) on a planar triangular grid with 10.4 • resolution in the x-direction.(b) The l 1 , l 2 and l ∞ error norms at different resolutions.The discrete divergence is calculated by first evaluating Eq. ( 11) at edge centers then applying operator (6).Numerical error is computed with respect to cell average.
In a finite-volume perspective, the divergence computed by the Gauss theorem should be interpreted as cell average rather than a pointwise value.However, it is worth noting that, in an equilateral triangle, the cell-center value can be viewed as a second-order approximation of the cell average.Therefore, the first order error term in Eq. ( 7) will also be present in the approximation of the cell average.Indeed, it is not difficult to check that where () stands for cell average.The leading error remains first order and also features a checkerboard pattern.In order to check empirically the impact of Eq. ( 10), numerical calculations have been performed using the vector field the divergence of which reads The discrete divergence is calculated by first evaluating Eq. ( 11) at edge centers and then applying operator (6).Numerical errors are calculated against the cell average given in Appendix A. Figure 4 shows the spatial distribution of the error and the convergence with respect to resolution.These results confirm the error analysis in Eq. ( 10).It could also be remarked that, if the operand of the divergence operator is interpreted as the average along the edge rather than the point value at the edge center, then the Gauss theorem will give the exact cell-averaged divergence without any error.However, it should be noted that in a C-grid discretization, the divergence operator is typically not applied to the horizontal velocity but to the mass flux (cf.Eqs. 3 and 4).Since the mass flux is not a prognostic variable but needs to be derived, an accurate edgemean is not available.In ICON and in many other models the interpolation of density (or equivalent variables) from neighboring cells to edges gives a second order mass flux on a regular grid.It can be shown analytically that if the edge-mean mass flux is approximated to m-th order, m being a positive even number, the divergence operator on an equilateral triangle will be of order m − 1, and the sign of the leading error depends on the orientation of the triangle.(Detailed derivation can be found in Appendix B.4 of Wan (2009), available from http://www.mpimet.mpg.de/en/science/publications/reports-on-earth-system-science.html.) In summary, the truncation error analysis shows that the divergence operator defined on the triangular C-grid yields a checkerboard error pattern.This appears to be an inherent property related to the cell shape and the placement of the velocity variables.The curl and gradient operators, on the other hand, are second-order accurate on the regular planar grid due to the symmetric geometry.The derivations are omitted in this paper.
Following the general procedure of studying a complex problem in a simpler but relevant context, the truncation error analysis presented here is performed on a planar grid consisting equilateral triangles.In an actual model built on a spherical grid, the spherical geometry and unavoidable grid irregularity introduce more terms to the truncation error.In addition, the divergence operator operates on the discrete velocity field produced by model numerics which contains numerical error.In that case Eqs. ( 7) and ( 10) are no longer accurate.On the other hand, the key features of the triangular C-grid that cause the grid-scale noise in the divergence operator, namely the asymmetric shape and the upward-and downward-pointing directions, stay unchanged.

Noise control
The checkerboard error pattern highlighted by the analysis in Sect.4.2 enters the hydrostatic model system because of the continuity Eq. ( 3) and the temperature advection term in Eq. (2) (cf. the discrete form in Sect.5.5).Grid-scale noise in the divergence operator typically causes noise in the divergence field and in temperature, which then affects the velocity field through the pressure gradient force.Such numerical noise, while less apparent in the barotropic tests reported in BR05 and Rípodas et al. (2009), significantly affects threedimensional baroclinic simulations.One possibility to avoid this problem is to improve the divergence operator by enlarging the stencil.For example, with a four-cell stencil (one triangle plus three of its direct neighbors, involving nine velocity points), one can construct a second-order divergence operator on a regular grid, and apply certain measures to approach second order in spherical geometry.Further discussions in this direction will be included in Part 2 of the paper.
Here we only point out that, if the transport of tracers in the same model is to be handled by finite volume methods using a single cell as control volume, a different stencil for divergence in the dynamical core will destroy the tracer-and-airmass consistency whose importance has been pointed out by, e.g.Lin and Rood (1996), Jöckel et al. (2001), Gross et al. (2002) and Zhang et al. (2008).In this paper, we rely instead on a carefully chosen numerical diffusion to suppress the checkerboard noise.Although it is generally not a preferred practice to use filtering or damping algorithms in numerical models, there is an interesting relationship between the discrete divergence and vector Laplacian on the triangular C-grid that can be exploited.
As in the ICON shallow water model described in Rípodas et al. (2009), the vector Laplacian is approximated in an intuitive manner by Here the subscript d denotes the discrete approximation, e the edge midpoint, and N e the unit normal vector associated to the edge.(The relationship between the unit normal vector N e of an edge and the outward unit vector n e of a cell that the edge belongs to is either N e = n e or N e = −n e .)The notations div, curl, grad n and grad τ denote the discrete divergence, curl, normal gradient, and tangential gradient operators defined in Sect.4.1.The divergence operator on a regular planar grid is defined by Eq. ( 6) in Sect.4.2.The fourth-order hyper-Laplacian, widely used for horizontal diffusion in atmospheric models because of its scale selectivity, is approximated by Like in the previous subsection, one can perform a Taylor expansion (but with respect to the edge midpoint e) and get The function H is the same as in Eq. ( 7).
Assuming that the dynamical core uses a diffusion coefficient K 4 and time step t, if Eq. ( 16) is multiplied by (− t K 4 ), applying the divergence operator, and retaining only two leading terms on the right-hand side yields Now we apply Eq. ( 7) to the second term on the right-hand side of Eq. ( 17), use a symbol E div,1 to denote the first order term in Eq. ( 7), in other words and let l = l/ √ 3 to denote the distance between neighboring cell centers.After some manipulation, Eq. ( 17) becomes This shows that, when the vector biharmonic operator (Eq.15) is used in the explicit numerical diffusion for horizontal wind, apart from the hyper diffusion one usually expects (i.e. the second term on the right-hand side of Eq. 19), there is an additional effect that compensates (at least to some extent) the leading error in divergence.This additional effect is similar to the divergence damping mechanism that has been adopted in many dynamical cores (e.g.Lin, 2004).Furthermore, if the coefficient K 4 is determined via a parameter τ * using the formula then t/τ * indicates the fraction of the checkerboard divergence error that can be removed after one time step.
In our numerical experiments it has been observed that the ratio t/τ * = 1 is very effective in removing the grid-scale noise and renders a reasonably stable model configuration.Furthermore, given that the icosahedral grid is not strictly regular in terms of cell sizes and edge lengths, we use the local edge length for Eq. ( 20) instead of the global mean, which further improves the effectiveness of the method.Examples in the shallow water model can be found in (Rípodas et al., 2009).It should be mentioned, however, that this approach does have a clear disadvantage.The two terms on the right-hand side of Eq. ( 19) are controlled by the same K 4 coefficient.Because the value of K 4 has to be chosen to achieve a sufficient compensation of the checkerboard error, there is no longer the freedom to choose the magnitude of the hyperdiffusion by physical arguments only.Since the characteristic damping time τ * corresponding to t/τ * = 1 is considerably shorter than usually seen in climate models, there is a danger of overly strong diffusivity in the triangular ICO-HDC.This is in fact our major concern regarding the viability of this dynamical core in long-term climate simulations, and a point to which special attention needs to be paid in the further development.

Discrete formulation of the dynamical core
We introduce now the discrete form of the primitive Eqs. ( 1)-(4) employed in the ICOHDC.

Horizontal interpolation
On a staggered horizontal grid, the normal wind and the relative vorticity are not co-located with mass (and temperature).Horizontal interpolation is thus necessary.The following interpolations are used in the ICOHDC: ψ c2e , distance-based linear interpolation of a scalar ψ from two neighboring cell centers to the midpoint of the shared edge; ψ v2e , linear interpolation of a scalar ψ from two vertices of an edge to its midpoint (i.e.arithmetic average); ψ e2c,lin , a bilinear interpolation from the three edges of a triangle to its circumcenter.The interpolation is performed in a local spherical coordinate whose equator and primal meridian intersect at the cell center; ψ e2c,aw , an area-weighted interpolation where A c is the cell area, and A c,e the area of a subtriangle defined by the cell center c and the two vertices of edge e (Fig. 5a).This interpolation is constructed from the finite-volume perspective for conservation purposes, assuming piecewise constant sub-grid distribution.The ψ e in Eq. ( 21) is assumed to represent the average value of a kite-like area whose diagonals are the edge e and its dual.ψ e2c,aw is to be understood as the average over the triangular cell c.A c,e in Eq. ( 21) and Fig. 5 is the overlapping area of cell c and the kite associated with edge e.

Vorticity flux and kinetic energy gradient
The first two terms on the right-hand side of Eq. ( 1), i.e. the absolute vorticity flux and the kinetic energy gradient, represent the combination of horizontal momentum advection and Coriolis force in vector invariant form.The kinetic energy gradient is calculated in our model by and the absolute vorticity flux by As mentioned earlier, the C-grid discretization predicts only the normal component of the horizontal wind.The tangential wind v τ needed by Eqs. ( 22) and ( 23) is reconstructed at edge midpoints using the vector radial basis functions (RBFs) introduced in Narcowich and Ward (1994).We use a stencil that involves four edges surrounding the target edge (Fig. 5b), the inverse multi-quadric kernel ψ(r) = 1/ 1 + ( r) 2 , and the shape parameter = 2. Test results have shown that this particular stencil used in our model is very insensitive to the choice of kernel function and shape parameter.Following BR05, it is assumed that the normal and tangential components form a right-hand system.For a more detailed description and testing of this algorithm, the readers are referred to the work of Ruppert (2007) and Bonaventura et al. (2011).Here we only point out that the RBF vector reconstructions have the nice property that they allow for straightforward extension to larger stencils and higher order approximations, as shown, e.g., by Bonaventura et al. (2011).However, this method may have numerical issues related to the ill conditioning of the interpolation matrix when the stencil is relatively large or in the case of highly irregular node distribution.In those cases, enlarging the shape parameter can help alleviate the problem.For the reconstruction of edge-based tangential velocity discussed above, because of the small stencil and thanks to the quasi-regularity of the icosahedral mesh employed, the ill conditioning problem does not arise.Nevertheless there are other reconstruction methods available in the literature that are potentially attractive due to their mimetic properties (e.g., Perot, 2000;Thuburn et al., 2009;Wang et al., 2011).

Pressure and layer thickness
The η coordinate of Simmons and Strüfing (1981), a terrain following coordinate near the earth's surface that gradually transforms into pressure coordinate in the upper troposphere, has been widely used in atmospheric GCMs.Here we only mention a few technical details for completeness and clarity: the pressure at layer interfaces (see Fig. 1) is given by Here p s stands for surface pressure.NLEV is the total number of vertical layers.A and B are predefined parameters (see, e.g.Roeckner et al., 2003).B = ∂p/∂p s is used in Eq. ( 4).The pressure thickness of the k-th model layer is denoted by p k = p k+1/2 − p k−1/2 .

Continuity equation
To compute the right-hand side of Eq. ( 3) the divergence operator is applied to the mass flux v ∂p/∂η , followed by an integral through the vertical column.This discretization does not introduce any spurious sources or sinks in the total air mass, as long as the mass flux has a unique value at each edge.In this paper, the air mass flux in the normal direction N e of an edge e is computed by v n (∂p/∂η) c2e . (25)

Horizontal advection of temperature
The horizontal advection of temperature at cell c in layer k is discretized in an energy-conserving form The mass flux divergence in this equation is the same as in the discrete continuity equation.The heat flux divergence is calculated by first interpolating temperature and layer thickness separately from cells to edges using the distance-based linear interpolation, multiplying by the normal wind, and then applying the discrete divergence operator.Because of the rather simple flux calculation and the inherent property of the divergence operator discussed earlier, the temperature advection is only first-order accurate.The limitation of this low-order scheme is discussed in Sect.6.1.2.

Vertical advection of momentum and temperature
The vertical advection terms in Eqs. ( 1) and ( 2) are discretized in the same way as in ECHAM following Simmons Here ψ is either temperature or a horizontal wind component.The vertical indices are illustrated in Fig. 1.The vertical velocities at layer interfaces are diagnosed by Eq. ( 4).Note that Eq. ( 27) can be derived by applying the idea of Eq. ( 26) to the vertical direction, then replacing the divergence operator by the central difference, and the cell-to-edge interpolation by the arithmetic average.For the advection of the normal wind in the ICOHDC, Eq. ( 27) requires layer thickness and vertical velocity at edge midpoints.These are obtained by the linear interpolation mentioned earlier in Sect.5.4.

Hydrostatic equation
The hydrostatic equation involves only vertical discretization, thus takes the same form in the ICOHDC as in SB81 and in the spectral core of ECHAM.The discrete counterpart of Eq. ( 5) reads for layer interfaces with k 1.Here φ c,s denotes the surface geopotential at cell c.The geopotential at full levels are given by where (30)

Pressure gradient force and adiabatic heating
After computing the geopotential by Eq. ( 29), the last term in the velocity Eq. ( 1) can be obtained by applying the normal gradient operator.
The pressure gradient term in the same equation is calculated by Using Eq. ( 31), the first part of the adiabatic heating term in the temperature equation ( 2) can be obtained by The remaining part is approximated by Both Eqs. ( 32) and ( 33) are derived following the energy conservation constraint in SB81.

A few remarks on the spatial discretization
In this section we have mentioned repeatedly the C-grid discretization of SB81, which is an energy-conserving scheme on regular latitude-longitude grids, designed for an early finite-difference version of the NWP model of the European Centre for Medium-range Weather Forecasts.In the baseline version of the ICOHDC, the kinetic energy gradient and the absolute vorticity flux discretized on the triangular grid (Eqs.22 and 23) do not guarantee energy conservation to machine precision, partly due to the tangential wind reconstruction using the RBFs.The other discrete formulae described in Sects.5.4-5.8, on the other hand, do help avoid spurious energy sources/sinks.Potentially, there might be another issue related to the RBF reconstruction.As discussed in Hollingsworth et al. (1983) and Gassmann (2013), certain inconsistencies between the discrete vorticity flux and kinetic energy gradient can trigger nonlinear instability that manifests itself as small-scale noise at high resolutions.So far we have not seen clear evidences of such instability in the test results of the ICOHDC (cf.Sect.6).

Time stepping scheme
In the baseline version of the ICOHDC we use a time stepping method similar to that of the ECHAM dynamical core.This consists of the leapfrog scheme, the Asselin (1972) filter, and the widely used semi-implicit correction scheme for linear gravity waves (see, e.g.SB81).The reference atmosphere used by the semi-implicit correction is isothermal (T r = 300 K), at rest (v r = 0) on a flat surface (φ r s = 0) with constant surface pressure (p r s = 800 hPa).To achieve computational efficiency, the resulting 3-D Helmholtz equation of divergence is decomposed into a series of 2-D equations, each corresponding to one vertical mode.The 2-D equations associated with phase speed higher than 30 m s −1 are numerically solved using the generalized minimal residual method (GMRES, Saad and Schultz, 1986).
In the standard model configuration, the time stepping scheme uses the following parameters: the Asselin coefficient is 0.1 following ECHAM.In the semi-implicit correction, coefficients of the gravity wave terms evaluated at the old and new time steps are set to 0.3 and 0.7, respectively, following the global forecast model GME of the German Weather Service (Majewski et al., 2002).Detailed formulation of the time integration algorithm is given in Appendix B.
As mentioned earlier in the introduction, this particular time stepping scheme is used here for the sake of a clean evaluation of the spatial discretization employed on the triangular icosahedral grids.The GMRES solver is chosen for its reliability and stability, since this is one of the standard solvers for large non-symmetric linear systems.Issues arising in parallel versions of this algorithm have been discussed in the literature by, e.g., de Sturler and van der Vorst (1995).In massively parallel simulations, especially when the dynamical core consumes a significant portion of the total computing time, linear solvers that require global communications will pose constraints on the computational performance.In such cases, explicit time stepping schemes are an option to consider.For example, the ICON nonhydrostatic model uses a two-time-level predictor-corrector scheme, with implicit treatment restricted to the vertically propagating sound waves.In the ICOHDC, we have implemented various explicit or semi-implicit two-time-level integration schemes for the purpose of achieving tracer-and-air-mass consistency, but these are considered beyond the scope of the baseline model.

Tracer transport
A group of upwind, conservative, flux-form semi-Lagrangian transport algorithms are implemented for tracer transport in the ICON triangular models.Options for horizontal advection include the first-order upwind scheme, and a triangular version of the "swept-area" approach following Miura (2007).The latter algorithm is second-order in time and either second-, third-or fourth-order in space, depending on the choice of reconstruction polynomial (linear, quadratic or cubic).The polynomial coefficients are estimated using a conservative weighted least squares reconstruction method (Ollivier-Gooch and Altena, 2002).Transport in the vertical is calculated with the Piecewise Parabolic Method (PPM, third-order, Colella and Woodward, 1984).The optional limiters employed include semi-monotonic and monotonic slope limiters (Barth and Jespersen, 1989) as well as the Flux Corrected Transport approach (FCT, Zalesak, 1979).Consistency between tracer transport and the discrete continuity equation can be achieved when the dynamical core uses a two-time-level time stepping scheme.These transport algorithms in the ICON models have gone through comprehensive testing, e.g., following the proposal of Lauritzen et al. (2012b).Some of the evaluation results are presented in http://www.cgd.ucar.edu/cms/pel/transportworkshop/2011/16-Reinert.pdf.More detailed results of the second-order (linear) horizontal transport can be found in Lauritzen et al., (2013)  1 .
In the present version of the ICOHDC, the discretization of horizontal temperature advection does not yet make use of the semi-Lagrangian transport schemes, but is calculated using Eq. ( 26) with a rather crude approximation of the heat flux at triangle edges (Sect.5.5).The limitation of the latter scheme is discussed in Sect.6.1.2.In the nonhydrostatic dynamical core Zängl et al. (2013), horizontal advection in the thermodynamic equation and the continuity equation is discretized using the edge-based potential temperature and density values estimated with the second-order Miura (2007) scheme.

Idealized dry dynamical core tests
In this section we present results of numerical simulations carried out to evaluate the new dynamical core.The goal here is to find out whether the present version of the ICO-HDC can correctly represent basic processes of the adiabatic atmospheric dynamics, and to analyze the sensitivity of the numerical solutions to horizontal resolution.During the development of the new model, we routinely perform a suite of idealized dry dynamical core tests with different levels of complexity.The simplest ones are 3-D extensions of the widely used shallow water tests 5 and 6 from Williamson et al. (1992).The barotropic cases are in some sense a sanity check for the 3-D formulation of the ICOHDC and its practical implementation in the code.These results can be found in Wan (2009), and are not repeated here.This section focuses on two deterministic baroclinic tests (Sect.6.1) and an idealized dry "climate" test (Sect.6.2).
Here we want to make the remark that model evaluation and inter-comparison is a delicate matter.Due to the complexity of the governing equations, the wide spectrum of dynamical regimes and waves that they support, and the various components that constitute (and affect the properties of) a discrete model, it is very difficult, if possible at all, to select one single metric to summarize the model performance.In this section, especially in the baroclinic wave test, we carry out simulations with both the ICOHDC and the spectral transform dynamical core of ECHAM at multiple resolutions, and compare the results in several ways including qualitative comparisons, difference norms, and dispersion errors.We think this collection of numerical results provides useful information for the climate and weather modelers who face Fig. 6.Evolution of the baroclinic wave in the Jablonowski and Williamson (2006a,b) test case, as shown by the surface pressure (unit: hPa, left column) and 850 hPa temperature (unit: K, right column) simulated by the new dynamical core at R2B5 (70 km) resolution.Note that in the left column, the two upper panels use a different color scale than the lower panels.Further details can be found in Sect.6.1.2. the question of which numerical method better suits their needs and expectations.
All ICOHDC results presented in this section are obtained using revision 6489 of the model code.The vertical grid is fixed at L31, which resolves the atmosphere from the surface to 10 hPa, as commonly used in the tropospheric version of the ECHAM model (e.g., Roeckner et al., 2006).The model time step is set to 600 s at R2B4, and reduced by half when the grid level is increased by one.

Deterministic baroclinic tests
The test case proposed by Jablonowski and Williamson (2006a,b, hereafter JW06) has been widely used in recent years for testing 3-D atmospheric dynamical cores.Inspired by the baroclinic instability theory, the deterministic test con-sists of two parts: a steady state test followed by a baroclinic instability test.

Steady state test
In the first part of the test, the dynamical core is initialized with a zonally symmetric, geostrophically balanced condition specified by analytical functions.Since this initial condition is a steady state solution of the primitive equations, a perfect numerical model would retain the initial state to machine precision.The spectral core of ECHAM can preserve the zonal symmetry in an arbitrarily long integration.Meanwhile, the model state evolves continuously (but very slowly) from the initial state because of the horizontal diffusion.In the ICOHDC, zonal asymmetries are triggered immediately after model initialization due to grid irregularity near the pentagon points (cf.Sect.3), resulting in wavenumber 5 patterns near 26.6 • N/S.Embedded in the dynamically unstable mean state of this test case, the perturbations amplify for more than 10 days, then reach a quasi-equilibrium state after 20 to 30 days (not shown).As the horizontal resolution increases, the magnitude of the numerical errors is reduced and the perturbations evolve less rapidly.These behaviors agree with our expectation, and are similar to the results of the GME model (also built on icosahedral grids) presented in JW06.

Baroclinic wave test
The second part of this test case focuses on the evolution of an idealized baroclinic wave in the Northern Hemisphere, triggered by an analytically specified large-scale perturbation in the wind field.Cyclone-like structures develop in the course of about 10 days, featuring lows and highs in the surface pressure field (Fig. 6, left column) and accompanying fronts in the lower troposphere temperature (Fig. 6, right column).This figure shows the ICOHDC simulation at R2B5 resolution which has an average grid spacing of 70 km between neighboring mass points.The key features of the simulated baroclinic wave evolution, including the slow development of the perturbations in the first 6 days and the subsequent exponential intensification, as well as the magnitude of the closed cells in surface pressure and the fronts in temperature, agree well with the reference solutions given by JW06.(at R2B3) to 17.5 km (at R2B7).The solution obtained on the coarsest grid (R2B3) is of unsatisfactory quality, in that the depressions are too weak, while the spurious perturbations at the rear of the wave train are too strong.This resolution is thus not recommended for future applications of the new dynamical core.The next solution, at R2B4, is significantly improved, although the first two low pressure cells are still somewhat weak, and there is an easily detectable phase lag in the propagation of the wave in comparison with the solution at R2B7.As the grid is further refined, the phase lag gets smaller, and the depressions become deeper.The two runs at the highest resolutions (R2B6 and R2B7) look very similar, and are hardly distinguishable from the reference solutions in JW06 by visual comparison.

(a) Convergence
To quantitatively assess the convergence of these numerical solutions, we follow JW06 and use the l 1 , l 2 and l ∞ differences norms of surface pressure as the metric.In the work of JW06 it was found that differences among solutions from four models using very different discretization methods stopped decreasing once the resolutions increased beyond a certain limit.Based on this observation, the uncertainty in their reference solutions was estimated.The corresponding uncertainties in the difference norms are shown by the yel-low shading in Fig. 8.When the difference norms fall below the uncertainty limit, the solution being tested is considered as having the same quality as the reference solution.
In Fig. 8 the norms of p s differences are shown between the R2B3 to R2B6 simulations and the R2B7 run (upper row), as well as between the ICOHDC simulations and a reference solution in JW06 from the National Center for Atmospheric Research Semi-Lagrangian dynamical core (NCAR SLD, Fig. 8 second row).Difference norms computed against the other reference solutions in JW06 are very similar hence not shown.Regardless of the choice of reference, panels in Fig. 8 clearly indicate a decrease in the difference norms when the horizontal resolution increases.Convergence of the numerical solution is achieved at R2B6.The R2B6 and R2B7 solutions are able to represent the baroclinic wave evolution within the uncertainty in the reference solution.

(b) Phase speed
In Fig. 9 the ICOHDC results are presented side-by-side with simulations from the spectral core at four spectral resolutions that have been used for the full model in various applications.The difference norms of the spectral model results with respect to the NCAR SLD solution are shown in Fig. 10.By   R2B5 are stronger than those of T85 (Fig. 9, third row), the T85 simulation captures the reference solution within the uncertainty while R2B5 does not.According to the snapshots, the errors in the lower-resolution spectral model results are mainly in the strength of the vortices and the spatial gradients, while in the ICOHDC the phase speed is also a major source of numerical error.Phase error is a typical problem associated with dynamical cores using second (or lower) order spatial discretization methods.It is also one of the main disadvantages of such models at medium and low resolutions in comparison with the spectral transform method.It is worth noting that in JW06 the finite-difference model GME has a similar phase problem, while the NCAR finite volume core (Lin, 2004), which uses the third-order piecewise parabolic advection algorithm, does not.Skamarock and Gassmann (2011) showed that in their models, replacing the second-order potential temperature transport by third-order schemes can significantly reduce the magnitude of the phase error in this test case and suppress its growth.In the hydrostatic dynamical core discussed here, horizontal temperature advection is computed using the first-order scheme described in Sect.5.5.A higher-order discretization, e.g., using the transport algorithms outlined in Sect.5.11, will probably help to improve the solution quality at R2B5 and lower resolutions.On the other hand, Fig. 9 also suggests that phase error in the ICOHDC becomes negligible at R2B6 (35 km).Since the ICON models are developed for high-resolution modelling, the phase error is not expected to be an obstacle in those applications of the new model system.

(c) Equivalent resolution in terms of solution quality
Since the ultimate purpose of developing the new ICON models is to use them in operational NWP and climate appli-cations, a natural question one would expect from the potential users, especially from those having been using ECHAM, is the equivalent resolutions between the ICOHDC and the spectral core in terms of solution quality.For reasons discussed in detail in Appendix C, we believe such relationships are difficult to identify a priori, but rather need to be established by comparing results from numerical tests.This is in line with, e.g., work of Williamson (2008a) who discussed the equivalent resolutions between a finite-volume model and a spectral transform model.
Baroclinic wave tests are carried out with the icosahedral grids listed in the left half of Table 2, and at the spectral truncations given without parenthesis in the right half of the table.Qualitative comparison as in Fig. 9 is used to identify resolution pairs (R5B3, T63), (R2B5, T106), etc., that produce similar results.It turns out that for these visually identified pairs, the ratios between the average grid spacing of different icosahedral grids match well with the ratios between the corresponding truncation wavenumbers.We then use this relationship to derive the wavenumbers given in parenthesis, which are not "standard" resolutions of the ECHAM model.
In principle it would be useful to verify the established equivalent resolutions in some quantitative manner, for example by calculating the difference norms of surface pressure in each pair, and comparing them with the uncertainty estimates (the yellow shading in Figs. 8 and 10).At the current stage, however, the difference norms between the mediumand low-resolution pairs would lie outside the uncertainty range unless the phase error in the ICOHDC was reduced.The verification is thus not done in this study.
It is interesting to note that, in the table, the average grid spacings of the ICON grids match well the zonal grid size at 60 • N on the Gaussian grids of the corresponding spectral resolution.This is probably not a coincidence, but a result of the fact that in this test case the baroclinic wave evolves and propagates near this latitude.Nonlinear terms in the primitive equations play a crucial role in the baroclinic instability development.In both models these terms are computed in grid-point space using similar discretization schemes following the work of SB81.It is thus not surprising that the equivalent resolutions we identified turn out to have similar spacing at 60 • N. In a different test case that features dynamical processes confined to, say, the tropics, the conclusions on equivalent resolutions may be different.

Held-Suarez test
After the adiabatic deterministic test cases discussed above, we consider here the dry "climate" experiment proposed by Held and Suarez (1994), in which the dynamical core is forced by Rayleigh damping of horizontal wind in the nearsurface layers as well as relaxation of the temperature field towards a prescribed, north-south and zonally symmetric radiative equilibrium.The original goal of this popular test was to evaluate the zonal-mean climatology obtained from the last 1000 days of a 1200 day simulation.However, a more comprehensive analysis of the inherent low-frequency variability was carried out in Wan et al. (2008), where an ensemble approach was proposed for the evaluation of the results.Here we follow this approach and perform ensembles consisting of 10 independent 300 day integrations.Each integration starts from the JW06 zonally symmetric initial condition with random noise of magnitude 1 m s −1 added to the wind field.(This choice is rather arbitrary.As long as the 10 integrations within an ensemble are independent, the conclusions drawn in this subsection are not affected by the initial condition.)Simulations are performed at resolutions R2B3, R2B4 and R2B5 using the same configurations as in the deterministic test cases.The zonal-mean climate states are diagnosed from the last 100 days of each integration.
Figure 11 presents the ensemble mean model climate at R2B5.Although simple by design, the Held-Suarez test is able to reproduce many realistic features of the global circulation.Baroclinic eddies cause strong poleward heat and momentum transport (Fig. 11d and c, respectively).The heat transport reduces the meridional temperature gradient in comparison to the prescribed radiative equilibrium (Fig. 11b, here in comparison to Fig. 1c in Held and Suarez, 1994).The meridional transport of angular momentum converges in the mid-latitudes, forming a single westerly jet in each hemisphere (Fig. 11a).The core regions of the jets are located near 250 hPa.The maximum time-and zonal-mean zonal wind is about 30 m s −1 .Easterlies appear in the equatorial and polar lower atmosphere, as well as in the tropics near the model top.The baroclinic wave activities concentrate in the midlatitudes, as depicted by the transient eddy kinetic energy and temperature variance (Fig. 11e and f).The single maximum of eddy kinetic energy in each hemisphere appears in the up-per troposphere near 45 • latitude, close to the core region of the westerly jet.Easterlies in the tropics show little variance.In each hemisphere, the maximum temperature variance appears near the earth's surface and extends upward and poleward.A second maximum of smaller magnitude occurs near the tropopause.These features of the simulated circulation agree well with results reported in the literature (e.g.Held and Suarez, 1994;Jablonowski, 1998;Lin, 2004;Wan et al., 2008).
Sensitivity of the ICOHDC results to horizontal resolution is revealed by Figs. 12 and 13.The contour lines show the differences in the ensemble average of the quantities displayed in Fig. 11, while the gray and light-blue shadings indicate where the differences are significant (at 0.05 and 0.01 significance levels, respectively) according to the Kolmogorov-Smirnov test (Press et al., 1992).Comparing R2B3 with R2B5, the increase in horizontal resolution leads to a substantial enhancement of the eddy activity in the mid-latitudes, stronger poleward transport, and consequently higher temperature in the Polar Regions as well as a poleward shift of the westerly jets.The differences between the R2B4 and R2B5 ensembles are much smaller (Figs. 13).Although one can still see enhancement in the eddy activities (Fig. 13d-f) and temperature differences in high altitudes/latitudes regions (Fig. 13b), the discrepancies are generally much smaller than between R2B3 and R2B5.

First results from the aqua-planet experiments
In the previous section we have evaluated the ICOHDC using dry dynamical core tests at various resolutions.On the whole, the new core produces results that agree reasonably well with those from the spectral core of ECHAM, as well as with the reference solutions available in the literature.In these test cases, the grid-scale noise discussed in Sect. 4 is effectively suppressed and has not yet brought obvious detrimental effects.One might argue that when moist processes are included in the model, condensational heating will act as a positive feedback, which will amplify the grid-scale noise and make the model unstable.To find out whether this is the case, we perform aqua-planet simulations following the proposal of Neale and Hoskins (2000).
For this exercise, the ICOHDC is coupled to the cumulus convection, large-scale condensation, turbulent mixing and radiation parameterizations of the ECHAM6 model (Giorgetta et al., 2012).Second-order horizontal diffusion is applied to the three uppermost model layers to enhance horizontal damping.This is a widely used technique in climate models to effectively dissipate upward propagating waves of various scales and avoid spurious reflection of the vertically propagating waves triggered by cumulus convection and other sub-grid processes.Radiative transfer calculation is performed every other hour as in ECHAM.The large-scale horizontal transport of water vapor, cloud liquid and cloud ice is computed using the second-order Miura (2007) scheme with monotonic FCT, while the vertical transport uses the PPM scheme with semi-monotonic slope limiter (Sect. 5.11).We refer to the resulting model configuration as the ICOsahedral Hydrostatic Atmospheric Model (ICOHAM).In the areas with gray and light blue shading, the differences are judged to be significant in the Kolmogorov-Smirnov test at 0.05 and 0.01 significance levels, respectively.Further details can be found in Sect.6.2.
Simulations are performed at R2B4 using the "Control" and "Qobs" sea surface temperature (SST) profiles of Neale and Hoskins (2000).The reference solutions are from ECHAM6 at T63.Both models are integrated for 1200 days using an 8 min time step and 31 vertical levels.The last 800 days are used in the analysis presented in this paper.Here we do not attempt to investigate the convergence of the aqua-planet experiments (APE) from either ICOHAM or ECHAM6, because both models are new, and neither has been tuned at many resolutions.The intention here is rather to have a first look at the main features of the model "cli-mate".A more comprehensive evaluation of the ICOHAM aqua-planet simulations is the topic of a separate paper.By comparing ICOHAM at R2B4 and ECHAM6 at T63, we are not suggesting that they are an equivalent pair.These resolutions are chosen because the T63 is currently the default Fig. 14.Time and zonal mean surface precipitation rate (unit: mm day −1 , solid black lines) simulated by ICOHAM in aqua-planet simulations at R2B4L31 resolution using the "Control" (left) and "Qobs" (right) SST profiles.The contributions from convective (dotted red lines) and large-scale (dashed blue lines) precipitation are also displayed.Further details can be found in Sect.7. Fig. 15.Wavenumber-frequency diagrams of tropical precipitation (meridionally averaged between 10 • S-10 • N) in aqua-planet simulations carried out with ICOHAM at R2B4L31 (left column) and ECHAM6 at T63L31 (right column).The color shading shows the logarithm of the power of the symmetric component of the unnormalized spectra, diagnosed using the methodology of Wheeler and Kiladis (1999).The black lines indicate the dispersion relationships of westward propagating equatorial Rossby waves, eastward propagating Kelvin waves, and inertia-gravity waves that can propagate either westward or eastward.The upper row shows results corresponding to the "Control" SST profile.The lower row corresponds to the "Qobs" profile.Further details can be found in Sect.7. horizontal resolution of ECHAM6, which has also been used in the CMIP5 simulations, while R2B4 is the resolution used by the ICOHDC and ICOHAM developers in the day-to-day routine tests.
The latitudinal variations of the simulated time-and zonalmean surface precipitation rate in ICOHAM are shown in Fig. 14.In the "Control" case the total precipitation rate peaks at the equator and at about 35 • latitudes, with the main contributors being cumulus convection and large-scale condensation, respectively.In the "Qobs" case, which has the same SST at the equator but weaker meridional gradients in the low latitudes, the tropical precipitation features two peaks, and the mid-latitude rainfall shifts slightly poleward.
Next, we follow Williamson (2008a) and consider the equatorial wave propagation characteristics.Figure 15 presents the wavenumber-frequency diagrams of tropical precipitation (meridionally averaged between 10 • S-10 • N), diagnosed using the methodology of Wheeler and Kiladis (1999).The quantity shown is the logarithm of the power of the symmetric component of the unnormalized spectra.It has been shown by Williamson (2008a) that the power of the background spectrum, usually used for normalizing the "raw" spectrum to identify spectral peaks, is sensitive to model resolution.In order not to lose such signals, we choose to show the raw spectra in Fig. 15, which are also meant to serve as a reference for future work.The normalized spectra are shown in Appendix D. The black lines that indicate the dispersion relationships of equatorial Rossby waves, Kelvin waves and inertia-gravity waves are the same as in Fig. 6 of Williamson (2008a).
In both the "Control" and the "Qobs" simulations, the tropical precipitation has higher power at lower frequencies.The Kelvin waves are more evident in the "Control" case (upper row in Fig. 15), while the Rossby waves show the opposite sensitivity.Interestingly, the Rossby waves in ICO-HAM show a clear peak at westward zonal wavenumber 5 in Fig. 15c.A question naturally arises whether this is an imprint of the icosahedral grid.It should be noted not only that the corresponding ECHAM simulation indicates a similar (albeit weaker) peak (Fig. 15d), but also that the APE Atlas has revealed wavenumber 5 features in the global circulation in many models that employ different types of grids and discretization methods (Williamson et al., 2011, Fig. 4.99).So far it is not yet clear whether and to what extent the icosahedral grid imprint is interacting with this mode.Comparing the four panels in Fig. 15, one can see that on the whole, the power of the waves in ICOHAM at R2B4 is comparable to that in ECHAM6 at T63.The differences between the ICOHAM and ECHAM6 results are much smaller than the sensitivity to SST.
In order to have a first look at the impact of horizontal diffusion on the APE simulations in both models, Fig. 16 shows the 250 hPa kinetic energy (KE) spectrum diagnosed from daily instantaneous output of the vorticity and divergence fields.Each curve in the figure is the average of  Rauscher et al. (2012, MPAS model).The spectra in ECHAM generally follow the n −3 slope from wavenumber 10 up to the truncation limit, as a result of empirical tuning of the order and damping time scale of the hyper-diffusion.The spectra in ICOHAM start to deviate from the n −3 slope at about wavenumber 20, qualitatively similar to the behavior of the NCAR finite volume model at 1.9 • ×2.5 • resolution as shown in Lauritzen et al. (2012a).If we follow Skamarock (2011) and define the effective resolution of the triangular model as the point at which the slope of the simulated spectrum becomes steeper than n −3 , then the ICOHAM R2B4 APE simulation has a effective resolution of about 1000 km, translating to 7 x where x is the grid spacing, which seems to fall into the typical range of 6-10 x as pointed out by Skamarock (2011) for models that use C-grid discretization.
Because of the growing interest in high-resolution modeling in recent years, dynamical core developers have been paying more attention to their models' ability to produce the observed transition of the KE spectrum from the n −3 slope in the inertial regime to a n −5/3 slope in the mesoscale regime, occurring at spatial scales of a few hundred kilometers (Nastrom and Gage, 1985;Lindborg, 1999).For example, Evans et al. (2012) showed that the CAM4 spectral element dynamical core, which uses fourth order hyper-viscosity, is able to resolve the transition when the horizontal resolution is increased to 0.125 • .They pointed out that the CAM finite volume dynamical core with second-order divergence damping has a clearly weaker divergent component of the simulated flow, and expect the version with fourth-order damping to behave similarly to the spectral element core.Takahashi et al. (2006) carried out a series of simulations with the spectral model AFES to empirically determine the appropriate relationship between the magnitude of hyper-diffusion and model resolution, aiming at correctly capturing the shape of the KE in both the inertial regime and the mesoscale regime.Their results suggest a scaling of n −3.22 0 (or x 3.22 , where n 0 is the truncation wavenumber, and x the grid spacing) for the diffusion coefficient.In the ICOHDC, the choice of a fourth-order diffusion with damping time equal to time step implies a scaling of x 3 according to Eq. ( 20), close to what is obtained by Takahashi et al. (2006).On the other hand, in terms of the absolute magnitude at each particular resolution, the diffusion in the ICOHDC/ICOHAM is considerably stronger than that in the spectral model ECHAM.The question whether ICOHAM can produce the KE transition when grid spacing is decreased, and if so, what is the critical grid size, remains to be answered by high-resolution simulations in the future.

Conclusions
In this paper we presented and evaluated a hydrostatic atmospheric dynamical core built on spherical triangular grids.The finite-difference discretization is based on the numerical techniques employed in the ICON shallow water model of Bonaventura and Ringler (2005) and Rípodas et al. (2009), as well as the vertical discretization of Simmons and Burridge (1981).The baseline version of the new dynamical core ICO-HDC uses leapfrog time stepping scheme, with additional semi-implicit correction to handle the fastest gravity waves.
The first outcome of this effort is an improved understanding of the numerical properties of the C-grid discretization on triangular grids.Through the truncation error analysis, it is shown that the discrete divergence operator defined on a single cell using the Gauss theorem is of first-order accuracy even when the triangles are equilateral.The leading error changes its sign from one cell to its immediate neighbors.This explains the grid-scale noise encountered in the development of the new dynamical core.In recent years, similar problems on triangular C-grids have been reported and investigated by other modelling groups (e.g.Le Roux et al., 2007;Danilov, 2010;Weller et al., 2012).Our analysis here provides more insight into the origin of the numerical noise from a different perspective.In addition to highlighting the source, the truncation error analysis also provides useful hints for finding a remedy to control the noise through numerical diffusion.Using the fourth-order hyper-diffusion with a timestep and grid-size dependent coefficient, the first order divergence error can be removed after each time step.The associated disadvantages, however, are the loss of freedom to choose the diffusion coefficient by physical argument, and a rather strong damping of the flow.
After the theoretical analysis, the ICOHDC is evaluated using idealized test cases of various complexity.We focused on the deterministic baroclinic instability test of Jablonowski and Williamson (2006a,b), carried out simulations at various horizontal resolutions, and compared the ICOHDC results with those from the spectral dynamical core of ECHAM and the reference solution from a NCAR model.In this test the ICOHDC results show the expected resolution sensitivity, and converge at R2B6 (35 km grid spacing).The R2B6 solution correctly captures the evolution of the dynamical instability, as well as the strong gradients in vorticity and temperature associated with the cyclone-like structures.The R2B5 solution does not yet reach convergence due to a phase error in the baroclinic wave, probably attributable to the loworder discretization method used for temperature advection.Longer, idealized "dry-climate" simulations are performed following the proposals of Held and Suarez (1994) and analyzed with the ensemble technique suggested by Wan et al. (2008), in which the ICOHDC correctly reproduces key features of the meridional heat and momentum transport by baroclinic eddies, and shows a similar resolution sensitivity in comparison to the spectral transform core of ECHAM.
As a first step of testing the interactions between the parameterized moist physics and the numerical schemes implemented in the new dynamical core, aqua-planet simulations (Neale and Hoskins, 2000) are carried out with the ICOHAM model at R2B4 (140 km) and with ECHAM6 at T63 (1.875 • ), both with 31 vertical layers extending from the earth's surface to 10 hPa.The ICOHAM model is able to reproduce the same equatorial wave propagation characteristics as in ECHAM6, including the sensitivity of such characteristics to the meridional SST profile.At this resolution, results from the new model do not show clear evidence of contamination by grid-scale noise.The 250 hPa kinetic energy spectra have less power than in the spectral model from wavenumber 20 onwards, as expected from the stronger diffusion applied.It is not yet clear how these features will change with horizontal In the following we use underlines to denote column vectors (matrices) containing discrete values of a quantity defined at all vertical layers at the same horizontal location, in other words Equations (B15-B18) can be seen as a simplified version of the formulae in the appendix of Simmons and Burridge (1981).In the ICOHDC we have followed ECHAM and chosen an isothermal reference state, whilst the reference temperature in Simmons and Burridge (1981)  where div and grad n are the discrete divergence and normal gradient operators, respectively.The semi-implicit time integration scheme applied in the hydrostatic model can be summarized as the following algorithm: 1. Apply the explicit leapfrog scheme to obtain v n n+1 expl , T n+1 expl and (p s ) n+1 expl .Calculate θ tt,expl D using Eq.(B20).

Appendix C A discussion on equivalent resolutions
In Sect.6.1.2,equivalent resolutions between the ICOHDC and the spectral transform dynamical core of ECHAM are established by comparing results from the baroclinic wave test.Here we make the attempt to explain why we believe such relationships cannot be easily estimated a priori.First we clarify that the concept of equivalent resolutions discussed here is to be understood as resolutions that yield numerical solutions of the same quality.Formally equivalent resolutions, i.e., resolutions that are expected to give equivalent results because of truncation error analysis, can only be defined a priori when the truncation error can be estimated, at least as an order of magnitude.We are not aware of reliable ways of doing this in the case of nonlinear problems.
The total degrees of freedom (DOF), defined as the total number of prognostic variables in the model, is not a good index for solution quality because it does not provide information about the order of accuracy of the whole suite of discretizations applied to the nonlinear governing equations, nor about the impact of other components such as the form and magnitude of numerical diffusion.There is, for example, no obvious reason why the ratio of total DOF between the resolution pairs in Table 2 is about 6 -7 to 1 (ICOHDC to spectral).
The DOF of a certain physical quantity, e.g., mass or velocity, is not necessarily a good metric either, because of the possible existence of redundant DOF and the actual reduction of DOF due to numerical diffusion.In Sect.6.1.2cwe made the comment that between the equivalent resolutions identified by comparing results of the baroclinic wave test (Table 2), the average grid spacings of the triangular grids match the zonal grid size at 60 • N on the Gauss grids of the corresponding spectral resolution.One might consider this as a hint that the number of mass points can be a useful index.However, it is worth noting that this match is observed in this particular test case, and under the condition that both models employ similar finite-difference discretizations for the nonlinear terms in the governing equations.Furthermore, the quasi-uniform distribution of grid points on the icosahedral grids and the clustering on the Gauss grids introduces further complication when using the number of mass points as a global measure when comparing resolutions.
Although one can define an expected effective or equivalent resolution using the abovementioned elements, we believe the actually achieved effective/equivalent resolution needs to be assessed a posteriori by comparing results from numerical simulations.

Fig. 2 .
Fig. 2. Schematic showing the stencils of (a) the divergence and curl operators, and (b) the normal and tangential gradient operators described in Sect.4.1.
Fig. 4. (a) Numerical error of the cell-averaged divergence of the velocity field defined by Eq. (11), calculated using Eq.(6) on a planar triangular grid with 10.4 • resolution in the x-direction.(b) The l 1 , l 2 and l ∞ error norms at different resolutions.The discrete divergence is calculated by first evaluating Eq. (11) at edge centers then applying operator (6).Numerical error is computed with respect to cell average.

Fig. 7 .
Fig. 7. Surface pressure (unit: hPa, left column) and 850 hPa temperature (unit: K, right column) at day 9 in the Jablonowski and Williamson (2006a,b) baroclinic wave test simulated by the new dynamical core at various horizontal resolutions.Further details can be found in Sect.6.1.2.

Figure 7 Fig. 8 .
Figure7shows the same fields as in Fig.6but after 9 days of integration, and at 5 different horizontal resolutions.The average grid spacing between mass points ranges from 280 km

Fig. 9 .
Fig.9.850 hPa relative vorticity (unit: 10 −5 s −1 ) at day 9 in theJablonowski and Williamson (2006a,b)  baroclinic instability test simulated by the new dynamical core (left column) and the spectral transform core of ECHAM (right column) at various horizontal resolutions.Further details can be found in Sect.6.1.2.

Fig. 10 .
Fig. 10.As in Fig. 8 but between simulations performed with the spectral transform dynamical core of ECHAM and the reference solution at T340 provided by the NCAR semi-Lagrangian model.
Figures 12 and 13 together show a clear trend of convergence in the ICOHDC results.

Fig. 11 .
Fig. 11.Zonal mean climate state simulated by the ICOHDC in the Held-Suarez test at R2B5 resolution.The quantities shown are ensemble averages of 10 independent integrations.Each ensemble member starts from the same initial condition but with random noise added to the wind field.Further details can be found in Sect.6.2.

Fig. 12 .
Fig. 12. Differences between the ensemble mean climate statistics in the Held-Suarez tests performed with the ICOHDC at R2B3 and R2B5 resolutions.Dashed contours indicate negative values.In the areas with gray and light blue shading, the differences are judged to be significant in the Kolmogorov-Smirnov test at 0.05 and 0.01 significance levels, respectively.Further details can be found in Sect.6.2.

Fig. 16 .
Fig. 16.250 hPa kinetic energy spectra in the aqua-planet simulations performed with ICOHAM at R2B4L31 resolution and with ECHAM6 at T63L31.The upper panel shows results obtained with the "Control" SST profile.The lower panel corresponds to the "Qobs" profile.The spectra are diagnosed from daily output of instantaneous vorticity and divergence fields.Each curve shown in the figure is an average of 800 snapshots.

Fig. A1 .
Fig. A1.As in Fig. 15 but showing the normalized spectra.Further details can be found in Sect.7.

4.
Substitute the newly obtained T n+1 and p n+1 s into Eq.(B12) to get the normal velocity at the new time step.5. Apply Asselin filter to v n , T and p s .6. Switch time indices: n − 1 ← n , n ← n + 1 , then go to step 1.

Table 1 .
Total number of triangle cells and edges in various grids with root division n r = 2, the average distance between neighboring cells, and the area ratio of largest to smallest triangles.Grid −1 is the icosahedron projected onto the sphere.

k − 1/2 k k + 1/2 p , φ , η ∂p ∂η T v n ζ Fig. 1.
Illustration of the triangular grid and the location of main variables.Vertical level indices are shown to the left of the sketch.The meaning of the symbols can be found in Sects. 2 and 3.

Table 2 .
Resolutions of the ICOHDC and the spectral core of ECHAM that produce similar results in the baroclinic wave test case.The grid size given in the left half of the table is the average distance between mass points on the triangular grids."dx" in the right half of the table refers to the zonal spacing of the Gauss grid.The degrees of freedom (DOF) and the total number of mass points n M are included in the table for the discussion in Appendix C. The DOF of the ICOHDC is defined as the total number of velocity and mass (temperature) points on one vertical level.The DOF of the spectral core is defined as the total number of spectral coefficients of vorticity, divergence, and temperature on one vertical level.The n M in the spectral model is that of the corresponding Gauss grid.