Non-orthogonal version of the arbitrary polygonal C-grid and a new diamond grid

. Quasi-uniform grids of the sphere have become popular recently since they avoid parallel scaling bottle-necks associated with the poles of latitude–longitude grids. However quasi-uniform grids of the sphere are often non-orthogonal. A version of the C-grid for arbitrary non-orthogonal grids is presented which gives some of the mimetic properties of the orthogonal C-grid. Exact energy conservation is sacriﬁced for improved accuracy and the resulting scheme numerically conserves energy and potential enstrophy well. The non-orthogonal nature means that the scheme can be used on a cubed sphere. The advantage of the cubed sphere is that it does not admit the computational modes of the hexagonal or triangular C-grids. On various shallow-water test cases, the non-orthogonal scheme on a cubed sphere has accuracy less than or equal to the orthogonal scheme on an orthogonal hexagonal icosahedron. A new diamond grid is presented consisting of quasi-uniform quadrilaterals which is more nearly orthogonal than the equal-angle cubed sphere but with otherwise similar properties. It performs better than the cubed sphere in every way and should be used instead in codes which allow a ﬂexible grid structure.


Introduction
Quasi-uniform grids of the sphere have become popular recently since they avoid parallel scaling bottlenecks associated with the poles of latitude-longitude grids.The predominant groups of quasi-uniform grid are hexagonal icosahedral, triangular icosahedral and cubed-sphere (Weller et al., 2009).There is also an octagonal grid (Rančić et al., 2008) that has not been used much but has similar properties to the cubed sphere and there are reduced or skipped lat-lon grids which are not much used except in conjunction with spectral transform models (Hortal and Simmons, 1991;White, 2003).The details of the grid are critically important for low-order finite volume methods that rely on super-convergence for accuracy (second order accuracy only for a sufficiently smooth grid).
The hexagonal C-grid has become popular since Thuburn (2008), Thuburn et al. (2009) and Ringler et al. (2010) worked out how to calculate the Coriolis term so as to get steady geostrophic modes.This has been named TRiSK.TRiSK has mostly been used on Voronoi tesselations of the sphere (e.g.Ringler et al., 2008) which are orthogonal (the primal and dual edges cross at right angles) and each shape has more than (or occasionally equal to) four sides.
C-grids based on primal cells with more than four sides in 2-D will have more than twice as many velocity degrees of freedom (dofs) as mass dofs and will therefore suffer from spurious computational modes (Staniforth and Thuburn, 2012).The hexagonal C-grid suffers from a branch of spurious Rossby modes (Thuburn, 2008) which do not interact correctly with the mass.The triangular C-grid does not have enough velocity dofs and so suffers from spurious divergent modes (Danilov, 2010;Gassmann, 2011).The spurious modes on triangles can be controlled by strong diffusion (Gassmann, 2011) or strong hyper-diffusion (Wan et al., 2013).The spurious modes on hexagons can be controlled using upwinded advection of potential vorticity (e.g.Weller, 2012) which does not destroy energy.However a more efficient discretisation would have the correct ratio of dofs and would not need to control spurious behaviour in the excess dofs.The correct ratio of dofs can be achieved by using grids of quadrilaterals, such as the cubed-sphere grid.However grids of the sphere using quadrilaterals are either non-orthogonal (e.g. the equal-angle cubed sphere, Fournier et al., 2004), have large variations of cell size (e.g. the conformal cubed sphere, Rančić et al., 1996) or are locally inhomogeneous (such as kite grids, Weller et al., 2012).This provides motivation for more nearly orthogonal quadrilateral grids of the sphere and a C-grid discretisation with the required mimetic properties on nonorthogonal grids.Thuburn and Cotter (2012) describe some desirable mimetic properties of atmospheric models (mimicking the properties of the continuous equations).Their mimetic properties are 1-6.Property 7 is clearly also desirable: 7. Second-order accuracy (or higher).
The TRiSK scheme (Thuburn et al., 2009;Ringler et al., 2010) gives properties (1-6) on orthogonal polygonal grids but it will be demonstrated that the discretisation of the perpendicular (perp) operator (for calculating dual grid fluxes from primal grid fluxes) is inconsistent (i.e.zeroth order accurate) even on the smoothest hexagonal icosahedral grids of the sphere.Extending TRiSK to non-orthogonal grids may ameliorate the lack of convergence of TRiSK since points other than the Voronoi generating points can be used as the cell centre in order to optimise aspects of the grid to improve accuracy of the perp operator.Thuburn and Cotter (2012) set out the mathematical constraints for mimetic C-grid discretisations on non-orthogonal grids but did not give an example of such a scheme.Subsequently, Thuburn et al. (2013) proposed a scheme suitable for grids whose duals consist of only triangles and quadrilaterals and whose duals are centroidal (primal vertices are at the centroids of the dual cells).However the results on cubed-sphere grids were much less accurate than those using similar resolution hexagonal-icosahedra.A variety of mixed finite-element schemes for grids of triangles or quadrilaterals have been proposed which give the above properties and second-order accuracy by constructing and inverting global mass matrices at every time-step (Cotter and Shipton, 2012;Cotter and Thuburn, 2013).Hollingsworth et al. (1983) described an instability that can grow when solving the primitive equations in 3D using the vector-invariant form of the momentum equation, conserving energy and enstrophy but not momentum.Gassmann (2013) found that this mode could grow when solving the fully compressible Euler equations on a hexagonal-icosahedral grid of the sphere using a C-grid discretisation and described how it can be controlled.It is possible that this mode grows more quickly when it interacts with the computational modes of the hexagonal C-grid but this is not proved and has not been demonstrated.If the discretisations described on various grids of quadrilaterals were extended to 3D, the behaviour of the Hollingsworth instability could be compared on hexagonal and quadrilateral grids.Gassmann (2013) found that this mode is triggered at the pentagons of the icosahedral grids.The cube corners of the cubed sphere grid have larger distortions that the pentagons of the icosahedral grid.Therefore it seems likely that this mode would also be triggered on a cubed-sphere grid.
A new diamond grid of quadrilaterals is introduced in Sect. 2 which is more nearly orthogonal than the equal-angle cubed sphere and nearly as uniform.The properties of the diamond grid are compared with those of the cubed sphere and orthogonal and non-orthogonal versions of the hexagonal icosahedron.In Sect.3, a more accurate non-orthogonal model is proposed that forgoes energy conservation for better accuracy than the scheme of Thuburn et al. (2013) and which can be used on grids with non-centroidal duals.The accuracy of the perp operator and the non-orthogonal correction is explored in Sect. 5 and the results of shallow-water test cases are presented in Sect.6.
2 Quasi-uniform grids of the sphere Seven types of grid are considered, some of which are displayed in Fig. 1.The grids are: 1.The Heikes and Randall (1995) optimised version of the orthogonal hexagonal-icosahedron (referred to as the HR grid).Neither this grid nor its dual are centroidal.

