In this continuation paper to Starius (2018), referred to as (S-E-P) below, we consider numerical methods for the shallow water equations on the sphere, using the Equator-Pole (E-P) grid system. We make use of a reduced grid system for an annular band around the equator, which is divided into segments, each covered with a lat-lon net. As in (S-E-P) the equatorial grid is complemented by two orthogonal polar grids obtained by suitably modified stereographic coordinates. The whole grid system is optimised by minimising a global measure of uniformity, namely the maximum of all deviation factors on the sphere, where as in (S-E-P), the deviation factor for a net rectangle is the greatest ratio of its side-lengths.

The methods of this paper give an accuracy in the total mass, which is considerably higher than needed for NWP. This was already shown for the smooth Rossby-Haurwitz wave example for the method in (S-E-P). Here it is also verified, with and without reductions, for non-smooth solutions by considering the Cosine bell no. 1 and the mountain problem no. 5 in Williamson et al. (1992), the former with second order and the latter with first order discontinuous derivatives in their solutions. We have also integrated the mountain problem for a time period of 5 years, with maximal relative errors in the total mass of order 105. For the total energy the corresponding quantities were of order 102.

If K denotes the number of segments on a hemisphere, then there are 2 K-1 segments on the whole sphere, since the two segments closest to the equator are merged together. We note that K = 1 is the case considered in (S-E-P). By reductions the uniformity is improved and the total number of grid points can be decreased by 20% or more, with no substantial change in accuracy. For high resolution also the computing time will be reduced by a similar percentage. High order one-dimensional centred interpolation formulas are used to connect the discretisations on adjacent segments. The kind of deviations from uniformity that arise by changing the number of grid points on parallels do not cause any problems, at least not for the methods described in this paper.

We now briefly comment on some other methods considered in the literature. Our main sources are the two overview papers Staniforth and Thuburn (2012) and Williamson (2007). There are many other grid techniques than ours for the sphere, among which the most popular are the cubed sphere, the icosahedral and the Yin-Yang grid systems. The paper (S-E-P) contains a more detailed overview than the one given below.

Cubed sphere grids (Sadourny, 1972) consist of quadrilaterals and have been used for finite volume methods in, e.g. (Putman and Lin, 2007; Ullrich et al., 2010), for discontinuous Galerkin methods in, e.g. (Nair et al., 2005; Bao et al., 2014) and for finite difference methods in, e.g. (Ronchi et al., 1996). Some drawbacks are lack of orthogonality and deviation from uniformity, particularly at the 8 corner-points of the cube, where angles up to 120o instead of the ideal 90o occur.

Icosahedral grids (Williamson, 1968) are generally triangular and have been used for finite volume methods in, e.g. (Pudykiewics, 2011; Chen et al., 2014) and for discontinuous Galerkin methods in, e.g. (Giraldo, 2006; Läuter et al., 2008). At the 12 corner-points of the icosahedron there are angles up to 72o instead of the ideal 60o.

The more recently developed Yin-Yang grid, in Kageyama and Sato (2004), is an overset grid system consisting of two congruent orthogonal subgrids based on spherical coordinates and having maximal deviation factor 2. This grid has been used by, e.g. Qaddouri and Lee (2011) and for a full 3d model by Allen and Zerroukat (2016). Both this grid system and ours use overset and orthogonal grids, which leads to an obvious similarity between the two techniques.

The great importance of the axis of the Earth for the dynamics of the weather is clear. However, the Yin-Yang grid does not have good rotational symmetry properties relative to this axis, not even near the poles, for details see, e.g. Section 1 in (S-E-P). At the end of Section 3.4, we compare the Yin-Yang grids and ours with respect to the total number of grid points, both with the same number on the equator.

The paper is organised as follows. In Section 2, we introduce the set of governing equations, some notation, and also some coordinate systems. This section is essentially a summary of Section 2 in (S-E-P).

Section 3 is devoted to reduced grid systems for annular bands around the equator. Two methods, I and II, to divide the band into segments and then to connect them by interpolation, are considered in Sections 3.1 and 3.2. Method I is quite general and easy to implement but requires more interpolation coefficients than method II. However, method II has limitations in its applicability, which is discussed in Section 3.2. In Section 3.3, two different measures of uniformity, used to generate suitable grid systems for the sphere, are compared. In Section 3.4, an E-P grid system is constructed with square net rectangles at mid-latitudes, which is often required for NWP models.

