Icosahedral Shallow Water Model (ICOSWM): results of shallow water test cases and sensitivity to model parameters

The Icosahedral Shallow Water Model (ICOSWM) has been a first step in the development of the ICON (acronym for ICOsahedral Nonhydrostatic) models. ICON is a joint project of the Max Planck Institute for Meteorology in Hamburg (MPI-M) and Deutscher Wetterdienst (DWD) for the development of new unified general circulation models for climate modeling and numerical weather forecasting on global or regional domains. A short description of ICOSWM is given. Standard test cases are used to test the performance of ICOSWM. The National Center for Atmospheric Research (NCAR) Spectral Transform Shallow Water Model (STSWM) has been used as reference for test cases without an analytical solution. The sensitivity of the model results to different model parameters is studied. The kinetic energy spectra are calculated and compared to the STSWM spectra. A comparison to the shallow water version of the current operational model GME at DWD is presented. The results presented in this paper use the ICOSWM version at the end of 2008 and are a benchmark for the new options implemented in the development of the ICON project.


Introduction
ICON (acronym for ICOsahedral Nonhydrostatic) is a joint project of the Max Planck Institute for Meteorology (MPI-M) and Deutscher Wetterdienst (DWD), the national weather service of Germany, for the development of new general circulation models. The project aims at unified general circulation models for climate modeling and numerical weather Correspondence to: P. Rípodas (maria-pilar.ripodas@dwd.de) forecasting on global or regional domains. The new model will be based on finite volume (for the continuity equation) and finite difference discretizations of the fully elastic, nonhydrostatic Navier-Stokes equations on geodesic, icosahedral, locally refinable grids. Various research institutes in Germany and elsewhere are also contributing to the project, among which are University of Postdam, Free University of Berlin and Los Alamos National Laboratory. Bonaventura (2004) discussed the current problems in NWP and climate modeling like mass conservation and monotonicity of tracer concentrations, local mesh refinement and the use of massively parallel computers for high resolution modeling. The ICON project joins DWD and MPI-M resources to face these problems in the development of new models.
As a first step, a shallow water model has been developed: the Icosahedral Shallow Water Model (ICOSWM). A first version of ICOSWM has been described in Bonaventura (2003Bonaventura ( , 2004 and . In Sect. 2 the main features of the model are given and the differences between the current version and the previous version in Bonaventura (2003Bonaventura ( , 2004 and  are highlighted. To test the results of ICOSWM, the standard shallow water test suite of Williamson et al. (1992) is considered. In particular, results for test cases 2 (global steady state nonlinear zonal geostrophic flow), 5 (zonal flow over an isolated mountain), and 6 (Rossby-Haurwitz wave) of the standard shallow water test suite are shown in Sect. 3. Convergence of model errors for different grid resolutions is considered. Model results for test cases 5 and 6, which have no analytical solution, are compared to high resolution runs of a variant

The Shallow Water Equations on the sphere
The vector invariant form of the shallow water equations on the sphere is considered here: Here v = (u,v) is the horizontal velocity vector (on the sphere), K = 1 2 (u 2 + v 2 ) is the kinetic energy per unit mass, ζ is the vertical component of the relative vorticity, f is the Coriolis coefficient, h * is fluid depth, h = h * + h s is the height of the free surface, h s is the height of the orography, g is the gravitational constant and k the unit vector in the radial outward direction.

The model grid and the discrete operators
The discretization method employed is defined as a special case of the Delaunay triangulation on the sphere, i.e. the icosahedral geodesic grid described e.g. in Baumgardner and Frederickson (1985). The main reasons for the choice of this type of grid is its quasi-uniform coverage of the sphere, which solves automatically the pole problem of regular latitude-longitude grids. Furthermore its hierarchical structure provides a very natural setting for local grid refinement on nested grid hierarchies. Finite element approaches based on such geodesic grids have been introduced in Cullen (1974), Giraldo (2000), Heinze and Hense (2002). Finite volume approaches were presented in Heikes and Randall (1995a), Ringler et al. (2000), Ringler and Randall (2002).
The icosahedral construction yields a Delaunay triangulation of the sphere to which a Voronoi tessellation is naturally associated (see e.g. Quiang et al. (2003) and the references therein for a complete description of Delaunay-Voronoi grid pairs on the sphere), which consists of convex spherical polygons (either pentagons or hexagons, see Fig. 1). The triangular Delaunay grid is chosen as the primal grid and the pentagon-hexagon Voronoi grid as is the dual grid for ICOSWM.
The mass and vorticity preservation properties in ICON are achieved by use of triangular Delaunay cells on the sphere as control volumes for mass and of the dual Voronoi cells (pentagons or hexagons) as control volumes for vorticity. The orthogonality of the primal and dual grid edges allows the use of simple approximations to the gradient, divergence and curl operators, in the framework of a C-type staggering of the discrete variables. This represents a major change with respect to the discretization employed e.g. in GME (operational global model of Deutscher Wetterdienst, Majewski et al., 2002), where an A grid approach was used and discrete variables were defined at the vertices of the Delaunay grid and the orthogonality of primal and dual grids was not exploited.
In order to develop an analog of the rectangular C-type staggering (see e.g. Arakawa and Lamb, 1981;Lin and Rood, 1997;Ringler and Randall, 2002;Sadourny, 1975) on the Delaunay grids, the mass points are defined as the circumcenters of the triangular grid cells, while the velocity points are defined for each cell edge as the intersection between the edges of the Voronoi and Delaunay cells (see Fig. 1). By construction, each of these points is equidistant from the centers of the Voronoi cells at the ends of that edge. In Fig. 2 i is a mass point, l is a velocity point and v is a vorticity point.
In a C grid discretization approach, the discrete prognostic variables considered are the value of the height field at the mass points (i), interpreted as a cell averaged value, and the velocity components normal to the triangle edges at the edge midpoints (l). The tangential velocity components, which are needed for the computation of the Coriolis force term, must be reconstructed.
The projection of the regular icosahedron on the sphere yields the so-called base grid or grid level −1. Two of the vertices of the icosahedron coincide with the poles of the reference spherical coordinate system. A first refinement step for which edges are divided in 2 or more, generally in n equal arcs, and connected by great circle arcs "parallel" to the edges of a parent cell, then results in the so called root grid, or grid level 0. Hence each cell of the base grid is divided in n 2 new triangular cells, or 4 cells if the original triangle edges are divided in two equal sections. From here on the grid construction allows only repeated bisection of triangular cell edges, yielding a hierarchy of computational grids numbered as grid levels 1, 2, etc. The number of cells quadruples at each refinement step. Note that this numbering of grid levels is different from that used in Bonaventura (2003Bonaventura ( , 2004 and . Figure 1 shows the triangles in red and hexagons/pentagons in blue of the root grid, or level 0 grid, resulting from an initial dyadic refinement step. Table 1 shows the number of mass points, velocity points and vorticity points of the different grid levels, again for a dyadic refinement of the base grid. Before introducing the discrete operators, some notation to describe the grid topology and geometry will now be introduced. Let i denote the generic cell of the Delaunay grid. Let E(i) then denote the set of all edges of cell i. The grid point associated to cell i will also be referred to as the cell center. The generic vertex of a cell, which is also the center of a cell in the dual grid, is denoted by v. E(v) denotes the set of all edges of the dual cell whose center is vertex v. The area of cell i is denoted by A i , while the area of the dual cell is denoted by A v . Let then l denote the generic edge of  −1  20  30  12  0  80  120  42  1  320  480  162  2  1280  1920  642  3  5120  7680  2562  4  20 480  30 720  10 242  5  81 920  122 880  40 962  6 327 680 491 520 163 842 a cell. It is to be remarked that this index can be assigned at the same time to an edge of the primal grid and the edge of the dual grid, which by construction intersects the primal grid edge at its midpoint. The number of edges is actually equal for both grids. The length of the edge l of a cell is denoted by λ l and the distance between the centers of the cells adjacent to edge l (i.e., the length of a edge of the dual cell) is denoted by δ l . At each edge, a unit vector N l normal to the edge l is assigned. T l denotes the unit vector tangential to the edge l, chosen in such a way that N l × k l = T l holds, where k l denotes the radial outgoing unit vector perpendicular to the tangent plane at the intersection of primal and dual edge l. Furthermore, for each cell edge, the unit vector pointing in the outer normal direction with respect to cell i is denoted by n i,l . Unit vectors n v,l are also introduced, as pointing in the outer normal direction with respect to the dual cell v. The corresponding tangential vectors t v,l are defined so that n v,l × t v,l = k l . Given the edge l of a cell, the adjacent cells are denoted by the indexes i(l,1) and i(l,2), respectively. The indexes are chosen so that the direction from i(l,1) to i(l,2) is the positive direction of the normal vector N l . Vertex indexes v(l,1) and v(l,2) can also be defined analogously, so that the direction from v(l,1) to v(l,2) is the positive direction of the vector T l . Discrete divergence and curl operators are now introduced in the context of the C grid staggering outlined above. Given a generic discrete vector field G on the sphere, its value at a velocity point can be represented as G l = g l N l +ĝ l T l , where g l ,ĝ l denote the normal and the tangential components, respectively. The operators are defined as acting on the set of values g l assigned at the edges of the Voronoi-Delaunay grid. The discrete divergence and curl operator can be naturally defined as follows: The Voronoi-Delaunay property of the grid (the normal direction to the edges of the primal grid cell is the tangential direction to the corresponding edge of the dual grid cell) has been used here in an essential way. By the same property, the discrete normal and tangential derivatives can also be approximated as where φ i , ψ v are discrete functions defined on the primal and dual grid cells, respectively. Since the velocity points are not equidistant from the adjacent Delaunay triangular grid cell centers, the difference operators described above are only first order accurate. However, grid optimization procedures can partly cure this problem by reducing the off-centering to rather small values. The grid generator for the ICOSWM model has several optimization options. In this paper only results for a grid optimized with the method suggested by Heikes and Randall (1995b) are shown.
In Table 2 the minimum, mean and maximum distances between mass points and between vorticity points for grid levels 0 to 6 for an Earth radius of 6.371229×10 6 m are shown for the case of a Heikes-Randall optimization of the icosahedral grid. The off-centering (%) of a velocity point is defined as where d vel-mass is the distance between the velocity point and one of the adjacent mass points and d 2massp is the distance between the two adjacent mass points. In Table 2 the maximum of the velocity point off-centerings for each grid level is also given. The off-centering is reduced to rather small values for the higher grid levels. In Bonaventura (2003), Bonaventura (2004) and Bonaventura and Ringler (2005), the off-centering is defined as twice the value in Eq. (7).

Reconstruction of a vector field from the normal components
In order to recover the full velocity vector from the normal velocity components prescribed at the velocity points in a C grid variable staggering, a reconstruction procedure is needed. This is essential for the discretization of the shallow water equations, especially for the representation of the Coriolis force terms. We will always be concerned here with a vector field that is reconstructed at the triangular cell centers and whose normal components are assumed to be known at the edges of the triangles. Two options are available in the ICOSWM model for the reconstruction, the Raviart-Thomas element of order 0 (RT0) and a Radial Basis Function reconstruction (RBF). In Bonaventura (2003Bonaventura ( , 2004 and , only the Raviart-Thomas reconstruction was available and used. The Raviart-Thomas technique was introduced in Raviart and Thomas (1977) and a complete description of the mathematical properties can be found in Quarteroni and Valli (1994).
In the following, the RBF vector reconstruction implemented in ICOSWM (Ruppert, 2007) is described. Radial basis functions (RBFs) are powerful tools for interpolating scattered data (Narcowich and Ward, 1994). In ICOSWM, they are used for providing accurate approximations to a two-dimensional vector field on the sphere obtained from its components normal to the edges of the icosahedral triangular grid. The interpolation is performed in a 3D Cartesian coordinate rotating with the Earth, whose origin is located at the Earth's center, and the z-axis aligns with the Earth's rotation axis and points to the north. In the following, any point on the sphere of radius a is indicated by the location vector x = (x,y,z).
The interpolation problem can be described as follows: given an arbitrary destination point x 0 , define a region surrounding it. In this region at location x j (j = 1,...,m), the projection of the vector field of interest (denoted by v) in the direction of the unit vector N j is know as v j · N j = v n j . Assume that within this region the vector field can be approximated by the function Here φ is a radial basis function whose value depends only on the great circle distance between x and x j , i.e., Commonly used RBFs include The parameter ε is often called the scale factor and defines a kind of influence radius of the RBF. In ICOSWM, we use RBFs that are monotonically decreasing with respect to the distance r. From Eqs. (8) and (9) one can derived a linear algebraic equation in which Formally solving Eq. (12) for c and substituting the solution into Eq. (8), we get in which For a particular application, once the stencil and the form of the RBF have been determined, the matrix −1 and vector φ 0 can be calculated and stored. Each time when the vector reconstruction is needed, the approximation can be obtained by applying Eq. 14.
In the ICON models, the RBFs are used for reconstructing the horizontal wind vector at triangle centers from the normal wind given at some triangle edges in the vicinity. Three stencils of different sizes have been implemented (Fig. 3). Numerical tests (Ruppert, 2007) have shown that the 3-point stencil produces in general first-order approximations, the 9point stencil second-order and the 15-point stencil third order accuracy. The influence of the RBF kernels, scale factor and stencil used in the reconstruction on the performance of ICOSWM is tested and shown in Sect. 3.

Spatial discretization
The discrete operators introduced in Sect. 2.2 are used to define the spatial discretization. The discrete normal component of the velocity with respect to a triangle edge will be denoted by u l , while the corresponding discrete tangential component will be denoted by v l . The associated discrete vector field will be denoted by U l = u l N l + v l T l , The spatial discretization of the continuity Eq. (1) is straightforward by integration on cell i and application of the divergence theorem: whereh * l denotes an average of the layer thickness values in the neighboring cells. The resulting numerical method conserves mass by construction. This is important as discrete conservation properties have long been identified as an important feature of global circulation models (see e.g. Arakawa and Lamb, 1981;Lin and Rood, 1997;Ringler and Randall, 2002;Sadourny, 1975).
The discrete momentum equation can be derived by taking the scalar product of Eq. (2) with the unit vector N l at a generic velocity point. Using the vector identity and the definitions given in the previous section yields the equation Here, v l is an approximation of the tangential velocity component, η v = ζ v + f v , where ζ v = curl(u) v andη l is an average of the absolute vorticity values at the ends of the cell edge. The tangential velocity at the edge is obtained as the projection in the tangential direction of the (inverse distance weighted) average of the reconstructed wind vectors at the cell centers adjacent to the edge. In Bonaventura and Ringler (2005), either potential enstrophy conserving or total energy conserving variants of the same method were proposed. Equation (20) in Bonaventura and Ringler (2005) specify how to calculate the edge average of the absolute vorticity in Eq. (17) in order to conserve potential enstrophy. Furthermore, in Bonaventura and Ringler (2005) a simpler formulation was also introduced, which is essentially equivalent to the potential enstrophy preserving scheme and produces indistinguishable results. In the latter formulation, the edge averaged value of the absolute vorticity is obtained by simple arithmetic average of the values of η v at the neighboring vertices. In the present paper this simpler formulation is employed, combined to the more accurate RBF reconstruction procedure described in the previous section.

Three-time-level semi-implicit time discretization
A three time level semi-implicit time discretization of Eqs. (17)-(16) based on the leapfrog scheme is given by Here, φ n+ 1 2 = φ n+1 + φ n−1 /2, h * l = h n l −h s l , and t is the time step.
The detailed outline of the practical implementation of the three time level semi-implicit discretization method is the following: -Substitute the expression for u n+1 l into Eq. (19). For each cell i, one obtains the discrete wave equation where all the explictly computed values are collected in the terms − 2 t f l +ζ n l v n l − t δ ν gh n−1 l − 2 t δ ν K n l .
-once h n+1 has been computed, it is back substituted in Eq. (18) to obtain final update of u n+1 .
The set of all Eq. (20) for each cell i yields a linear system in the unknowns h n+1 i . Its matrix is sparse, symmetric, positive definite and diagonally dominant, which allows for efficient solution even when using relatively simple solvers.
Once the values of h n+1 i have been computed, they are back substituted in Eq. (18) to obtain the final update of the discrete velocities.
Asselin time filtering (see e.g. Asselin, 1972) has to be applied to filter computational modes of the leapfrog discretization, so that quantities at time level n are filtered as follows: where is a coefficient independent of the time step and of the resolution. Results of ICOSWM are shown in Sect. 3 for different values of (0.2-0.03).
In this paper all the results shown are for the threetime level semi-implicit scheme, while in Bonaventura (2003Bonaventura ( , 2004 and  a two-time-level semi-implicit time discretization is used. The three-timelevel scheme is computationally more efficient and makes ICOSWM more comparable to the NCAR STSWM and GMESWM, that also employ a three-time-level scheme.

Numerical diffusion
An explicit diffusion term can be added to the right hand side of the prognostic variable equations to remove small-scale noise and improve the stability of the numerical scheme. Due to non-linear interactions, there is a transport of energy to small scales and consequently the kinetic energy spectrum grows at the small scales. To avoid it, a linear fourth-order diffusion is applied optionally to the velocity field. There is also an option to apply diffusion to the mass field, but it is not used here because it is considered unphysical as the continuity equation does not contain a turbulent mixing term in the full atmospheric equations. In addition, it could lead to a violation of mass conservation for our implementation of the diffusion coefficient.
To implement the numerical diffusion, the discrete vector Laplacian must be defined at the velocity points. It is based on the relation ∇ 2 v = ∇ (∇ · v)−∇ ×(∇ × v) and involves all the basic discrete operators defined in Sect. 2.2.
The discrete vector Laplacian in ICOSWM of the discrete velocity field U l , is defined as For a fourth-order diffusion, with the form of −k 4 ∇ 4 U, where k 4 is the so called diffusion coefficient, the discrete vector Laplacian operator is applied twice.
In ICOSWM, Eq. (23) is used for the relation between the diffusion coefficient and a characteristic damping time τ in which the grid scale noise is removed, where dgl is the dual grid length or distance between mass points and pgl is the primal grid length or distance between vorticity points. A uniform planar grid of equilateral triangles is considered to derive Eq. (23) (Wan, 2009).
The ICOSWM grid is not homogeneous, the value of the primal and dual grid length varies from point to point, and therefore a location dependent diffusion coefficient is used.
The characteristic damping time τ is a good parameter to represent the amount of diffusion applied in ICOSWM, although in the literature usually is the value of k 4 what is documented. Table 3 shows the minimum, mean and maximum values of the diffusion coefficients used in the model for different grid levels for a characteristic damping time of 2 h. The corresponding diffusion coefficients for a different characteristic damping time parameter τ can be obtained multiplying the values in Table 3 by the factor (2 h/τ ).
The inhomogeneities of the grid edge lengths, which increase with resolution, translate to ratios between maximum to minimum diffusion coefficients of ∼2.7 to ∼3.4 for grid levels 2 to 6, respectively.

Results of shallow water test cases
The standard shallow water test suite of Williamson et al. (1992) is a very useful benchmark for the model development process. This test suite comprises a number of idealized tests which are representative of some main features of large scale atmospheric motion. This section presents results for the steady state zonal geostrophic flow (test case 2), the zonal flow over an isolated mountain (test case 5) and the Rossby-Haurwitz Wave (test case 6). The time steps are set to 1440, 720, 360, 180, and 90 s for grid levels 2, 3, 4, 5, and 6, respectively to yield similar Courant numbers at different grid levels. All the results presented here are obtained with the semi-implicit three-time level scheme and Heikes-Randall optimized grid.
The normalized errors l 2 , and l ∞ in Williamson et al. (1992) are used to test the model quantitatively. For the case of the height field, the expressions for the l 2 , and l ∞ errors are where λ and θ are the longitude and latitude of the grid points, h is the model output, h T is the true solution if there is an analytical solution and a reference solution if not, and I is a discrete approximation to the global integral

Williamson's test case 2 with zonal flow
Test case 2 of the standard shallow water suite of Williamson et al. (1992) is a steady state solution of the non-linear shallow water equations. It consists of a solid body rotation with a balanced geostrophic height field. The spherical coordinate poles are not necessarily coincident with earth rotation axis. We denote α the angle between the coordinate and the rotational axis. We consider here only the case α = 0: the poles of the coordinate system, that are two of the vertices of the original icosahedron, coincide then with the rotation axis. For this test case, an analytical solution is available, so that approximate convergence rates can be computed by applying the numerical method at different resolutions.
Convergence results for different sets of model parameters after 10 days simulation are presented below.

Sensitivity to the characteristic damping time
We start to test the influence of using different characteristic damping times for the numerical diffusion. Several experiments have been run with different characteristic damping times (28, 12, 2, and 1 h) and also without diffusion (τ =∞).
The tests have been performed with Asselin filter parameter 0.1 and g9p-0.5 RBF reconstruction.
A comparison of the convergence of the normalized errors for the height and vorticity fields at day 10 for the different experiments is shown in Fig. 4. The results for the different characteristic damping times are plotted in different colours and in the order shown in the figure. As in the other convergence figures shown in this work, the black line represents a second order convergence and the l 2 and l ∞ normalized errors are represented with solid and dash-dotted lines respectively. The behavior of the l 1 errors (Williamson et al., 1992) is similar to the l 2 errors and it is not plotted.
The convergence plot for the wind field is not shown because no significant sensitivity to the characteristic damping time is observed in the normalized errors for this field. A second order convergence is observed for both l 2 and l ∞ (see Table 4).
In Fig. 4 (bottom) we observe that the normalized errors for the height field increase with decreasing characteristic damping times (increasing diffusion). The effect is larger in the case of the l ∞ error. A second order convergence is observed for both l 2 and l ∞ .
In Fig. 4 (top) we observe the opposite effect in the case of the vorticity field. The errors are reduced with increasing diffusion coefficients (decreasing characteristic damping times), especially for the higher grid levels. In a geostrophic balance, the vorticity is proportional to the Laplacian of the height field, and therefore any small scale noise present in the height field is amplified in the vorticity field. Increasing diffusion reduces this noise, reducing at the same time the normalized errors for the vorticity field. For the vorticity field second order convergence is only achieved for the higher grid levels when high diffusion coefficients are used (characteristic damping times of the order of 2 h). The positive effect of a larger diffusion in the vorticity errors is larger than the negative effect in the height errors. Therefore a 2 h characteristic damping time seems to be a good choice. Smaller characteristic damping times are not recommended because the experiment with τ =1 h has larger errors for all the variables than the experiment with 2 h characteristic damping time. In the case of τ =1 h, the diffusion is smoothing too much. Table 4 presents the numerical values of the normalized l 2 and l ∞ errors after 10 days for the different variables from grid level 2 to 6 for the experiment with 2 h characteristic damping time. Figure 5 shows the error fields for height (bottom) and vorticity (top) for τ =∞, i.e. no explicit diffusion, g9p-0.5 RBF reconstruction, Asselin parameter 0.1 and grid level 6. The height errors show a clear wavenumber-5 pattern due to the icosahedral grid.
In the vorticity errors a wavenumber-5 pattern can be observed, together with small scale noise. Error spikes near the special points (vertices of the original icosahedron) are a consequence of the irregularity of the grid. Applying explicit horizontal diffusion, the small scale noise and the effects of the irregularity of the grid are reduced, decreasing the l 2 and l ∞ vorticity errors.
This is a good test case to test the effect of the optimization of the ICOSWM grid, because the errors are small and show a clear wavenumber-5 pattern due to the grid. Two experiments have been run with the non-optimized grid, with τ =∞ and τ =2 h and the same parameters as in the previous experiments.
When no numerical diffusion is applied, and the grid is not optimized, the l 2 and l ∞ errors for the height field and the l 2 errors for the wind field after 10 days show the same convergence rate as in the case of the optimized grid, but the errors are smaller than in the case of the optimized grid. As an example, l 2 height errors after 10 days and grid level 6 are 0.5483e-5 and 0.6955e-5 for the non-optimized and optimized grids. The l ∞ errors for the wind field are slightly larger in the case of the non-optimized grid. But the l 2 and l ∞ errors for the vorticity field after 10 days are definitely larger for the non-optimized grid. l 2 vorticity errors after 10 days and grid level 6 are 0.3164e-2 and 0.1244e-2 for the non-optimized and optimized grids. The l ∞ vorticity error for grid level 6 is also ∼3 times larger in the non-optimized grid. Plotting the errors for the vorticity (not shown), a wavenumber-5 pattern together with small scale noise is observed. The shape of the original icosahedron is more visible than in Fig. 5 (top). As already mentioned, the vorticity field is more sensitive to small scale noise, and the optimization of the grid has a positive effect on it. In the process of the optimization, the off-centering is reduced to very small values. Also the range of values of the triangular (and pentagon/hexagon) areas for a given grid level is reduced with the optimization, increasing the homogeneity of the grid in this sense. But the range of values of the lengths of the primal (and dual) grid edges for a given grid level increases with the optimization process. This latter effect could be the reason why some errors increased when the grid is optimized.
When numerical diffusion is applied, the vorticity errors decrease in both cases. The l 2 error in the non-optimized grid now is only ∼2 times larger in the non-optimized grid than in the optimized one. The l ∞ error is still ∼3 times larger in the non-optimized grid. The numerical diffusion eliminates small scale noise to a large extent, but some large errors along the edges of the original icosahedron and the Equator are not significantly reduced and consequently the effect of the numerical diffusion is smaller in the l ∞ error.
We can conclude that the optimization has in general a positive effect.

Sensitivity to the wind reconstruction
Some experiments have been run to test the sensitivity of ICOSWM to the way the wind field is reconstructed from the normal components of the wind to the center of the triangular cells. For the RBF reconstruction, some parameters can be chosen. The radial basis functions used (the kernel), the stencil and a scale factor. For a detailed description of these parameters see Ruppert (2007).
Experiments with the following wind reconstructions have been run: In all the experiments the characteristic damping time is 2 h and the Asselin filter parameter is set to 0.1. Figure 6 shows the convergence results for wind (bottom) and vorticity (top). The convergence for l 2 and l ∞ is shown with solid and dash-dotted lines respectively. As in the other convergence plots, each experiment is identified with a colour and the experiment results are plotted in the order shown in the figure.
Most of the experiments give similar results, meaning that a variety of RBF options can be chosen without changing the ICOSWM performance. The experiments that lead to larger errors help to determine the range of values of the RBF options that are optimal for the reconstruction. The g9p-0.2 experiment gave for all variables significantly larger l 2 and l ∞ errors. The RBF parameters in this experiment are considered inadequate because of the larger errors. The imq9p-0.2 experiment also results in larger wind and vorticity l 2 errors and considerably larger vorticity l ∞ errors. This result can be explained considering that the scale factor defines a kind of influence radius of the RBF. A small scale factor means that the point where the wind is reconstructed can be outside the influence radius of the more remote stencil points, reducing the accuracy order of the reconstruction. The fact that the convergence lines for experiments with Gaussian kernel and different scale factors are parallel, means that it is a reasonable approach the use of a constant scale factor for the different resolutions.
It is remarkable that the experiments with bigger stencils for the reconstruction (RBF 15-and 9-point stencil) do not have better results than the experiments with RBF 3-point stencil and Raviart-Thomas reconstruction. In fact the experiment with a 3-point stencil RBF reconstruction yields slightly better results. In this test case, the wind field is smooth, so a bigger stencil for the reconstruction does not improve the results.
We can conclude that in the case of the RBF reconstruction, the results do not depend significantly on the kernel and the stencil chosen. There is a range of scale factors that give similar good results. For the Gaussian kernel and 9-point stencil values from 0.5 to 1 seem to be adequate, and 0.5 is a good selection for the inverse multiquadratic kernel with 9-point stencil.

Sensitivity to the Asselin filter
Some experiments have been run for different Asselin filter parameters. All of these experiments use the semi-implicit three-time-level scheme, characteristic damping time of 2 h and g9p-0.5 RBF reconstruction. Asselin filter parameters 0.2, 0.1, 0.08, 0.05, and 0.03 have been considered.
The model simulations become numerically unstable for Asselin parameter 0.05 and 0.03 at grid levels 5 and 6. This means that we need an Asselin filter parameter bigger than 0.05.
The normalized errors l 2 and l ∞ for the height, wind and vorticity fields after 10 days are not exactly the same for all the experiments, but the differences are too small to be seen in a convergence plot. Thus a comparison of the convergence on accuracy for the different experiments is not shown.
We can conclude that there is no important effect of the Asselin filter, but it must be larger than 0.05 for numerical stability in the model version documented here.
Some experiments have been tried looking for a combination of parameters that make the model stable for smaller Asselin parameters. It has turned out that the problem is related to the reconstruction of the wind field at the mass points.
There is already available a RBF reconstruction of the tangential velocity at the velocity points in a new version of the model that makes the model stable for very small (or even zero) Asselin filter parameters. It must be stressed that the present work shows the results at a stage (end of 2008) of a project in development. When the tangential velocity is reconstructed at the velocity points, the kinetic energy is calculated at the velocity points and then interpolated to the mass points. It means that a different spatial discretization of the rotation and kinetic energy gradient terms is used. Hollingsworth et al. (1992) show that a computational instability can happen because of a non-cancellation of certain terms in the linearized form of the momentum equation. With the new discretization (RBF reconstruction at the velocity points) the cancellation happens, the instability is not present and no Asselin filter is needed.

Williamson's test case 5
In test case 5 of Williamson et al. (1992) the initial state consists of a zonal flow impinging on an isolated mountain of conical shape. The surface or mountain height h s is given by where h s 0 =2000 m, R=π/9, and r 2 = min R 2 ,(λ − λ c ) 2 + (θ − θ c ) 2 . The center of the mountain is located at λ c =3π/2 radians, θ c =π/6 radians.
The imbalance in the initial state and the presence of the mountain lead to the development of a Rossby gravity wave which propagates all around the globe. This test is relevant to understand the response of the numerical solution to orographic forcing and it has been a common benchmark since the development of the first spectral models.
No analytical solution is available for this test case and a reference model is used to evaluate the errors of ICOSWM. As reference the NCAR STSWM has been used. The spectral resolution for the reference model is T426, the time step is 90 s, the diffusion coefficient is 4.97×10 11 m 4 s −1 and no Asselin filter is applied. The reference solution is available at http://icon.enes.org/swm/stswm/node5.html.
For a variety of model parameters 15-day runs have been done. The spectral reference solution is interpolated by bicubic interpolation from the corresponding Gaussian grid (that has a resolution of 31.25 km at the Equator) to the ICOSWM grids at different grid levels. The difference between the ICOSWM output fields and the interpolated reference solution is used to calculate the l 2 and l ∞ normalized errors.

Sensitivity to the characteristic damping time
Experiments with different characteristic damping times (28 h and 2 h) and without explicit diffusion have been run. All of them use Asselin filter 0.1 and g9p-0.5 RBF reconstruction. Figure 7 shows the height (bottom) and vorticity (top) fields after 15 days for the case of a characteristic damping time of 2 h and grid level 6. In the height plot the mountain height is represented by black contour lines at intervals of 400 m. Superimposed is the NCAR STSWM reference solution as contour black lines. There is a good agreement between the ICOSWM and NCAR STSWM solutions. The general pattern of the vorticity field is very similar in both models.
Except for grid level 2, no significant influence of the characteristic damping time is observed in the height and wind errors in the three experiments. In the case of the vorticity field the l 2 errors are slightly smaller with larger diffusion coefficients for all the grid levels (see Fig. 8). Table 5 shows the numerical values of the normalized errors after 15 days for the different variables from grid level 2 to 6 for the experiment with 2 h characteristic damping time.
A second order convergence for the height field is only achieved for the coarser resolutions. The l 2 errors for the wind and vorticity fields show a convergence rate slightly larger than first order. The l ∞ errors lose convergence for the higher resolutions. Figure 15 in Tomita et al. (2001) shows the temporal evolution of the height and wind error norms for test case 5. The    Table 5 for grid level 6. Grid level 7 in Tomita et al. (2001) and grid level 6 in ICOSWM are equivalent in the sense that they correspond to the same number of bisections performed on the original icosahedron, although the grids differ not only due to the optimization process but on the fact that the grid in Tomita et al. (2001) is an Arakawa-A type grid with all the variables defined at the vertices of the triangular cells. Figure 9 shows the difference with respect to the reference STSWM solution for the vorticity field for grid level 6 (bottom) and grid level 5 (top) with characteristic damping time 2 h, g9p-0.5 RBF reconstruction and Asselin parameter 0.1. The same colour table is used in both cases. Although Fig. 7 (top) shows a good agreement between the vorticity fields in ICOSWM and the NCAR STSWM reference solutions, the difference of both fields show discrepancies at some points. The difference map does not show a wave number-5 pattern, as in test case 2.
The number of grid points with large errors is considerably reduced moving from grid level 5 to grid level 6, reducing the l 2 normalized error. But the largest error is only slightly reduced, and the convergence rate of the l ∞ error is very small. The orography in this test case is not smooth and the representation of this mountain in the spectral model might involve Gibbs phenomena at the edge of the mountain that are not present in the ICOSWM model. Some discrepancies in the lee side of the mountain between the two models seem to be propagated and amplified with time.
One of the purposes in test case 5 is to investigate the global conservation properties. We check the conservation of total energy TE and Potential enstrophy PENS defined as where E p0 denotes the potential energy in the initial state. The definition of total energy in Eq. (28) is the same as in Tomita et al. (2001) and Stuhne and Peltier (1999) but differs from that of Williamson et al. (1992). In Eq. (28) the initial potential energy is subtracted because the interest is to test the variation of the total energy with respect to the available energy. The normalized deviations from the initial values for TE and PENS are then defined as where TE 0 and PENS 0 are the initial values. Figure 10 shows the temporal evolution of TE and PENS for the cases without numerical diffusion and with characteristic damping time of 2 h, for grid level 6. The maximum values of TE and PENS are of the same order of magnitude than the corresponding values in Fig. 17 in Tomita et al. (2001) for grid level 7.

Sensitivity to the wind reconstruction
As in test case 2, some experiments with different wind reconstructions have been run to evaluate the impact of the reconstruction technique in the model results. For simplicity, and considering the previous results for test case 2, only experiments with Raviart-Thomas reconstruction and RBF with Gaussian kernel and different stencils and scale factors have been performed. Figure 11 shows the l 2 (solid lines) and l ∞ (dash-dotted lines) normalized errors after 15 days for the height (bottom) and wind (top) fields for different grid levels. The l 2 and l ∞ errors with Raviart-Thomas and RBF 3-point stencil reconstruction are very similar, the convergence line for both experiments would be indistinguishable and the Raviart-Thomas experiment results are not presented in Fig. 11.
The two experiments with RBF-9 points stencil reconstruction and different scale factors also show very similar results. The 15-points stencil generally produces larger errors than a 9-points stencil. The difference in the errors is reduced with increasing grid level and at grid level 6 all the experiments have similar errors. In the case of the height field (and the vorticity field, not shown here), the errors using the 15-or 3-point stencil are very similar. In the case of the wind field, the 15-point stencil reduces the wind errors for the higher grid levels compared to the 3-point stencil experiments. From these results we can conclude that a 15-point stencil for the reconstruction does not improve the results and is not recommended as it is computationally more expensive.

Sensitivity to the Asselin filter
Experiments with different Asselin filter parameters (0.2, 0.1, 0.08, and 0.05) have been run, with semi-implicit threetime-level scheme, g9p-0.5 RBF reconstruction and characteristic damping time of 2 h.
The convergence plots for l 2 (solid lines) and l ∞ (dashdotted lines) after 15 days for the height and wind fields are shown in Fig. 12 bottom and top, respectively. The normalized errors for the higher resolutions are significantly reduced with decreasing Asselin filter parameter, especially for the height l 2 and l ∞ normalized errors and the wind l 2 normalized errors. In the reference model the Asselin filter parameter is set to zero and the ICOSWM solution is closer to the reference model when a small Asselin filter parameter is used.
The experiment with Asselin parameter 0.05 can not be run for grid level 2. Smith and Dritschel (2006) report that they found a limit for the minimum value of the Asselin filter parameter that could be used in their model. This value is related to the mean short-scale gravity wave speed and depends on the time step. This minimum Asselin parameter increases with increasing time step. For grid level 2, the time step is 1440 s and 0.05 is slightly below the minimum Asselin parameter reported by Smith and Dritschel (2006).
No sensitivity to the Asselin filter parameter is observed in the case of the vorticity field normalized errors.
Following these results, 0.05 would be the best choice for the Asselin filter parameter.

Williamson's test case 6
In test case 6 of Williamson et al. (1992) the initial state consists of a Rossby-Haurwitz wave of wavenumber-4. This type of wave is an analytic solution for the barotropic vorticity equation and has also been widely used to test shallow water models, since the analysis in Hoskins (1973) supported the view that wavenumber-4 is stable also as a solution of the shallow water equations. However, some recent work presented in Thuburn and Li (2000) has shown that the Rossby-Haurwitz wave of test case 6 is actually unstable as a solution of the shallow water equations, since small random perturbations in the initial state or small numerical errors result in long term disruption of the wavenumber-4 pattern. This was shown to be the case for a wide range of numerical models, including spectral transform models. Therefore, the usefulness of the Rossby-Haurwitz wave of wavenumber-4 as a benchmark for the solution of the shallow water initial value problem is limited to time ranges shorter than those sometimes considered in the literature. We choose a run time of 10 days for the study of the convergence of the solutions, although the ICOSWM solution is still stable and the wavenumber-4 is well kept after 14 days. Figure 13 shows the height field at day 14 for an ICOSWM solution with a characteristic damping time of 2 h, Asselin filter parameter 0.1, and g9p-0.5 RBF and grid level 6. The NCAR STSWM reference solution (T511) is superimposed in the figure as black contour lines for comparison. The wavenumber 4 pattern is well kept after 14 days. No phase delay is observed in the solution. In Bonaventura and Ringler (2005) a slight phase delay was observed in the solution produced by the Icosahedral C-staggered grid model. The fact that in Bonaventura and Ringler (2005) a simple twotime-level semi-implicit time discretization is used, while ICOSWM uses a three-time-level semi-implicit time discretization, can be the reason why the phase delay is not observed in ICOSWM. It must also be noted that the NCAR STSWM uses a three time-level scheme.
The NCAR STSWM is used as a reference to evaluate the normalized errors of the ICOSWM. The spectral model resolution is T511, the time step is 90 s, the diffusion coeffi- cient is 3.4×10 12 m 4 s −1 and no Asselin filter is applied. The reference solution is available at http://icon.enes.org/swm/ stswm/node5.html. The corresponding Gaussian grid has a resolution of 26 km at the Equator. The explicit diffusion coefficient used in test case 6 to produce the STSWM reference solution is larger than in test case 5 to reduce the small scale noise in the vorticity field.
Convergence tests after 10-day runs with a variety of model parameters are shown.

Sensitivity to the characteristic damping time
Experiments with different characteristic damping times (28 h, 2 h, and 0.5 h) and without explicit diffusion have been run. All of them use Asselin filter 0.1 and g9p-0.5 RBF reconstruction.
No important variability is observed for the l 2 and l ∞ errors of the wind field. In the case of the vorticity, the l 2 and l ∞ errors are reduced with larger diffusion coefficients (smaller characteristic damping times), see Fig. 14. In the case of the height field, the l 2 and l ∞ errors increase slightly with increasing diffusion but the (negative) effect is smaller than the (positive) effect of the increasing diffusion on the vorticity error norms. Table 6 summarizes the numerical values of the normalized errors after 10 days for the different variables from grid level 2 to 6 for the experiment with 2 h characteristic damping time.
The l 2 normalized height and wind errors show approximately second order convergence. But the l ∞ height and wind errors lose convergence for the higher grid levels. For the vorticity field, both the l 2 and l ∞ errors lose convergence for the higher grid levels. Figure 15 shows the differences of ICOSWM with respect to the STSWM reference for the vorticity field after 10 days and grid level 6, with characteristic damping time Fig. 14. Vorticity field convergence test for case 6 after 10 days. Tests for different characteristic damping times. Solid lines for l 2 errors and dash-dotted lines for l ∞ errors. Table 6. Normalized l 2 and l ∞ errors after 10 days. Case 6. Characteristic damping time 2 h, Asselin filter parameter 0.1, g9p-0.5 RBF reconstruction. 2 h, g9p-0.5 RBF reconstruction and Asselin parameter 0.1. The difference map shows a wavenumber-4 pattern and it is not related to the grid. A characteristic damping time of 0.5 h is a good choice because of the smaller vorticity errors. For test case 2 a characteristic damping time of 1 h or smaller turns out to be a bad choice because it leads to a degradation of the results for all the variables. Thuburn and Li (2000) concludes that in test case 6 a generation of small scales and a cascade towards unresolved scales takes place. Therefore numerical models need a dissipation mechanism to avoid noisy structures in the solution, particularly in the potential vorticity field. Increasing the numerical diffusion in ICOSWM to a characteristic damping time of 0.5 h leads to a less noisy vorticity field and to smaller vorticity errors. Figure 16 shows the evolution of the relative changes of total mass for the runs with no diffusion and with characteristic damping time of 2 h for grid level 6. Figure 17 shows the temporal evolution of T E (top) and PENS (bottom). Although the model is not formally preserving potential enstrophy and total energy, these quantities are nearly conserved. It can be observed that the conservation properties of the model are not significantly affected by the amount of diffusion applied.

Sensitivity to the wind reconstruction
The different reconstructions for the wind field as in test case 5 have been tried with the aim of testing the influence of the reconstruction on the model results. In all the cases the Asselin filter parameter is set to 0.1 and the characteristic damping time to 2 h. Figure 18 shows the height (bottom) and wind (top) l 2 (solid lines) and l ∞ (dash-dotted lines) normalized errors for different grid levels for the different experiments.
For the height and vorticity field (not shown here), starting from grid level 3, the l 2 errors increase if more points are used for the reconstruction. The differences decrease with increasing grid level, and at grid level 6 all the experiments have similar l 2 height and vorticity errors. The l ∞ height and vorticity errors for the experiments with 3-point stencil are smaller than the errors for the experiments with 9-and 15point stencil at grid level 6. In the case of the wind field, the l 2 errors at grid level 6 for the experiments with a 3-point stencil are the largest. The wind is a diagnostic variable, to calculate the normalized wind errors, the wind has to be reconstructed from the normal velocity components and a bigger stencil yields smaller errors. Again we can conclude that the used of a 15-point stencil is not recommended because it is more expensive and does not improve the model results. It is interesting to see that the experiments with a 3-point stencil have better results than the ones with a 9-point stencil, with the exception of the wind l 2 error at grid level 6.

Sensitivity to the Asselin filter
Experiments with different Asselin filter parameters (0.2, 0.1, and 0.08) have been run, with semi-implicit three-timelevel scheme, g9p-0.5 RBF reconstruction and characteristic damping time of 2 h. The convergence plots for l 2 (solid lines) and l ∞ (dashdotted lines) after 10 days for the height field are in Fig. 19.
The height l 2 and the wind (not shown here) l 2 normalized errors for the higher resolutions are reduced with decreasing Asselin filter parameter, the effect being larger in the case of the height field. On the contrary, in the case of the vorticity field, the errors decrease with increasing Asselin filter parameter, but the difference is very small and it is hardly seen in a convergence comparison plot.
The model solutions get numerically unstable for Asselin parameters of 0.05 and smaller for higher grid levels. Thus an Asselin filter parameter larger than 0.05 is needed for stability reasons.
As mentioned in Sect. 3.1.3, a RBF reconstruction of the wind at the velocity points has been developed on a new version of the code which allows for running the model without explicit diffusion and no Asselin filter. Using the same model parameters in both versions of the code leads to very similar normalized errors. When using the new code with the wind reconstruction at the velocity points together to a very small Asselin filter parameter (0.01), leads to smaller l 2 height error norms, particularly for the higher resolutions, increasing the convergence rate of the l 2 height error norms.

Comparison to GMESWM
To evaluate the results of the ICOSWM model, a comparison with GMESWM, the shallow water version of GME (operational global model of Deutscher Wetterdienst, Majewski et al., 2002) has been considered. The GMESWM model uses a non-staggered icosahedralhexagonal grid. The prognostic variables are the height and zonal and meridional velocity components at the centers of the hexagons/pentagons. The number of mass points is equal to the number of vorticity points for a given resolution.
The ICOSWM model uses a C-staggered grid, the prognostic variables are the height at the centers of the primary cells (the centers of the triangles) and the normal velocity components with respect to the cell edge. The vorticity is calculated at the centers of the dual grid (hexagons and pentagons) and the number of mass points is different from the number of vorticity points for a given grid level.
Convergence plots comparing the l 2 height and vorticity errors of ICOSWM and GMESWM for test cases 5 and 6 are shown in Figs. 20 and 21.
In GMESWM the number of equal intervals into which each side of the original icosahedral triangles is divided, n i , is used as a parameter for specifying the resolution of the grid. GMESWM runs for n i =32, 64, 96, and 192 corresponding to a spacing between grid points of about 240, 120, 80, and 40 km, with diffusion coefficients 8×10 15 , 1×10 15 , 4.22×10 14 , and 5×10 13 m 4 s −1 , respectively, have been considered. The Asselin parameter used in the GMESWM runs is 0.03. The ICOSWM results in Fig. 20 correspond to a run with a characteristic damping time of 2 h, g9p-0.5 RBF reconstruction and Asselin filter 0.1. In Fig. 21 the ICOSWM results with characteristic damping times of 2 h and 0.5 h are shown.
For test case 5 (Fig. 20), ICOSWM shows a better height l 2 convergence rate and a similar vorticity l 2 convergence rate as GMESWM.
For test case 6 ( Fig. 21), ICOSWM shows a better height l 2 convergence rate and smaller normalized l 2 height errors. With a characteristic damping time of 2 h, ICOSWM shows a worse vorticity l 2 convergence rate than GMESWM, but with a characteristic damping time of 0.5 h, ICOSWM has similar results as GMESWM. It has been seen that test case 6 needs more numerical diffusion to eliminate small scale noise, which is more noticeable in the vorticity field. The mean diffusion coefficient for grid level 6 with characteristic damping times of 2 h and 0.5 h are 3.1599×10 12 and 1.2640×10 13 m 4 s −1 . The latter is the same order of magnitude than the one used in GMESWM.  We can conclude that in general the ICOSWM results are better than those of GMESWM. It could be due to the Cgrid formulation and the Heikes-Randall grid optimization. To test the effect of the grid optimization, ICOSWM has also been run for test case 6 without optimization with g9p-0.5 RBF reconstruction, Asselin filter 0.1 and characteristic damping times of 2 h and 0.5 h. The normalized errors for the different variables obtained (not shown here) are very similar to the ones obtained with the optimized grid. The errors in test case 6 are larger than in test case 2, hiding the benefits of the optimization observed in test case 2. Without optimization, ICOSWM has still better results than GMESWM.

Kinetic energy spectra
In this section the kinetic energy spectra for some experiments are shown. They are also compared to the kinetic energy spectra of the variant of the NCAR STSWM used as a reference to evaluate the ICOSWM output fields. The kinetic energy spectra are calculated using Eq. (32) (Eq. 3.2 in Jakob- Chien et al., 1995) where ζ m n and δ m n are the spectral coefficients of the divergence and the vorticity, n is the total wavenumber, m is the longitudinal wavenumber and a is the radius of the earth.
To calculate the spectra, the divergence and vorticity fields are first interpolated to a Gaussian grid (T426 or T511). Then the spectral divergence and vorticity coefficients are calculated. Figure 22 shows the kinetic energy spectra at day 15 for the NCAR STSWM model and ICOSWM with different characteristic damping times (three-time-level semi-implicit scheme, Asselin parameter 0.1, g9p-0.5 RBF reconstruction) for grid level 6.

Test case 5 kinetic energy spectra
There is a good agreement between the kinetic energy spectra of both models. With decreasing characteristic damping time (increasing diffusion coefficients), the kinetic energy of the highest wavenumbers decreases in ICOSWM.

Test case 6 kinetic energy spectra
Divergence at initial time is zero in test case 6, as is the initial tendency of the divergence field. This initial state is a solution of the barotropic equations and in a barotropic model the solution would remain divergence free. In a shallow water model, some divergence is produced. The energy of the even wavenumbers for n>20 in the NCAR STSWM model is much smaller than in the ICOSWM, indicating that less divergence is produced by the NCAR STSWM model. Apart from this, there is a good agreement between the kinetic energy spectra of both models.
The kinetic energy spectra show a n −3 dependence for the odd wave numbers till wavenumbers of about 160. The kinetic energy spectrum of the NCAR STSWM evolves with time towards a n −3 dependence for wavenumbers 20<n<160 (not shown here), reaching an equilibrium state between day-10 and day-15. It is not clear if the n −3 dependence is developed only for this range of wavenumbers or if more time is needed to develop it for larger wavenumbers or if the n>160 wavenumbers are just affected by the numerical diffusion. In the latter case, we could not use the NCAR STSWM kinetic energy spectrum as the true kinetic energy spectrum and use it as a reference to define the effective resolution of ICOSWM as it is proposed in Skamarock (2004).
Decreasing the characteristic damping time (increasing diffusion coefficients) in ICOSWM, the kinetic energy of the even total wavenumbers decreases for n>20, but for a characteristic damping time of 2 h, it is still much larger than for the NCAR STSWM. Decreasing the characteristic damping time also reduces the energy of the higher (odd and even) total wavenumbers.
To estimate the impact of the interpolation from the ICOSWM grid to the Gaussian grid on the calculated ICOSWM kinetic energy spectra, the spectrum of the initial state of test case 6 is considered. Figure 24 shows the spectra for the STSWM and ICOSWM models at initial time. The divergence is zero and the vorticity is a linear combination of the spherical harmonics Y m=0 n=1 and Y m=4 n=5 (Williamson et al., 1992). The spectra of both models show peaks at the expected total wavenumbers n=1 and n=5. The energy for the other total wavenumbers should be zero. Values for STSWM are typically ∼10 −29 , corresponding to machine precision, and ∼10 −8 for ICOSWM. We can conclude that the accuracy of the ICOSWM spectra is limited to ∼10 −8 due to the interpolation from the icosahedral grid to the Gaussian grid.

Conclusions and outlook
The Icosahedral Shallow Water Model (ICOSWM) is described and results for tests cases 2, 5, and 6 of Williamson et al. (1992) are presented for a variety of model parameters. For test cases 5 and 6 the NCAR STSWM is used as reference.
ICOSWM uses a grid optimized with the method suggested by Heikes and Randall (1995b). The benefits of the optimization can be observed in test case 2, where the error models are small enough. In test case 6, with larger errors, the benefits of the optimization are hidden.
ICOSWM simulations for test cases 5 and 6 have been compared to the Shallow Water Version of the current operational model at DWD (GMESWM), ICOSWM shows better results than GMESWM, probably because of the C-grid formulation. The ICOSWM normalized height and wind errors for test case 5 (day 15) and test case 6 (day 10) are very similar to those shown in Tomita et al. (2001), although the figures in Tomita et al. (2001) do not allow to do a more quantitative comparison.
The kinetic energy spectra of ICOSWM has been calculated and a good agreement with the kinetic energy spectra of the NCAR STSWM model is observed.
An instability is observed in ICOSWM that is overcome with a relative large value of the Asselin Filter parameter (0.05-0.08). It seems that it is related to the reconstruction of the wind at the center of the triangular cells. The current work presents the results of the ICOSWM at the end of the year 2008, but the model is being further developed. A RBF reconstruction of the wind at the velocity points has been developed on a new version of the code which allows for running the model without explicit diffusion and no Asselin filter. Hollingsworth et al. (1992) show that an instability can happen because of a non-cancellation of certain terms in the linearized form of the momentum equation. With the new discretization (RBF reconstruction at the velocity points), the cancellation happens, the instability is not present and no Asselin filter is needed. Using the same model parameters in both versions of the code leads to very similar normalized errors. When using the new code with the wind reconstruction at the velocity points together to a very small Asselin filter parameter (0.01), leads to smaller l 2 height error norms in test case 6, particularly for the higher resolutions, increasing the convergence rate of the l 2 height error norms.
In the framework of the ICON project a hydrostatic dynamical core has been developed (Wan, 2009) and a local grid refinement option and a non-hydrostatic dynamical core are under development. Among other things, these new models include the above-mentioned new implementation of