2.
A non-orthogonal version of the hexagonal icosahedron with dual vertices moved from the Voronoi generating points to the centroids of the polygons, making the primal grid centroidal.

3.
The centroidal equal angle cubed sphere (with dual vertices at the primal cell centroids).

4.
A diamondised version of the cubed sphere.The diamond grid is constructed by replacing each edge of the cubed sphere with a primal cell whose vertices consist of the two vertices of the original edge and the cell centres either side of the edge (bottom right of Fig. 1).The dual vertices are then placed at the primal cell centres to make the primal grid centroidal.Versions of 2-4 above but with centroidal dual grids rather than centroidal primal grids.So once the duals are defined, the primal vertices are moved to be at the centroids of the dual grids.These grids are used since the non-orthogonal algorithm defined by Thuburn et al. (2013) is only consistent if the dual grid is centroidal.

Discussion
The diamond grid is topologically different from the cubed sphere and different from the dual of the cubed sphere although it still suffers from the problem of having 3 quadrilaterals meet at one vertex at 8 locations in the grid.The panels of the cubed sphere and diamond grids are shown in Fig. 2. The diamond grid for an equal angle cubed sphere grid at a cube corner is also shown in Fig. 2.This shows that, near the cube corners, the diamond grid cells become rectangles with aspect ratio √ 3 and, along the edges, the cells are kite shaped.In the limit of high resolution, it would be possible to construct the diamond grid to be orthogonal but then both the primal and dual grids would be highly non-centroidal, leading to large truncation errors.Instead diamond grids with either centroidal primal or dual grids are used.
The skewness and non-orthogonality of coarse versions of these grids are shown in Fig. 1.The non-orthogonality is the difference between the angle between the primal and dual edges (in degrees) and 90 • .The non-orthogonality is shown from black (orthogonal) to blue (nonorthogonal) on the primal grid edges.The skewness of edge e, s e , measures the departure from the edge centre of the primal-dual edge cross-over point: where x v and x w are the primal vertices at either end of edge e.The skewness of the primal and dual meshes is different but the skewness of the primal mesh is shown from red (no skewness) to yellow on the dual edges.In Fig. 1, the diamond grid is more nearly orthogonal and less skew than the cubed sphere, a result that holds at all resolutions considered.The Heikes and Randall (1995) (HR) grid (top left of Fig. 1) is orthogonal and optimised to minimise skewness.This grid has recently been revisited by Heikes et al. (2013).The HR optimisation minimises the error of discretising a Laplacian but the value at the dual vertex is not a second-order approximation of the primal cell average because the dual vertex is not at 6 the centroid.In moving the dual vertex to the centroid of the primal cell (top right of Fig. 1), the grid becomes centroidal but non-orthogonal and also the skewness is increased.An alternative is the centroidal Voronoi grid (Ringler et al., 2008) which is orthogonal but more skew than the HR grid.Using non-orthogonal grids opens up many more options for optimising a combination of the orthogonality, skewness, uniformity and centroidality of both the primal and dual grids.
However this has not been done.Some of the properties of the grids at different resolutions are shown in Tables 1 and 2. The non-centroidality of the primal is defined as the distance between the dual vertex and the cell centroid of a primal cell divided by the square root of the primal call area: By making the primal centroidal, the dual may become less centroidal.The centroidal hexagonal grid (Table 1) has non-orthogonality of less than 1 • and skewness similar to the orthogonal version whereas the hexagonal grid with a centroidal dual (Table 2) has much larger nonorthogonality -up to 13 • .The centroidal cubed-sphere (Table 1) has non-orthogonality increasing with resolution up to 30 • for the resolutions considered and maximum skewness of 0.25 at the corners.The ratio of maximum to minimum grid spacing reaches 1.78 for the grids presented in Table 1 which is larger than the asymptotic value of 1.3 given by Staniforth and Thuburn (2012).This is because we are measuring the cell centre to cell centre distance rather than grid edge length.In moving the dual vertices to the primal cell centroids, cell centres have become closer together at the cube corners.The cubed sphere with a centroidal dual (Table 2) has similar properties but is slightly less skew at the corners.Both diamond grids are more orthogonal than the cubed sphere in the mean and maximum (less that 9 • ) and the skewness and non-centroidality are also smaller.The diamond grid is slightly less uniform that the cubed sphere ( ∆xmax ∆x min < 2.09 for the diamond grid whereas ∆xmax ∆x min < 1.8 for the cubed sphere) but otherwise does not appear to suffer from any deficiencies relative to the cubed sphere.Again, like the cubed-sphere, the cell centre to cell centre distances vary more than the edge lengths.(The maximum to minimum edge length for the diamond grid should approach √ 3 ≈ 1.7.)The impacts of the different grid structures on the accuracy of the perp operator (for estimating the velocity perpendicular to the normal velocity at each edge) and on the non-orthogonal correction will be seen in Sect. 5 and on the solution of the shallow-water equations in Sect.6.