In Section 4, we investigate our E-P grid systems numerically, both with and without reductions. Centred finite differences of high order are used for the spatial discretisations and the explicit fourth order Runge-Kutta method for the time integration. First some notation for various errors is introduced and then, in Section 4.1, smoothing by addition of a ’hyper-diffusion’ term is considered. For a discretisation method of order 2p a smoothing term of the same order is used, in contrast to the method in (S-E-P) where the order is always equal to 4. In Section 4.2, the rotation of the Cosine bell is considered and for this scalar advection example both the total mass and energy are obtained with much higher accuracy than for the height of the fluid. In Section 4.3, our results for an example by Läuter et al. (2005) show that the reduction techniques work quite well, the accuracy is very much the same for K = 1,2,3,4, particularly for the 6th order method. In Section 4.4, we consider zonal flow over an isolated mountain, example no. 5 in Williamson et al. (1992), with non-smooth solutions. The errors in the total mass shown in Table 6 are clearly small enough for NWP. We mention that no ‘corrections’ or ’fixers’ of any kind have been used.

In (S-E-P) further numerical examples are given, and also investigations of grid imprinting and comparisons between formal and computational orders of discretisation errors are considered.


Coordinate systems and governing differential equations

In this section, we introduce some notation and give a frame of reference for the paper. It is essentially a summary of Section 2 in (S-E-P).


Spherical coordinates and the shallow water equations

The transformation between Cartesian and spherical coordinates on the sphere is