The non-orthogonal C-grid discretisation
We present a discretisation of the rotating, non-linear shallow-water equations in vectorinvariant form in which the continuity and momentum equations are: where φ is the geopotential (hg, fluid depth times gravity), v is the horizontal velocity, v ⊥ = k×v where k is the local unit vertical vector, ζ = f +ξ is the absolute vorticity, where f = 2k•Ω is the Coriolis parameter associated with rotation Ω and ξ = k • (∇ × v) is the relative vorticity.

Notation
The notation has some minor differences from Thuburn and Cotter (2012).The primal (solid) and dual (dashed) grids from Thuburn and Cotter (2012)  denoted by v or w.These definitions and some of the finite volume approximations are given in Table 3.

Discretised momentum and continuity equations
The prognostic variables of the shallow-water equations of a C-grid are usually cell average geopotential, φ i , and the normal component of the velocity at the cell edges, u e = v • pe .However on the non-orthogonal C-grid, the prognostic velocity variable is V e (Thuburn and Cotter, 2012).We consider split space-time discretisation and so the discretisation of the temporal derivatives is considered separately.The spatially discretised continuity equation for φ i and momentum equation for V e can be written: The discretisation of each of the terms will now be described, including the H operator and the graident along a dual edge, ∇ d .

Perp operator, ⊥
The perp operator, ⊥ , in Eq. ( 6) is the discrete operator described in Thuburn et al. (2009) for mapping from fluxes on the primal grid to fluxes on the dual grid.This operator ensures that the divergence of the mapped fluxes is a convex combination of the divergence of the fluxes on the primal grid.This is why the discrete perp operator acts on the flux, U e , (in the p direction) and maps to the d direction, despite p and d not being at right angles.This operator is inconsistent (zeroth order accurate) even on an orthogonal grid but does not always prevent convergence with resolution of the shallow-water equations (see Sects.5.1 and 6).The perp operator is defined to be: e ∈EC(i),EC(j) where primal cells i and j are the cells either side of edge e p , EC(i) means the edges of cell i and Thuburn et al. (2009) derived the weights w ee : where the vs are the vertices in a walk between edges e and e .If the walk starts in the p ⊥ direction and if n ei = 1 then the sign is positive.

Mapping from primal cell averages to edges
In order to ensure energetic consistency, the mapping of φ from cells to edges (φ i → φ e ) must use the same weighting as will be used for calculating the kinetic energy (Eq.10): where cells i and j are either side of edge e.This is the reverse of linear interpolation and ensures an exact transfer between kinetic and potential energy.However, higher-order upwind interpolations can be used instead, foregoing this form of energy conservation in favour of a smoother geopotential field.Here, we will present results using CLUST (Weller, 2012) with a blending of 50 % between linear and linear-upwind differencing which gives smoother advection of geopotential (see Sect. 3.6).

PV and curl on dual cells
From Ringler et al. (2010), q v is discretised as: A iv φ i is the conservative mapping of φ from the primal to the dual grid and, from Thuburn and Cotter (2012), the curl is discretised as:

Mapping pv from dual cells to edges
The potential vorticity at the edge, q e , is interpreted as the pv at the primal and dual edges.It is interpolated from surrounding q v values from an upwind-biased stencil using CLUST which was developed for mapping pv from vertices to edges of the polygonal C-grid Weller (2012).
The CLUST blending coefficient between linear differencing and linear upwind used is 0.5.It is essential to use CLUST for this work rather than a conventional high-order upwind or monotonic advection scheme such as quadratic-upwind, linear-upwind or a TVD scheme such as van Leer.The conventional schemes have switching between upwind on either side of an edge when the flow is aligned with the edge.When used for interpolating pv on the dual of the C-grid, this leads to errors.The advantage of CLUST is that it blends smoothly with linear when the flow is aligned with the edge so there is no switching.APVM (Ringler et al., 2010, which is equivalent to Lax-Wendroff) could also be used as it does not involve switching but CLUST was found to control the grid-scale pv noise better (Weller, 2012).

Energy conserving Coriolis flux averaging
The averaging between q e (φ e U e ) ⊥ and (q e φ e U e ) ⊥ in Eq. ( 6) is necessary for the Coriolis force to be energetically neutral (Ringler et al., 2010).Without this averaging (if this term is simply represented as q e (φ e U e ) ⊥ ) the pv evolves exactly as if it were advected by the fluxes U ⊥ e by the advection scheme used to map pv from dual cells to edges.Thus high order, monotonic advection of pv can be obtained.

Gradients along dual edges
The gradient of φ along dual edge e, integrated along the edge, between cells i and j is: This gradient is also used for the kinetic energy.

Kinetic energy
The kinetic energy in primal cell i is defined as: where dist (x e , x i ) is defined in Sect.3.13.Ringler et al. (2010) weighted the edge contributions by 1 2 rather than A ie /A e .Then, in order to achieve energetic consistency, they use A i = 1/4 e A e which gives the correct area on Voronoi grids for which the edges bisect the lines between the Voronoi generating points.Weller et al. (2012) suggest the weighting as in Eq. ( 10) for non-Voronoi grids.For energetic consistency, the same weighting must be used for mapping the geopotential from cells to edges (see Sect. 3.4).

Divergence operator
The divergence of a vector v e defined at edges in cell i is given by Gauss's divergence theorem: Thuburn and Cotter (2012) suggest operator H to transform V e into U e so that: The H operator is described further in Sect.3.11 below.

H operator
The H operator transforms from the set of V values to the U values and is therefore referred to as a non-orthogonal correction.

Symmetric H
Thuburn and Cotter (2012) proved that energetic consistency is achieved if H is symmetric and positive definite but they did not suggest the form of H for non-orthogonal grids.Thuburn et al. (2013) suggest an H operator with the desired properties for grids whose dual consists of quadrilaterals or triangles and for which the dual grid is centroidal: where the stencil of edges e consists of the edges of the dual sharing one vertex with e and f e = 4 when dual edges e and e are edges of the same quadrilateral and f e = 6 when e and e are edges of the same triangle.All Voronoi grids and many other grids of polygons have dual grids consisting of only triangles and the dual of the cubed sphere grid consists of quadrilaterals and 8 triangles.So this operator should cover most of the grids anyone would be interested in as long as the primal vertices are moved to the dual cell centroids.The centroidal dual constraint is necessary to ensure first-order accuracy.However, for grids with triangular duals, the offdiagonal terms of H do not vanish as the grids tends towards an orthogonal grid.The accuracy of solutions using this operator will be presented in Sect.6.
where T i = e∈EC(i) p e p e .T i is a 3×3 tensor so not computationally expensive to invert and is only dependent on geometry and so can be pre-computed.It can be shown that this operator will exactly reconstruct a uniform velocity field and for non-uniform velocities, it is a least-squares fit.It is similar to Perot's reconstruction (Perot, 2000) but a comparison of the two methods has not been done.Alternatively, to reconstruct dual cell velocity, v v from V e we can use: where We can then interpolate the v v s from the dual cells to the dual edges using any centred or upwind-biased interpolation.We will show results in Sect.6 using mid-point interpolation.The resulting velocity on edges is referred to as v e , with the prime because this is not the final velocity that will give us U e .We require U e = p e /d e V e for an orthogonal grid ( p = d ⊥ ) so we can correct v e to give exactly V e in the direction d For nearly orthogonal grids this operator is close to diagonal so the requirement for positive definiteness should easily be met for most grids.However this operator is not symmetric (unless of course it is diagonal).In fact the requirement for symmetry is not consistent with the requirement that H ee = 0 for all e = e if p e • d ⊥ e = 1, since if we also have p e • d ⊥ e = 1 then we may have H e e = 0.
The H operator described in this subsection is referred to as the asymmetric/diagonal H operator (or just the asymmetric operator for short) since it is asymmetric for a non-orthogonal grid and diagonal for an orthogonal grid.The relative merits and accuracy of the symmetric and asymmetric/diagonal H will be assessed in Sect.6.