{x=acosθcosλ,y=acosθsinλ, 0λ<2π, π2θπ2,z=asinθ, 
where a is the radius of the sphere, λ the longitude, and θ the latitude. The advective form of the shallow water equations in spherical coordinates is given by
where Ψ is the geopotential, f=2Ωsinθ the Coriolis parameter, Ω the rotation rate of the sphere, and u and v the eastward and northward velocity components, respectively.


Standard and modified stereographic coordinates

Consider first the standard stereographic coordinates. The relation between these and spherical coordinates can be expressed as

{xst=a q(θ)cosλ,yst=a q(θ)sinλ,
where q(θ)=2tan(φ/2), with the colatitude φ equal to 12πθ and 12π+θ, on the northern and southern hemisphere, respectively. We also need the map factor mf(θ)=2/(1+αsinθ), where α = 1 on the northern and −1 on the southern hemisphere.

Since the velocities in spherical coordinates, u=a cosθ dλ/dt,v=a dθ/dt, are not uniquely defined at the poles, we use instead the stereographic velocities U=dxst/dt,V=dyst/dt. Observe the use of the minus sign in the definition of U, intended to get nicer formulas, namely

(uv)=A(UV) and (uv)=A(UV),whereA=(sinλcosλcosλsinλ),
on the northern and southern polar regions, respectively. In (4), U and V have been scaled by the factor 1/mf(θ). The mapping (4) is also used in Starius (2014), where it is derived by orthogonal projection of the vector (u, v) onto the equatorial plane.

Stereographic projection is only suitable for grid generation close to the poles. Here we will use a modification of (3), already recommended in Section 2.3 of (S-E-P), so that the lat-lon and our new modified stereographic grid points will coincide on the meridian corresponding to λ = 0, and automatically also on the meridians λ=12π,π,32π, provided the number of grid points on parallels in the last segment of the reduced grid is divisible by 4. The assumption for λ = 0 leads to new coordinates (x,y) related to (xst,yst) by

{xst=2a tan(x/2),yst=2a tan(y/2). 

The modified coordinates (x,y) give rise to more uniform polar grids and will also simplify the connection between these grids and the last segments of the equatorial grid system. The shallow water equations expressed in the coordinates defined by (3) and (5), and the stereographic velocities U and V given in (4), are written down in full detail in Section 2.3 of (S-E-P).


Reduced latitude-longitude grid systems covering an equatorial band

In this Section, we consider the possibility to utilise a reduced grid system for an annular band around the equator, cf. Gates and Riegel (1963), Kurihara (1965), Tolstykh and Shashkin (2012), and Starius (2014). We also mention the reduced Gaussian grids for spectral models in Hortal and Simmons (1991), reported to decrease the computing time by about 20%, and used by ECMWF in their operational model. Further, the papers Layton (2002) and Li (2018) are worth mentioning in this context. The band is divided into segments, each covered with a lat-lon net. Let n be the given number of grid points on the equator and K the number of segments on a hemisphere. Here we choose to have invariance under reflection in the equator, where the closest two segments are merged together, which implies that the total number of segments is 2 K-1 and the total number of subgrids is 2 K + 1. We recall that K = 1 is the case studied in (S-E-P).

We will later define latitudes θ1<θ2<<θK such that the i:th segment Si is bounded by θi1 and θi, i=1,2,,K, where θ0=0. Also the number of grid points ni on a parallel in Si has to be defined and satisfy n1>n2>>nK, and the longitudinal steps are defined by Δλi=2π/ni, where i=1,2,,K. We note that n1=n and when only one longitudinal step is used, for example for K = 1, it is called Δλ.

We now discretise the segments introduced above. Assume that all the latitudes θi are multiples of a latitudinal step Δθ. Consider for each i, with 1iK, the parallels θi1+kΔθ,k=0,,si, where si=(θiθi1)/Δθ, and let Sid be the corresponding lat-lon net. Note that we have let Sid and Si+1d both contain the parallel corresponding to θi, which will ease the presentation in Sections 3.1 and 3.2. When the discrete segments have been defined, we omit the parallel θi from Si+1d, i = 1,…K-1, cf. Fig. 1.

Fig. 1

A small part 0λ45o and θ0 of a reduced grid for n = 160 and K = 3.

In this section, we will often use the fact that for a parallel with latitude θ and angular longitudinal step Δλ, the distance along the parallel between two consecutive grid points is Δλcosθ.


Method I: reduced latitude-longitude grid system

Here the number of segments K will be given by the user. The deviation factor at the equator is σK=Δλ1/Δθ, where Δθ will be the same for all discretised segments and determined later. We consider the segment Sid and require that the distances between grid points on the parallel θi1 are equal to Δλ1, that is Δλicosθi1=Δλ1, where i=2,,K. Further, the maximal deviation factor for the segment Sid is max(Δλicosθi1/Δθ , Δθ/(Δλicosθi) ), which is obviously minimised, with respect to Δλi/Δθ, when the two deviation factors in the max above are equal. We now have the equalities

Δλicosθi1=Δλ1 and Δλicosθi1Δθ=ΔθΔλicosθi,
and substituting Δλi in the second equality by Δλ1/cosθi1 and using σK=Δλ1/Δθ, gives

The transition from SKd to a polar grid gives that σK=mf(θK)=2/(1+sinθK), where mf(θ) is the mapfactor in Section 2.2, cf. (12) in Section 3.2 of (S-E-P). This equality and the right hand side of (7) imply

which can be solved by, e.g. Newton’s method. We then determine σK=2/(1+sinθK) and by using cosθi=cosθi1/σK2=σK2i we find θi, i=1,2,,K−1. Note that (8) implies that θK12π, and hence σK1, as K.

Now Δθ=Δλ1/σK can be evaluated and for practical reasons modified by defining the integer m to be 12π/(Δλ1/σK) rounded to the nearest integer and then redefining Δθ=12π/m. With Δθ known we now round the latitudes θi, i=1,2,,K, to their closest multiples of Δθ. Finally we set ni=ncosθi1, round it to the nearest integer, and define Δλi=2π/ni,i=2,,K.

In Table 1, the overlap between the last segment in the reduced grid and the polar grid is 3, which only affects the number of parallels in this segment. An increase in the number of segments leads to a decrease in the total number of discretisation points and in the maximal deviation factor. In Fig. 1, a reduced grid with n = 160 and K = 3 is illustrated. More realistic resolutions for NWP are for instance n = 1000, K = 10, which gives p = 20.7% and s10=9, and n = 2000, K = 10, which gives p = 20.8% and s10=16. We recall that the above procedure requires only one-dimensional interpolation and that the coefficients in the formulas should be computed and stored during preprocessing. The number of sequences of interpolation coefficients is n1+2(n2++nK1)+nK, which for high resolution is negligible compared to other storage requirements. A less general method using considerably fewer interpolation formulas is briefly described in Section 3.2.

A way to avoid scalability problems for the couplings between the reduced grid and the northern polar grid, say, is to let the former reach a high latitude. Since θK increases slowly with K (for instance K = 10 gives only θ10=65.8o) we will now consider a modification of method I, which is suitable for high latitudes. We choose K and determine the first K segments and the corresponding maximal deviation factor σK, and then complement this by the desired number of additional segments, each optimised with respect to σK. More specifically, additional latitudes θK+1,θK+2, can be determined by the equation cosθi=σK2i,i>K, and then also ni=ncosθi1 and Δλi=2π/ni,i>K, where the ni have been rounded in a suitable way.

We conclude this subsection with two examples. For initial K = 3, we have σ3=1.0988 and the corresponding latitudes are: 34.0, 46.8, 55.1, 61.9268, 67.0595, 71.1657, 74.4912, 77.2050, 79.4305, 81.2615. If we want the reduced grid to reach the polar circles, then we add two segments, so that K = 5, and reach latitude 67o. The maximal deviation factor for the whole reduced grid system is still σ3. In a second example, with initial K = 1 and σ1=1.1892, two segments have to be added, so that K = 3, for the reduced grid to pass the polar circles. There is obviously a trade-off between uniformity and the possibility to use a small number of segments.


Method II: reduced latitude-longitude grid system

This method is essentially applicable with a small number of segments and special choices of n. However, its range of applicability can be extended in ways described below. The notation and assumptions of Section 3.1 will be used here, and only a brief outline of the method will be given.

The method is based on the following simple idea. In a consecutive manner we replace q longitudinal mesh intervals in segment Sid with q−1 intervals in the following segment Si+1d. As for method I we require that Δλicosθi1=Δλ1, where i=2,,K, and further that each ni, with i < K, is divisible by q. This implies that Δλi+1/Δλi=ni/ni+1=q/(q1). Further, Δλicosθi1=Δλi+1cosθi together with (7) leads to Δλi+1/Δλi=cosθi1/cosθi=σ2, where we write σ instead of σK, since σ is now independent of K. We choose q and determine σ and ϑ0, where ϑ0 is the width of the reduced grid, by using

σ2=qq1 and 21+sinϑ0=σ.

We now determine θi by cosθi=σ2i,i=1,2, until θi>ϑ0. Further, we can choose K from |θKϑ0|=mini|θiϑ0| and then redefine ϑ0=θK. For low resolution it might be necessary to replace K by a smaller value in order to get a reasonable number of parallels in the segment SKd. Since n1=n is known we can use ni+1=(q1)ni/q,wherei=1,,K-1, and then define Δλi=2π/ni,i=2,,K.

The above works only if n is divisible by qK1, where K has our desired value, in which case the number of sequences of interpolation coefficients required is only q−1, provided symmetry is used to reduce the number of sequences whenever possible.

If we get a number ni0, with i0<K, not divisible by q, then λ=2π will not correspond to grid points in Si0+1d. Since we now cannot have equidistant grid points in all our discretisations, we switch to a modification of method II, which allows that. From segment Si0+1d onward, we define ni+1=(q1)ni/q, for i=i0,,K and then round them to integers and define Δλi+1=2π/ni+1, for i=i0,,K. This transition to method I is very simple but requires more sequences of interpolation coefficients.

For high resolution, we indicate an improvement of the above modification of method II. If ni0+1 is rounded to an integer divisible by q, then the interpolation between Si0+1d and Si0+2d can be done with the q−1 interpolation formulas. We still use method I for the interpolation between Si0d and Si0+1d. This can be repeated by choosing ni0+3 divisible by q, using method I between Si0+2d and Si0+3d, and then the q−1 interpolation formulas between Si0+3d and Si0+4d, etc.

We conclude this subsection by considering an example with q = 6, which preliminarily gives σ=6/5=1.0954. We have ϑ0=55.7 and the following sequence of latitudes determining the segments: 33.6, 46.0, 54.6, 61.2, 66.3, 70.4, 73.8, 76.6, 78.8, 80.7. Thus according to the above we choose K = 3, redefine ϑ0=θ3, and obtain a new σ=2/(1+sinθ3)=1.1019. If n is divisible by 62, then method II works perfectly well for K = 3. However, if we want K > 3 and n is not divisible by 63, we can use a modification as described above for the segments following S3d.


Comparison between two measures of uniformity

The measure of uniformity most commonly referred to in the literature is the ratio of maximum to minimum grid length, which we call γ in this section. Our main interest here is to investigate which measure of uniformity is best suited for grid generation. Our maximal deviation factor has already been used for this purpose, and we therefore now turn our focus to the ratio γ by studying two examples.

Example (i): Let Sd be a segment on the northern hemisphere bounded by the latitudes θ1<θ2 and consisting of a lat-lon net with the longitudinal step Δλ. We define σ as the deviation factor for net rectangles at θ1, that is σ=Δλcosθ1/Δθ, where Δθ is the latitudinal step, not yet determined. The ratio of maximum to minimum grid length on Sd is

γ=γ(σ)=Δλcosθ1min( Δθ,Δλcosθ2 )=σmin(1,σcosθ2/cosθ1 ),
from which it follows
γ= {cosθ1/cosθ2,σcosθ1/cosθ2,σ,σcosθ1/cosθ2, 

The minimum of γ is cosθ1/cosθ2, which is attained for all σcosθ1/cosθ2. The result is unsatisfactory since we want the grid to be uniquely defined.

For comparison we mention that the minimum of the maximal deviation factor for this segment is cosθ1/cosθ2, which is attained if and only if Δλ/Δθ=1/cosθ1cosθ2. We get square net rectangles at the latitude θ for cosθΔλ=Δθ, which implies


This formula will be used in the next Section 3.4. ⊗

Example (ii): In this example the ratio γ will be minimised over the whole sphere. We choose the E-P grid system with K = 1, which we optimised in (S-E-P) by using the maximal deviation factor. Now instead we minimise γ. Let the equatorial grid cover the band |θ|ϑ0, which will be determined later. Our northern polar grid will go down to ϑ0 and the corresponding map factor is mf=2/(1+sinϑ0). According to (12) in (S-E-P) the greatest and smallest grid lengths in our polar grids are Δθ and Δθ/mf=Δθ h(ϑ0), respectively, where h(ϑ0)=(1+sinϑ0)/2. For our polar grids the ratio γ and the maximal deviation factor are the same.

The ratio γ for the whole sphere can now be written as

γ=Δλmin( Δλcosϑ0,Δθh(ϑ0) )=1cosϑ01min( 1,h(ϑ0)/(σcosϑ0) ),

where σ=Δλ/Δθ. We now consider the equation h(ϑ0)/(σcosϑ0)=1, with a unique root θ*(σ) in the interval (0π/2). Elementary computations show that

cosθ*(σ)=4σ/(4σ2+1) and θ*(σ) is increasing.

We look at two cases, namely


In both cases γ(σ,ϑ0)1/cosθ*(1)=1.25, with equality if and only if σ = 1 and ϑ0=θ*(1)36.87o. This result was also reported in (S-E-P) but without mentioning the condition σ = 1, which we regret.

The sole reason why lat-lon coordinates are not suitable for the whole sphere is the convergence of meridians when we are moving towards a pole. Therefore, we find it inappropriate to use square net rectangles at the equator instead of rectangles with longer longitudinal than latitudinal sides. ⊗

The examples (i) and (ii) above convincingly show that the maximal deviation factor is a better tool for generation of suitable grid systems on the sphere, than the ratio of maximum to minimum grid length.


Square net rectangles at mid-latitudes

Grids for NWP are generally constructed to be uniform (having square net rectangles) at mid-latitudes, θ=±π/4, a condition we will study in this section. We will also compare the total number of grid points used in the Yin-Yang grids and ours, both with the same number of grid points on the equator, and by considering the mid-latitude condition above.

We first use method I of Section 3.1 to determine the segment Si0d with θi01<π/4<θi0 and then use (9) with θ=π/4, which leads to the equation


For our ’mid-latitude method’ we use notation with bars to distinguish the variables from those used for methods I and II. In order to avoid over-determined systems of equations we replace equation (8) by (12) and use the left part of (7) to get cosθ¯i0=σ¯2i0, which together with (12) gives σ¯4i0+2=1/2. We then determine σ¯ and θ¯i,i=1,,K, by solving the equations


We conclude this subsection with an example, in Table 2, where two reduced lat-lon grids are considered. Choose n = 3600 and K = 10 in method I, which corresponds to Δλ=0.1o on the equator and about 11 km between neibouring grid points everywhere on the sphere.

For the Yin-Yang grid we let the Yin subgrid have square net rectangles at θ=±45o. If the Yang grid is congruent with the Yin grid the total number of grid points will be (2·2·3/16)n2, which for n = 3600 is about 6.873·106. This grid system only satisfies the mid-latitude condition in a crude sense. By replacing this grid by one of ours in Table 2 we get a reduction in the number of grid point of about 34%. Moreover, we get good rotational symmetry properties relative to the axis of the Earth, much better uniformity, and perfect possibilities to satisfy the mid-latitude condition in suitable ways.

Provided we accept less uniformity than in Table 2, a smaller number of segments can be used, e.g. K = 3 or 5. If we also want the reduced grid to reach high latitudes it can be complemented by additional segments by the procedure described at the end of Section 3.1.

For a full lat-lon grid, satisfying the mid-latitude condition and with n = 3600, the total number of grid points is about 9.164·106.


Numerical experiments

In this section, we investigate the E-P grid system numerically, both with and without reductions in the equatorial grid. Test examples no. 1 and 5 from Williamson et al. (1992) and time dependent flow from Läuter et al. (2005) are considered, cf. also (S-E-P) where other numerical examples and investigations, such as grid imprinting, are studied. We recall that the number of segments in the equatorial grid is 2 K-1, that the total number of grids in the E-P system is 2 K + 1, and that K = 1 corresponds to the grid system used in (S-E-P). Method I, intended for reduction and connection between segments, has been used in all the tables below. However, in Sections 4.3 and 4.4 we give some short comparisons between methods I and II.

The spatial discretisation is obtained by using centred finite difference approximation of orders 6 or 4 on non-staggered grids, for which formulas can be found in Section 2.4 of (S-E-P). The classical explicit fourth order Runge-Kutta method is used for the time integration, which requires interpolation between grids for each of the four stages in a time step.

Explicit time integration seems to work quite well for the shallow water equations. However, for realistic NWP, 3d models are used and with vertical spacing much smaller than the horizontal ones, which would lead to very small time steps for explicit methods. This is why HEVI(horizontal explicit and vertical implicit) schemes are of interest in the atmospheric community.

In the tables below, 2p is the approximation order of the spatial discretisations of the underlying differential equations. The order of the interpolation formulas connecting the polar and equatorial grids is 2p + 1, and between discretised segments the order used is 2p + 2. The overlapping number is denoted by s (=p), cf. Section 3.3 in (S-E-P).

The geopotential Ψ is defined as gH, where g is the acceleration of gravity and H the height of the fluid. We will often use the relative error norm 2(H), introduced in Williamson et al. (1992). Let 2(H,d) denote this norm after d days. For the vector (2(H,1),2(H,2),,2(H,T)), where T is the integration time in days, we define a mean value by


The global conservation properties for mass and energy will also be studied below. The integral formulas can be found, e.g. in Starius (2014) and their numerical computation is described in (S-E-P). The relative errors for the total mass after each day are stored in the vector (rm(1),rm(2),,rm(T)), with mean value r¯m(T), defined by (14). The quantity r¯e(T) for the total energy is defined analogously.



Smoothing is essential for the centred finite difference methods to work properly and will be achieved by adding a so called ‘hyper-diffusion’ term to each equation. Since ‘hyper-diffusion’ has no known physical meaning, we consider our smoothing a purely numerical stabilisation process.

Let Lh be the five-point operator associated with minus the Laplacian 2, and defined by


To the first, second, and third equation in (2), we add the terms cLhpu,cLhpv, and cLhpΨ, respectively. In the equations for the polar grids, U and V are used instead of u and v. Since our grid systems are almost uniform, it is natural to use Lh for the whole sphere. We observe that LhpΨ=O(Δθ2p), for any smooth function Ψ. In the temporal discretisations the smoothing terms will be multiplied by Δt. The smoothing constant c0 will often be determined by minimising 2(H,T),¯2(H,T),r¯m(T) or r¯e(T). We have used the MATLAB® routine fminbnd, for the numerical minimisation.

In an earlier work we have used Lh as a five-point difference approximation to minus the Laplacian in spherical coordinates, without any improvement for almost uniform grids.

We point out that the optimisations in this paper, are essential in order to ease comparisons between different methods.


Rotation of the Cosine bell

In this example we consider a solid body rotation of the Cosine bell, according to test problem no. 1 from Williamson et al. (1992). The advecting wind is given by

where α is the angle between the axis of the solid body rotation and the polar axis, and u0=2πa/12 m/day.

The E-P grid system will be used for K equal to 1 and 3 and the smoothing parameter c will be determined by minimising 2(H,12),r¯m(12) and r¯e(12). The results are given in Table 3, first for 2p=6 and then for 2p=4.

Despite the solution having discontinuous second order derivatives, the accuracy of H is considerably higher for the case 2p=6 than for 2p=4. The actual approximation orders can still be fairly equal. The accuracy for the total mass must be caused by cancellation of errors since it is much higher than for H. In this specific example, the same can be said about the total energy. Minimisation of 2(H,12) leads to slightly more accurate values of H, for K = 1 than for K = 3, but the values of r¯m(12) are much smaller for the smoother case K = 3, which also require fewer grid points.

For the rest of this subsection we only consider the most important case 2p=6. The minimisation of r¯m(12) leads to much greater values of c than for the other alternatives, which results in greater values for 2(H,12) and r¯e(12), see for example the case K = 1. For this reason we reject the use of r¯m(12) for the determination of c and note that r¯e(12) works better for this example if we focus on the other errors 2(H,12) and r¯m(12) and not on the values of c themselves, which are much smaller compared to the ones corresponding to 2(H,12).


Time dependent flow, the Läuter example

The following time dependent flow is a solution of the shallow water equations and was introduced in Läuter et al. (2005), cf. also Pudykiewics (2011),

where u0=2πa/12 m/day, k1=133681 and k2=10m2/s2. We will compute Ψ*=gH*, where H* is the height of the fluid. We note that the solutions are periodic in time with period 1 day, which means that the problem corresponding to (15) is periodic in all its variables. In Kreiss and Oliger (1972) it is shown, for a simple periodic advection model, that the discretisation error increases linearly in the number of time periods. This behaviour in time of the discretisation error is the best we can hope to achieve.

We make some comments concerning Table 4. For order 2p=6 the errors decrease somewhat from the case Δt=300 to 200. This indicates that there is no real cancellation between the spatial and temporal local discretisation errors. For order 2p=4 there seems to be cancellation, and more so for Δt=300 than for Δt=200. For 2p=6 the accuracy obtained for K = 1,2,3,4 is very much the same, particularly for Δt=300. For each case the error grows very close to linearly in the number of days.

The third row of Table 4 will be (300,3,12.3,6,2.1010e9,4.1408e9,6.2201e9), if method II is used instead of method I. For this very smooth example there is no real difference in accuracy between the two methods.

In Table 5, r¯e(15) is minimised in order to investigate whether it may be used instead of ¯2(Ψ*,15), since the latter quantity is not easily available for realistic NWP. The values of c in the second to fourth rows in Table 5 are much smaller than those in Table 4, but the corresponding errors do not differ very much, except for the two cases K = 3, 4 and T = 15 days. However, the minimisation of r¯e(15) might still give some insight in certain situations. In the next section, we will look at another way to determine values of c.

The solutions (15) are very smooth and have very small variation, which explains the high accuracy for the geopotential, and for the total mass and energy cancellation of errors has further increased the accuracy, see Table 5.


Zonal flow over an isolated mountain

Test problem no. 5 in Williamson et al. (1992) is concerned with zonal flow impinging on a mountain on the sphere, with a circular base. The centre of the base is denoted by (λc,θc) and the initial conditions are

where u0=20 m/s and h0=5960 m. The geopotential corresponding to the conical mountain is
and r=min((λλc)2+(θθc)2,R). We have computed Ψ=Ψ*+Ψs, which occurs in the two momentum equations, whereas Ψ*=ΨΨs appears in the continuity equation. After an initial transient phase there seems to be fairly well-behaved solutions, but with discontinuous first order derivatives on the circle limiting the base of the mountain and thus less smooth than the Cosine bell in Section 4.2. Contour curves for the total height H are depicted in Fig. 2.

Fig. 2. 

Contour curves for the total height H of the fluid for the mountain problem, no. 5 from Williamson et al. (1992): n = 360, K = 3, 2p = 6, c=105, T = 15 days. Method I.

It is fairly easy to experimentally determine a ’stability interval’

where cmin guarantees that there is enough smoothing for stability and cmax that we stay inside the stability region of the time integrator, in the left part of the complex plane. We observe that for implicit methods cmax can be infinite. For the present example we can choose cmin=0.1×106 and cmax=18×106 for 2p=6, and cmin=0.5×106 and cmax=140×106 for 2p=4.

In Table 6, the conservation properties are about the same for 2p=6 and 2p=4, which is natural since they are caused by cancellation of errors in the height of the fluid rather than by accuracy. We recall that q is the number of longitudinal mesh intervals that is replaced by q−1 intervals in the following segment, cf. Section 3.2. Only q−1 interpolation formulas are required and repeated periodically. This can be seen as a more regular way to interpolate between segments than the one of method I, which has resulted in greater accuracy for the total mass. The values of c in Table 6 have been obtained by using some tests.

In order to further investigate the reliability of our Equator-Pole grid systems, we have integrated the mountain problem for 5 years, using methods of order 6. In Table 7, the maximal relative errors in magnitude of the total mass appear in the fourth column and the corresponding quantity for the total energy in the last one, both after 5 years. The values of c are taken from Table 6.



The Equator-Pole grid system introduced in (S-E-P), for the solution of the shallow water equations on the sphere, consists of one equatorial and two polar grids. Here we generalise this by permitting the equatorial grid to be reduced, that is divided into segments. We give two methods of reduction, called method I and method II, in Sections 3.1 and 3.2, respectively. The reduction technique seems to work quite well, as can easiest be seen in Table 4.

By using reductions the total number of grid points on the sphere can be decreased by 20% or more, and for high resolution also the computing time by a similar percentage. The change in accuracy for the height of the fluid is negligible, and for the total mass the accuracy is increasing with increased uniformity of the grid system, cf. Table 3. High uniformity might have additional advantages, for example in conjunction with comprehensive wave propagation. The different wave speeds computed will, among other things, depend on the uniformity of the grid system used and its resolution. Low uniformity or resolution might lead to interaction between waves with no analogue in the continuous model.

We have demonstrated high accuracy in the total mass for fairly simple problems, both with and without smooth solutions. For realistic NWP models the situation can be more complicated, e.g. water can contribute to the atmosphere’s mass and give rise to far less smooth problems than those considered here. We mention that Allen and Zerroukat (2015) have derived a two-dimensional generalisation of the shallow water equations, which includes moist.

We use high order centred finite difference approximations of first order spatial derivatives, in the shallow water equations. This requires only a fraction of the work per grid point compared with some other methods of the same order. This fraction decreases rapidly with the order of approximation. Further, one can very easily increase the order of approximation by increasing the order of approximation for first order derivatives. For methods of order 4, 6, 8 and 10 the corresponding spatial stencils contain only 9, 13, 17 and 21 grid points, respectively. The higher the order of the spatial approximation, up to a bound depending on the resolution, the better treatment of fronts in the solutions, see for instance Table 2 in (S-E-P). For the time integration, we have used the explicit fourth order Runge-Kutta method.

The Equator-Pole grids presented here only focus on uniformity, except in Section 3.4, and are mainly intended to be basic grid systems. To make them more flexible, for instance when taking jet streams and orography into account, they need to be modified, perhaps by using overlapping techniques.