Full velocity field at cell edges
The full velocity field at cell edges is needed in CLUST and for post-processing such as calculating errors and error metrics.Due to the inconsistency of the perp operator, Eq. ( 16) is used to reconstruct the full velocity field.

Spherical areas and distances
All distances (or lengths) are great circle distances so that the distance between points x v and where a is the magnitude of both x v and x w .The areas on the surface of the sphere are calculated to be consistent with the distances to retain the correct mimetic properties.Thus the area of a triangle with points x, y and z on a sphere of radius a is: Areas A ie , A iv , A i and A v are composed of the sums of triangles A ive , as shown in Fig. 3.

Semi-implicit solution technique
The momentum and continuity equations are solved semi-implicitly using Crank-Nicolson for the implicit terms.Two outer iterations are used so that, for the first iteration, the explicit terms are solved with Euler-explicit and the second iteration uses 50-50 weighting of the old and new values so that the explicit terms are second order, with the same weighting as the implicit terms.
For the simplest possible implementation, H is split into diagonal and off diagonal elements: H = H d + H off and only the diagonal terms are included in the implicit part.The off diagonal terms are lagged corrections.Between times-levels n and n + 1 this gives: where values at time level are at time-level n for the first iteration and they are the most up to date value (but not implicit) for the second iteration.Using just two iterations, the explicit scheme is the Heun scheme which is weakly unstable.However the instability is not seen in the simulations undertaken.An additional explicit step would remove the instability if needed.Equations ( 21) and ( 22) are solved simultaneously for φ n+1 by substituting V n+1 from Eq. ( 22) into Eq.( 21) (taking the Schur complement) to form a Helmholtz equation which is solved using a conjugate gradient solver with incomplete Cholesky pre-conditioning (OpenFOAM, cited 2013).
The normal modes (eigen vectors) and corresponding amplification factors (eigen values) of the model for the linearised shallow-water equations More explicit updates are not needed because, after 2, the fields do not change for this small time-step (to within machine precision) so the time-stepping is effectively Crank-Nicolson.Crank-Nicolson time-stepping is neutrally stable and the symmetric H operator is energy conserving and therefore the eigen values using the symmetric H have magnitude 1 for both grids in Fig. 4. Use of the asymmetric H does not formally guarantee energy conservation but, on the grids with centroidal duals, the eigen values still have magnitude 1. However using the centroidal grids, some eigen values have magnitude greater than one, implying instability.(The symmetric H gives the same eigen value magnitudes on the centroidal grids but this is not shown because the symmetric H is inconsistent on grids without centroidal duals.)We can conclude that the asymmetric H should not affect stability but that the use of a non-centroidal dual grid may affect stability.
The zero frequency modes confirm the existence of steady geostrophic modes.The accuracy of the TRiSK perp operator (Eq.7) is assessed by reconstructing the solid body rotation velocity field of test case 2 of Williamson et al. (1992) from the normal components of this velocity field at each edge.The maximum errors in the reconstructed velocity in comparison to the analytic profile for each grid at each resolution are shown in Fig. 5.Total degrees of freedom (dofs, number of cells plus number of edges) is used as a proxy for both resolution and computational cost although it is not directly proportional to either.It can be argued that models with a non-ideal ratio of dofs will have lower effective resolution since some dofs must be slaved to others in order to avoid computational modes.However, for consistency with other studies (e.g.Weller et al., 2012;Thuburn et al., 2013), we will stick with using total dofs as a measure of resolution.
The errors of the TRiSK perp operator on the centroidal grids in Fig. 5 saturate at quite low dof count and are highest on the cubed-sphere.It is surprising how much more accurate the perp operator on the centroidal hexagonal grid is in comparison to the orthogonal version.
On the grids with centroidal duals, there is no convergence with resolution of the TRiSK perp operator.
The perp operation can also be evaluated using a least-squares fit (Eq.16) which is also shown in Fig. 5.The least squares fit is much more accurate than the TRiSK perp and converges with resolution but does not satisfy the important mimetic property that the divergence on the dual is a convex combination of the divergence of the primal.Therefore if the shallow-water equations are solved using the least squares operator, steady geostrophic modes would not be maintained which leads to considerably less accurate solutions (not shown).

Accuracy of the H operator
The accuracy of the symmetric and asymmetric H operators (Sect. 3.11)  to cell centre direction for the solid body rotation velocity field of test case 2 of Williamson et al. (1992).The maximum errors for each of the non-orthogonal grids at each resolution are shown in Fig. 6.The centroidal primal grid is particularly beneficial for the hexagonal grid, presumably because this grid is so close to orthogonal.The diamond grid with a centroidal dual has insufficient convergence with resolution using both H operators but it is particularly bad for the symmetric H.This is contrary to the analysis of Thuburn et al. (2013) that the symmetric H is first-order accurate on centroidal grids.But it should be noted that the grid is changing as resolution increases.In particular, the centroidal dual diamond grid is getting less orthogonal.
On the centroidal dual grids, the asymmetric H is more accurate than the symmetric H.

Results of shallow-water test cases
The time-steps used for each test case, each grid and each resolution are shown in Table 4. Time-steps are chosen to give similar advective and gravity wave Courant numbers for each grid.However the Courant numbers chosen are different for each test-case, as described in the sub-sections describing the test-case results below.
6.1 Solid body rotation (Williamson et al., 1992, test case 2) Height errors and height contours after five days for the Williamson et al. (1992) solid body rotation, test case 2 are shown in Fig. 7 for coarse versions of the grids, each with similar numbers of total dofs.The height errors on the orthogonal, hexagonal HR grid are the lowest but the errors on the centroidal hexagonal and centroidal diamond grids are also low.These are the grids that are the most orthogonal but they do not necessarily have the lowest H errors (Fig. 6) or the lowest perp errors (Fig. 5).All of the cubed sphere grids and the diamond grids with centroidal duals have much higher errors, regardless of the H operator used.The centroidal dual hexagonal grid has high errors along the lines of non-centroidal primal cells.The version of H used for this makes little difference.The disadvantage of the asymmetric H is that it spoils the energy conservation properties of the spatial discretisation.In order to judge the extent of this problem, the normalised energy change for symmetric and asymmetric/diagonal H for the simulations shown in Figs.7 are shown in Fig. 8. (The kinetic energy is defined as in Eq. ( 10) and the energy change is normalised as described by Williamson et al. (1992) by dividing by the initial total energy.) Positive changes are solid and negative changes are dashed.The symmetric H does indeed have marginally better energy conservation.The centroidal hexagonal and centroidal diamond grid have energy conservation very similar to the orthogonal hexagonal grid which formally conserves energy.The good energy conservation of the simulations using the asymmetric H are consistent with the results of the linear stability analysis in Sect. 4. However the energy conservation of the centroidal cubed sphere is less good, consistent with the high truncation errors seen if Fig. 7.
The motivation for the grids of quadrilaterals is to avoid the computational modes of the hexagonal C-grid.The computational Rossby modes manifest themselves as grid-scale enstrophy.This is controlled using upwind advection of pv (CLUST is used here on all grids and for both model versions).The solid-body rotation test on the orthogonal and centroidal hexagonal grids loses total enstrophy (bottom left of Fig. 8), related to the existence of (controlled) computational modes.The cubed-sphere and diamond grids do not have these computational modes and their enstrophy conservation is about an order of magnitude better.Enstrophy conservation on the grids with centroidal duals is better but there are more instances of increasing enstrophy (solid lines) rather than decreasing enstrophy (dashed lines) which is consistent with growing grid-scale noise.
Convergence with degrees of freedom for the solid-body rotation test on all grids of the 2 and ∞ error norms of geopotential and velocity is shown in Fig. 9. Degrees of freedom (dofs, number of cells plus number of edges) is used as an approximate measure of resolution.Time steps are chosen (see Table 4) to maintain an advective Courant number of about 0.14 and a gravity wave Courant number of about 0.7, apart from at the lowest resolutions which need a shorter time-step in order to represent the scale independent Coriolis term accurately.Larger Courant numbers could have been chosen but the largest time-step cannot go much above 3600 s on the coarsest grid in order to maintain accuracy of the Coriolis term.
Convergence with resolution of all error metrics of this shallow-water test case is much better than the convergence with resolution of the perp operator alone.It seems that having the right divergence on the dual is more important for accuracy than convergence of the perp operator.In Fig. 5, the non-orthogonal HR grid has the lowest errors of the perp operator whereas solving the shallow-water equations with the same initial wind field, the orthogonal HR grid has the lowest errors, implying that other aspects of the discretisation are controlling the errors of the shallow-water model for this test case, not the perp operator.
The orthogonal HR grid has the best convergence with resolution.On the centroidal cube and diamond grids, the errors are low but the convergence of even the 2 (φ) error norm slows to less than first order at high resolution.This is solved on the centroidal dual grids with both H operators.Of the grids with centroidal duals only the hexagonal grid has at least first order convergence of ∞ (φ).For all cubed-sphere and diamond grids, the ∞ (φ) errors stop converging with resolution after about 10,000 dofs.
A least squares fit is used to calculate the velocity for the error norms and so the inconsistent perp errors do not appear directly in the error norms.As a consequence, the 2 (u) error norms all converge with second order and the ∞ (u) between first and second.
From this sub-section, we have learnt the following: -The centroidal grids give lower errors but the centroidal dual grids give better convergence with resolution.
-The asymmetric H on the centroidal dual grids gives similar accuracy to the symmetric H and does not degrade accuracy on an orthogonal grid.
-Making the hexagonal grid centroidal increases the error a little in comparison to the orthogonal HR grid.
-The cubed-sphere grids have higher errors at all resolutions considered and lower order of accuracy in comparison to the hexagonal grids.-The centroidal diamond grid has better convergence with resolution and lower errors in comparison to the centroidal cubed-sphere grid.
-The accuracy of the solution of the shallow-water equations for this test case is not directly related to the accuracy of the perp operator.
6.2 Mid-latitude mountain (Williamson et al., 1992, test case 5) The Williamson et al. (1992) flow over a mountain (test case 5) does not have an analytic solution and so numerical solutions are compared with results of a version of the NCAR spectral transform shallow-water model (STSWM Hack and Jakob, 1992) revised by Pilar Rípodas from Deutscher Wetterdienst (Rípodas et al., 2009) and run at T426 resolution and using a time-step of 90 s and a hyper-diffusion coefficient of 4.96 × 10 11 m 4 s −1 .The spectral model results are interpolated from the native spectral model grid (640 × 1280) onto the computational points of the C-grids using the bicubic interpolation code available also from Deutscher Wetterdienst.As resolution increases, the errors for this test case become very sensitive to time-stepping errors (J.Thuburn, personal communication, 2012) due to the shock of the initialisation.Therefore small gravity wave Courant numbers are used for all grids, as shown in Table 4.
The height contours and errors in comparison to the reference solution for some midresolution results are shown in Fig. 10.For this test case, the differences between using the centroidal and centroidal dual grids and between the symmetric and asymmetric H are tiny.Additionally, considering the different resolutions used for each grid type, the errors using each grid are very similar and the errors do not impact significantly on the total height pattern.The errors are plotted piecewise constant on every cell and so it is clear that there is no obvious grid imprinting or grid-scale noise apart from on the cubed-sphere using the asymmetric H, for which the cube edges are visible in the error fields.
The convergence with resolution for all grids of the 2 and ∞ errors of geopotential are shown in Fig. 11.The symmetric and asymmetric H results give nearly identical results and all grids show remarkably similar errors with convergence between first and second for 2 (φ) and first order for ∞ (φ).Due to the more undulating nature of the flow, the mountain test case does not suffer due to lack of grid alignment with the flow and so the problems, particularly with the cubed-sphere grid, do not show up.

Galewsky et al. (2004) unstable jet
The Galewsky et al. (2004) barotropically unstable jet is challenging because the instability can be released prematurely by truncation errors related to lack of grid alignment (Marras et al., 2014), grid inhomogeneity, or asymmetries in the discretisation.The test case is used with an initial perturbation and without viscosity.There is no analytic solution and so the results are compared with the STSWM reference model run at T426 using a time-step of 30 s and a hyperdiffusion coefficient of 4.97 × 10 11 m 4 s −1 .The relative vorticity after 6 days for the reference solution and for high resolution versions of the orthogonal HR grid and the centroidal grids are shown in Fig. 12.The results of the spectral model are interpolated onto HR grid 8 and the vorticity is plotted piecewise constant in exactly the same way as for the C-grid results.All of the C-grids use the asymmetric/diagonal H and the results using the symmetric H and/or centroidal duals are visually identical.All of the C-grid model runs in Fig. 12 use a time-step of 240s but this results in different ratios of initial advective to gravity wave Courant number for each grid (Table 4) since the initial jet is to the north of the smallest cells of the cubed sphere and diamond grids and the advective Courant number is based on the normal velocities and the initial, zonal velocities are not normal to any of the edges of the diamond grid.However the results are not very sensitive to the time-step.
The unstable jet using both version of the hexagonal grid is very similar to the reference solution except that the hexagonal grid results have spurious vorticity stripes upstream of steep gradients caused by phase errors of the vorticity when using the energy conserving version of TRiSK.In contrast, the results using the cubed-sphere or the diamond grid contain dramatic wave number 4 patterns which are not showing any signs of lessening with increasing resolution.This is in contrast to the results of Thuburn et al. (2013) for this test case, who use a higher-order advection scheme for pv and do not use an energy conserving Coriolis operator.If this test case is indicative of models that do not work well in 3-D as weather or climate forecasting models, then the cubed-sphere or diamond grids should not be used with this low-order differencing scheme.

Conclusions
A new C-grid discretisation of the shallow-water equations suitable for non-orthogonal grids has been proposed.Unlike the scheme of Thuburn et al. (2013), the new scheme does not rely on the dual grid being centroidal.This has advantages since centroidal grids, rather than grids with centroidal duals, often lead to lower errors.The new scheme formally loses energy conservation of the spatial discretisation but in tests, the energy conservation is very similar to the scheme of Thuburn et al. (2013).The new scheme will extend to three dimensions and so can be used for non-orthogonal grids over orography.
It has been demonstrated that the TRiSK perp operator is inconsistent (zeroth order accurate) but that this does not prevent convergence with resolution of shallow-water test cases.The perp operator leads to fluxes on the dual grid which have divergence which is a convex combination of the divergence on the primal and so an aspect of the perpendicular velocity is exactly correct.This helps the accuracy of the shallow-water test cases.
A new diamond grid of the sphere has been proposed which consists of quadrilaterals, is more nearly orthogonal and nearly as uniform as the equal angle cubed sphere.It mostly out-performs the cubed sphere in the tests undertaken.
The grids of quadrilaterals do not admit the computational modes of the hexagonal C-grid and hence enstrophy is better conserved.(Growth of the computational mode can lead to enstrophy increase whereas control of the computational mode can lead to enstrophy decrease.)However they do not outperform the hexagonal C-grid in any way and the hexagonal-icosahedral grid gives more accurate results in most test cases.However the lack of computational modes could be more beneficial in 3-D where the computational modes of the hexagonal C-grid could interact with, for example, the Hollingsworth instability.

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | are shown in Fig. 3 with the surface normal vectors, lengths and fluxes.Edge e of the primal grid has length p = |p|, normal vector p and tangential vector p ⊥ .Edge e of the dual grid has normal vector d and tangential vector d ⊥ .Here we restrict our attention to a low-order finite volume discretisation so that the volume (or area) flux across edge e is U e = v e • p and the circulation along dual edge e is V e = v • d ⊥ .Lower case variable names indicate values sampled at a point whereas upper case names are integrated values.Primal cells are indexed or denoted by i or j and dual cells are indexed or Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | However, numerical tests by Ringler et al. (2010) and further unpublished work have found the energy conserving version to generate more accurate solutions Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | without serious oscillations in pv.The non-energy conserving version is used by Thuburn et al. (2013).

)
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | can write down U e = H (V e ): U e = v e • p e (17) = d ⊥ e d e • p e V e + v e • p e − d ⊥ e d ⊥ e • p e .(18) Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | 23)are found using the method ofWeller et al. (2012) which involves multiple model runs with different initial conditions in order to evaluate the matrix that represents the model.The model uses f = 2k • Ω where Ω = (0, 0, 0.1) s −1 , H = 1 m and g = 1 m s −2 .The corresponding eigen values give us the amplification factors of the normal modes and are shown in Fig.4for coarse versions of the cubed sphere and diamond grids.Crank-Nicolson time-stepping with a timestep of 1 s is used with the Coriolis term treated explicitly and updated 4 times per time-step.
are assessed by reconstructing the velocities normal to each edge from the velocity component in the cell centre Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | U e v e • p e = HV Flux through primal face e V e v e • d ⊥ e Circulation along dual edge e

Fig. 1 .
Fig.1.Orthogonality and skewness of some grids.Non-orthogonality (from black, orthogonal, to blue) is shown on the primal edges whereas skewness of the primal grid (from red, no skewness, to yellow) is shown on the dual edges.

Fig. 2 .Fig. 3 .↓Fig. 7 .
Fig.2.The panels of the cubed sphere primal grid (dashed) and the diamond primal grid (grey) and the grid boxes for an equal angle cubed sphere grid.In the plane limit for the diamond grid, b = a sin 30 and c = a cos 30 implying c/b = √ 3.
It would seem logical to have an H operator with vanishing off-diagonal terms as the grid becomes orthogonal, regardless of the cell shapes.It would also be desirable to use an operator that is at least first order accurate on any grid.The first/second order operator for reconstructing cell centre vectors from normal components at edges in OpenFOAM (cited 2013) is:

Table 1 .
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Properties of the orthogonal hexagonal HR grid and the centroidal grids.

Table 3 .
Variables used in the discretisation and some of their finite-volume representations.Extensive quantities are upper case.Perpendicular with same magnitudeA i , A v Area of cell i, v A ivOverlap area between primal cell i and dual cell v A eArea associated with edge e A ieFraction of A e in cell i A iveFraction of A ie in dual cell v n ei sign (p e • (x e − x i )) Edge orientation indicator n ev sign (d e • (x e − x v ))