In this article, we consider numerical methods for the shallow water equations on the sphere. Our particular approach is based on a system of three grids, together covering the entire sphere, which we call the Equator–Pole (E-P) grid system and which is required to have an optimal uniformity property. The basic grid is an ordinary spherical latitude-longitude grid placed in an annular band around the equator, $-45\xb0\le \text{latitude}\le 45\xb0$, which constitutes about 71% of the sphere. This equatorial grid is complemented by two other grids in the polar regions, obtained by using suitably modified stereographic coordinates. The grids are orthogonal and also almost uniform, by which we mean that for a net rectangle the maximal ratio of its side-lengths, the deviation factor, is only slightly greater than 1, for our grid system always less than 1.19. The connection between the grids is achieved by overlapping and the use of centred interpolation formulas. Overlapping techniques, sometimes called overset or Chimera grid systems, appeared originally in Starius (1977a, 1980), Volkov (1968) and later for instance in Browning et al. (1989) and Henshaw and Schwendeman (2006). We focus on keeping the overlapping to a minimum so that the cost for the connections will only be a small fraction $O(1/N)$ of the total computing time, where N is the number of gridpoints on the equator.
Our proposed E-P grid system is based on a family of grid systems depending on two parameters and is chosen to minimise the maximal deviation factor mentioned above over these two parameters. A simplified overview of the E-P grid system can be seen in Figs. 1 and 2. Our grid system shows similarity with Phillips (1959).
We will 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). Let us first note that according to Staniforth and Thuburn (2012), most operational global NWP models are based on full latitude-longitude grids and semi-Lagrangian schemes (Robert, 1982; Lauritzen et al., 2010). By replacing a full latitude–longitude grid with our E-P grid system, the so called pole problem disappears and the number of gridpoints would be reduced by a factor of about two-thirds $(2/\pi )$.
There are many other grids and grid systems for the sphere, among the most popular are the cubed sphere and the icosahedral grids. Cubed sphere grids (Ronchi et al., 1996; Sadourny, 1972) consist of quadrilaterals and have been used for finite volume (Putman and Lin, 2007; Ullrich et al., 2010), discontinuous Galerkin (Bao et al., 2014; Nair et al., 2005), and finite difference methods (Ronchi et al., 1996). Some drawbacks of this grid system are lack of orthogonality, and deviation from uniformity. At eight points on the sphere, corresponding to the corner-points of the cube, there are angles up to 120° instead of the ideal 90°. Icosahedral grids (Williamson, 1968) are generally triangular and have been used for finite volume (Chen et al., 2014; Li and Xiao, 2010), and discontinuous Galerkin methods (Giraldo, 2006; Läuter et al., 2008). At 12 points there are angles up to 72° instead of the ideal 60°. We see no reason here to consider in detail the implications of these irregularities, consult (e.g. Peixoto and Barros, 2013). In Pudykiewics (2011), it is reported (in the Section 1) how important it can be to perform some kind of suitable grid regularisation, in order to improve the convergence rate of various quantities.
The great importance of the polar axis for the dynamics of the weather is indisputable, therefore we should require that grid systems have good symmetric properties relative this axis, as in Fig. 2. In Staniforth and Thuburn (2012), Section 3.5.2 is a grid system the authors call modified Yin–Yang, that at first glance bears resemblance to our E-P grid. It consists of an ordinary lat(θ)–lon(λ) grid, with $-45\xb0\le \theta \le 45\xb0$ and $0\xb0\le \lambda \le 360$, and two polar grids, which can be obtained by moving the spherical quadrilateral ql of the equatorial grid corresponding to $-45\xb0\le \theta \le 45\xb0$ and $0\xb0\le \lambda \le 90$ to the north and to the south by $90\xb0$, respectively. The ql is not equilateral, the ratio between the longest and shortest edges is $\sqrt{2}$, which is also reflected in the shapes of the mesh rectangles. We refer to Figure 8(b) in Staniforth and Thuburn (2012), where the lack of symmetry relative the polar axis is obvious. For a three grid system these polar grids should not be used, since there are much better alternatives, cf. Section 2.3 and Fig. 1.
We also want to mention the interesting, overset Yin-Yang grid system (Kageyama and Sato, 2004; Qaddouri and Lee, 2011), developed in 2004. In a future paper, we plan to compare this with our E-P grid system. For the latter we will also experiment with letting the equatorial grid be a suitably reduced grid (cf. Gates and Riegel, 1963; Kurihara, 1965), in order to make it broader and to increase the uniformity. We admit that our E-P grid system does not quite share the classical and philosophical elegance of the Yin–Yang grid, but considering the above it might be more functional.
The article is organized as follows. In Section 2, we first express the shallow water equations in spherical and later in modified stereographic coordinates. Further, some centred difference approximations for first order derivatives are given.
In Section 3, the E-P grid system is defined and details concerning the underlying minimax problem are given. Further, overlapping and interpolation between grids with minimal overlapping are considered, depending on what we call strict overlapping, cf. Section 3.4.
Section 4 is devoted to numerical experiments. To evaluate our grid system we only consider the well-known high order centred difference methods, because of their low complexity. We present two groups of experiments, namely advection problems and problems for the non-linear shallow water equations. In the first group, we consider smooth deformational flow, and solid body rotations of the Cosine bell. The second group contains test problems 2, 3, and 6 from Williamson et al. (1992). For problems 2 and 3 we investigate whether any harmful grid imprinting occurs, and the answer is no. For problem 6, we study conservation properties for mass and energy.
The development of new methods generally involves many people and research publications. Thus, we will not state that our E-P grid systems is superior to, for example, the cubed sphere grid for NWP. However, the following list of advantages might be of interest: (i) Orthogonality and better uniformity, (ii) 3 grids instead of 6, (iii) the ideal latitude–longitude grid is used for a major part of the sphere, (iv) no harmful grid imprinting, cf. the test problems 2 and 3 in Section 4, (v) in Table 5, the Rossby–Haurwitz wave, the total mass is correct to about seven decimal digits, which is obviously much more than needed for NWP, (vi) possibility to let the equatorial grid be a reduced grid to increase the uniformity and by connecting to the polar grids sufficiently close to the poles, scalability problems can be avoided, and (vii) the computational results for the E-P grid system compare favourably with other alternatives, cf. Section 4.
Efficient parallel implementation of global methods for NWP are of course very important and are mostly performed by experts in HPC. The standard decomposition of the three grids in our E-P system is straight forward. The scalability considerations for the E-P and the Yin–Yang grid systems are very similar. In Kageyama and Sato (2004), the latter is used for mantle and geodynamo simulations on a parallel computer of distributive memory type, with 5120 vector-type processors. The conclusion of the authors in Kageyama and Sato (2004) is: ‘The Yin–Yang grid is suitable for massively parallel computers’. In this article, only the simplest version of an E-P grid system is considered in detail. However, the generalisation mentioned in (vi) above is probably more interesting and can completely avoid scalability problems in relation to the overlapping. Some additional aspects of implementation are given at the beginning of Section 4.
Throughout this article, we assume that the solutions, with the velocity components appropriately expressed, are smooth on the entire sphere, in local Cartesian coordinates.
The transformation between Cartesian and spherical coordinates is
In our numerical experiments the polar regions $-\pi /2\le \theta \le -\pi /4$ and $\pi /4\le \theta \le \pi /2$, will be covered with grids corresponding to somewhat modified stereographic coordinates, considered in detail in Section 2.3 below. Let us first look at the standard stereographic coordinates, for which many formulas can be found in, for example Lanser et al. (2000). We express the relation between the two systems of coordinates as
Since the velocities in spherical coordinates, $u=a\hspace{0.17em}\text{cos}\hspace{0.17em}\theta d\lambda /dt$, $v=ad\theta /dt$, are not uniquely defined at the poles, we use instead the stereographic velocities $U=-d{x}_{st}/dt,V=d{y}_{st}/dt$. Observe the insertion of the minus sign in the definition of U, intended to get nicer formulas, namely
For stereographic coordinates the distances between consecutive gridpoints are decreasing, when we move towards the equator along a meridian, and the reduction factor from a pole to a parallel with latitude θ is $1/{m}_{f}(\theta )$. From the pole to the equator the distances between pairs of neighbouring gridpoints decrease by 50%. Stereographic projection is only suitable for grid generation when used quite locally. We therefore recommend a modification of (3) so that the spherical (1) and our new modified stereographic gridpoints will coincide for $\lambda =0,\frac{1}{2}\pi ,\pi $, and $\frac{3}{2}\pi $. This leads to more uniformity for the polar grids and will also simplify the connection between the equatorial and the modified stereographic grids. Since these will have common gridpoints on the meridians referred to above, it will be easy to define overlapping quantitatively. Further, the number of gridpoints in the polar regions will decrease and the amount of interpolation between grids can be somewhat reduced.
Let us introduce new coordinates $(x\prime ,y\prime )$, instead of $({x}_{st},{y}_{st})$, by
By choosing $x\prime =\phi $ in the formula above and making a similar choice for y_{st}, we arrive at
We now turn to the shallow water equations expressed in the coordinates defined in (3) and (6) and the stereographic velocities U and V, given in (4). From the shallow water equations in stereographic coordinates, cf. Lanser et al. (2000) or Browning et al. (1989), we easily obtain the equations in our modified coordinates
As earlier α = 1 on the northern hemisphere and −1 on the southern.
In our numerical experiments, given in Section 4, the following discretisations will be used. Replace first order derivatives by centred equidistant difference approximations of the wanted order 2p, say, in (2) and (7).
With ${x}_{j}=jh$, for some meshlength h, and ${f}_{j}=f({x}_{j})$, we write down a couple of approximations, for first order derivatives, with error terms
We assume that the same physical meshlength is suitable everywhere, and that thus the grids should be as uniform as possible.
In the following the three grids in the E-P system will be referred to as: (i) The equatorial grid, (ii) the northern polar grid, and (iii) the southern polar grid. We will often consider only the northern hemisphere, because of symmetry with respect to the equator. Further, quadrilaterals on the sphere with right angles will be called rectangles or squares.
The E-P system will be uniquely determined by specifying the number of gridpoints, N say, on circles of parallel in the equatorial grid. We define $\Delta \lambda =2\pi /N$ and assume, because of symmetry, that 4 is a divisor of N. Further, denote the latitudinal meshlength by $\Delta \theta $, and let θ_{0} be the latitude on the northern hemispheres, where the grids meet. In the present section $\Delta \theta $ and θ_{0} will be determined by minimising the maximal deviation factor, defined in detail in Section 3.1, for the E-P system.
On a parallel with latitude θ, the distance between two consecutive gridpoints is $\Delta \lambda \hspace{0.17em}\text{cos}\hspace{0.17em}\theta $, provided the radius of the earth is used as the unit of length. In the λ and θ directions, the distances from a gridpoint to its closest neighbours should have a ratio as close to 1 as possible. For this reason we define the deviation factor for latitude θ as
We require σ to be greater than 1, since this will lead to a smaller maximal deviation factor, see Fig. 1. An immediate consequence is that $d(0)=\sigma $.
Let us now in detail consider the equatorial grid covering the region $|\theta |\le {\theta}_{0}$, where the value of θ_{0} will be chosen later. Because of monotonicity we have
Let us first define ϱ as the quotient $\ge 1$ of the areas of two net rectangles and ${\varrho}_{max}$ as its maximal value in a grid. Further let us again consider stereographic projection, which is a conformal mapping. The latter can be expressed for instance by saying that the mapping is locally a similarity transformation, or that the Jacobian of the transformation is an orthogonal matrix, multiplied by a scalar. This implies that the deviation factors are close to 1, but the value of ${\varrho}_{max}$ can be considerable, e.g. ${\varrho}_{max}=4$ for a hemisphere. Conformal mappings are often too inflexible for grid generation (cf. Starius, 1977b). In our modified stereographic coordinates the orthogonality is preserved and the uniformity properties are improved, in the sense that the values of ϱ become smaller. We do not actually use ϱ in the design of the E-P grid system.
For an orthogonal net only the deviation factors and the quotients ϱ are needed to measure uniformity. In our equatorial net the meshlength $\Delta \theta $ is kept constant, which implies that ${\varrho}_{max}={d}_{\mathit{\text{minmax}}}^{2}$. Thus it suffices to minimise the maximal deviation factor in our equatorial net.
For the modified stereographic coordinates in Section 2.3, we choose $\Delta x\prime =\Delta y\prime =\Delta \theta .$ On the meridians $\lambda =0,\frac{1}{2}\pi ,\pi $, and $\frac{3}{2}\pi $, this will give equidistant gridpoints with distance $\Delta \theta $. The distance in the orthogonal direction will then be about $\Delta \theta /{m}_{f}(\theta )$, because for $y\prime $ of order $\Delta y\prime $, the stereographic y_{st} will differ only slightly from $y\prime $, cf. (6). At meridians $\pi /4+k\pi /2,k=0,1,2,3$ the net rectangles in the modified grid are almost squares and the deviation factors on a circle of parallel are maximal for $k\pi /2,k=0,1,2,3$. Thus, the maximal deviation factor in the northern polar grid is
The optimal value of θ_{0} can therefore be determined by the equation, cf. Equations (11) and (12),
Let us briefly mention the frequently used ratio of maximum to minimum gridlength as a measure of uniformity. For our E-P grid system this ratio is minimised by using the equation $1/\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta}_{0}=2/(1+\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta}_{0})$ with solution $\text{cos}({\theta}_{0})=0.8$, $180{\theta}_{0}/\pi \approx 36.87\xb0$. The ratio is 1.25, which is small compared to other grids, cf. Section 3.3 in Staniforth and Thuburn (2012). For orthogonal meshes we prefer our own deviation factor (9), partly because it leads to smaller polar grids.
We require that the overlapping will make it possible to connect the different grids by centred interpolation formulas. The equatorial grid will be extended towards the poles by gridpoints on one or a few parallels, depending on the order of the scheme, and the polar grids will be extended towards the equator. We recall that on the meridian λ = 0, the step is $\Delta \theta $ in both grids, and further that they have one gridpoint in common, with latitude θ_{0}. We extend the equatorial grid so that it will cover the region $|\theta |\le {\theta}_{1}$, with
In Fig. 2, the E-P system, for N= 40 and with the overlapping included, is orthogonaly projected on the equatorial plane. We notice, for example, that for N= 1000 the overlapping zone would hardly be visible.
When we discretise on the grids by using one of the formulas in (8), values for the dependent variables are needed at some points outside the grid, but belonging to the stencil used. These points will be called interpolation points because the corresponding values will be determined by interpolation. In this article, we use centred bi-variate Lagrangian interpolation of order $2p+1$, where p is given by the order 2p for our centred discretisation. Points for which discretisations of the differential equations are used, will be called discretisation points.
The interpolation points for the equatorial grid are gridpoints on the parallels with latitude ${\theta}_{1}+k\Delta \theta ,k=1,\dots ,p$. The southern hemisphere is treated analogously. On the northern polar grid we proceed as follows. We approximate the circle of parallel corresponding to $\theta ={\theta}_{2}-\Delta \theta $ by a closed polygon having its corners at gridpoints of the northern polar grid, and with interior $\text{angles}\ge 3\pi /4$. This can be done in a step wise manner, in which each step consists of choosing among three appropriate grid points the one closest to the parallel. The points chosen will be interpolation points whose interior gridpoints are defined to be the discretisation points. By using the latter together with the stencil for our discretisation, the rest of the interpolation points on the northern polar grid can be determined. This procedure gives us good control over where the interpolation points of the northern polar grid will be placed in the equatorial grid, which is important when we want to minimise the overlapping. Because of symmetry with respect to the equator, no preparations concerning connections and interpolation formulas are needed for the southern polar grid.
We say that the overlapping is strict if the interpolation formulas only use discretisation points in other grids. A strict overlapping is called minimal if its overlapping number s is minimal.
Section 4.3.1 includes an example of an unstable scheme, using the E-P grid system, which violates the requirements of strict overlapping. We therefore recommend the use of strict overlapping for the kind of connections used in this article. This concept was introduced already in Starius (1980), but differently formulated and without a name.
Using our computer program we found that for our fourth and sixth order methods, minimal overlapping requires s= 2 and s= 3, respectively. Greater overlapping numbers s are of course possible. However, more than minimal overlapping can be meaningless, as we will see in Section 4.2.2.
In this section, we try to numerically evaluate the E-P grid system by considering a number of test examples, most of which are taken from Williamson et al. (1992). The spatial discretisations are given in Section 2.4, and for the time integration the classical explicit fourth order Runge–Kutta method is used, for which interpolations between grids are needed for each of the four stages. It is assume that the errors caused by each time integration below, not only ours, are small enough to be negligible in our comparisons between different methods for the SWE.
Methods for NWP must be able to use many thousands of processors in an efficient way, which favours the use of explicit time integration, since parallelisation can be considerably better and more easily implemented than for semi-implicit methods. Further, it might be desirable to resolve gravitational waves, which would make both ordinary explicit and exponential time integration methods (Clancy and Pudykiewics, 2013; Gaudreault and Pudykiewicz, 2016) more interesting.
Some further comments about the implementation. The variables U, V, and $\Psi $ on the southern and northern polar grids are represented by square matrices, corresponding to squares enclosing the polar circular regions. Gridpoints in the squares but not in the polar circular regions can be ‘enion in the computations, if this is considered to simplify the programming. By representing factors in Equation (7) such as
In all the tables below, 2p is the order of approximation of the discretisations for (2) and (7). The order of the interpolation formulas connecting the grids is $2p+1$. The overlapping number for the grid system is denoted by s, as in Section 3.3. Modified stereographic coordinates are always used, except in Table 3.
The geopotential $\Psi =gH$, where g is the acceleration of gravity and H the height of the fluid. The relative error norms ${\ell}_{2}(H)$ and ${\ell}_{\infty}(H)$, introduced in Williamson et al. (1992), will be used.
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 in (2) and (7). Since hyper-diffusion has no known physical meaning, we consider our smoothing a purely numerical process, used for carefully stabilizing the scheme. We do not aim at making our centred scheme less oscillatory by smoothing. Actually, the oscillatory behaviour decreases in the order of the centred scheme, cf. Table 2.
Let L_{h} be the five-point operator for the Cartesian form of the Laplacian, that is
To the first, the second, and the third system in Equation (2), we add the terms $-c{L}_{h}^{2}u,-c{L}_{h}^{2}v$, and $-c{L}_{h}^{2}\Psi $, respectively. In (7), U and V are used instead of u and v. We observe that ${L}_{h}^{2}\Psi =O(\Delta {\theta}^{4})$, for any smooth function $\Psi $. Note that in the discretisations the smoothing terms will be multiplied by $\Delta t$. The smoothing constant $c\ge 0$ is approximately determined by numerically minimising ${\ell}_{2}(H)$. In Section 4.2.2 ${\ell}_{\infty}(H)$ is minimised instead. We have used the MATLAB^{®} routine fminbnd.
We have also tried to let L_{h} correspond to a five-point difference approximation to the Laplacian in spherical coordinates. Not surprisingly, this did not improve the results, which seems natural since our grids are almost uniform.
We recall that the advective equation $dH/dt=0$ is linear, generally with variable coefficients, and describes conservation of mass.
In this first example two vortices will be created during the time integration and their positions can be given by a parameter, which makes it possible to place their centres in our overlapping zones. The details can be found in Nair and Jablonowski (2008). We give a short summary here, partly because we are going to use some of the formulas in this section also in Section 4.3.2.
Let $(\lambda \prime ,\theta \prime )$ be a rotated spherical coordinate system whose North pole, $(\lambda \prime ,\pi /2)$, is located at $({\lambda}_{p},{\theta}_{p})$ in another spherical system $(\lambda ,\theta )$. By Nair and Jablonowski (2008), the transformation between the two spherical coordinate systems can be written
The vortices will be created with centres at the poles of the system $(\lambda \prime ,\theta \prime )$ and will later be transformed to the system $(\lambda ,\theta )$. The normalized tangential velocity of the vortex is defined to be $V={v}_{0}\frac{3\sqrt{3}}{2}\frac{\mathrm{tanh}(\rho )}{{\mathrm{cosh}}^{2}(\rho )}$, where $\rho ={\rho}_{0}\hspace{0.17em}\text{cos}(\theta \prime )$ is the radial distance of the vortex, and ${\rho}_{0}=3$. In the initial condition (cf. Nair and Jablonowski, 2008), γ = 5, and the velocity ${v}_{0}=2\pi a/T$, where T is the total time of the simulation, set to be 12 days. For this deforming vortex case the angular velocity ω varies and is defined to be $\omega (\theta \prime )=V(\rho )/(a\rho )$, and we let ω = 0, for $\theta \prime =\pm \pi /2$.
In Nair and Jablonowski (2008), the zonal and meriodinal velocities are given, in the system $(\lambda ,\theta )$, by
Further set ${\lambda}_{p}=\pi $ and $\alpha =\pi /2-{\theta}_{p}$, where the parameter α is the opening angle between the North poles of the two spherical systems, seen from their common origin. We mention that in the system $(\lambda ,\theta )$ the case α = 0 means that the centres of the vortices are at the poles and the flow is along circles of parallel, and $\alpha =90\xb0$ means that the centres are on the equator. Since our grids are connected at latitudes $\theta =\pm \pi /4$, the value $\alpha ={45}^{o}$ is of particular interest to us, because the centres of the vortices will be located on parallels corresponding to these latitudes, cf. Fig. 3.
In Table 1 the ${\ell}_{2}(\Psi )$ errors, the ones we have minimised, are actually smallest for $\alpha =45\xb0$, which indicates that the connections between grids work quite well.
We conclude by some convergence results for $\alpha =45\xb0$ and use (2p, N, $\Delta t,{10}^{8}c,{\ell}_{2}(\Psi )$) for the corresponding computation. From (4, 240, 300, 1.1, 2.1457e − 5), (4, 360, 200, 0.7, 4.3405e − 6) we get the order 3.94(4) and (6, 240, 150, 0.080, 1.4892e − 6), (6, 360, 100, 0.030, 1.4129e − 7) give the computed order 5.81(6). The above computed orders of approximation are a bit too small, which probably depends on lack of smoothness or resolution at the centers of the vortices.
In this example, we consider a solid body rotation of the Cosine bell, according to test problem no. 1 in Williamson et al. (1992). The advecting wind is given by
In rows 4 and 5 of Table 2, we have used more overlapping, s = 3, than was needed, which led to no real increase in accuracy. This is one of the reasons we recommend minimal overlapping.
In Ronchi et al. (1996), an interesting gridding technique is introduced by combining the cubed sphere decomposition with overlapping grids. For the present problem their spatial discretisation is a centred fourth order difference scheme, used with N = 360 and $\Delta t=600$s. The main difference compared to our method, is in the choice of the grid system. For $\alpha ={45}^{o}$ corner points of the cube are inside the support of the bell during about one-third of a revolution. Unfortunately, for this case cancellations of errors occur, which distorts our comparison between the methods. We mention that the $1000\xb7{\ell}_{\infty}(H)$ error in Ronchi et al. (1996) is 7.33 and ours is 6.29.
For $\alpha =90\xb0$ the $1000\xb7{\ell}_{\infty}(H)$ error in Ronchi et al. (1996) is 10.17 and ours is 7.45. For given resolution, the error in Ronchi et al. (1996) can only be marginally decreased, by our or other grid systems, since the track of their Cosine bell is along a meridian passing through the midpoints on four edges, cf. figures in Staniforth and Thuburn (2012) and Williamson (2007), and thus displays the best uniformity possible, for the Cubed sphere grid. The mesh covering this track is almost uniform.
We conclude this example by some convergence results. Despite the solutions are only in ${C}^{1}(\mathit{\text{sphere}})$ the errors, ${\ell}_{\infty}(H)$ in Table 2, are decreasing fairly rapidly in the formal orders of approximation. As usual, lower bounds for the errors are settled by the continuous problem and the used grid and its resolution. We have minimised $c\to {\ell}_{2}(H)$, with N = 240, $\Delta t=900$s and later with N = 360, $\Delta t=600$s and got the computed orders 2.17, 2.19 and 2.33, corresponding to the formal orders 4, 6 and 8, respectively.
In the present section, we study the time integration of the complete inviscid shallow water equations.
We now consider test problem no. 2 from Williamson et al. (1992), which constitutes a steady state solution to the shallow water equations. The advecting wind is the same as for the Cosine bell. The parameter α is the angle between the North poles of the systems $(\lambda ,\theta )$ and $(\lambda \prime ,\theta \prime )$, both having the centre of the earth as their origin. For the latter system the pole $(\lambda \prime ,\pi /2)$ coincides with the North pole of the earth. The solution corresponds to solid body rotation and zonal flow in the system $(\lambda \prime ,\theta \prime )$. To express the Coriolis parameter f and the geopotential $\Psi $ in the coordinates $(\lambda ,\theta )$, the relation $\text{sin}\hspace{0.17em}\theta \prime =\hspace{0.17em}\text{sin}\hspace{0.17em}\theta \hspace{0.17em}\text{cos}\hspace{0.17em}\alpha -\hspace{0.17em}\text{cos}\hspace{0.17em}\theta \hspace{0.17em}\text{cos}\hspace{0.17em}\lambda \hspace{0.17em}\text{sin}\hspace{0.17em}\alpha $ is useful. We have $f=2\Omega \hspace{0.17em}\text{sin}\hspace{0.17em}\theta \prime $ and $\Psi ={\Psi}_{0}-(a\Omega {u}_{0}+{u}_{0}^{2}/2){\hspace{0.17em}\text{sin}\hspace{0.17em}}^{2}\theta \prime $, where ${\Psi}_{0}=2.94\xb7{10}^{4}$ and ${u}_{0}=2\pi a/12$ m/day.
Numerical results are presented in Table 3, both for stereographic and modified stereographic coordinates. For $\alpha ={45}^{o}$ there are significant differences. This case is more sensitive than the other two, with $\alpha ={0}^{o}\text{and}{90}^{o}$, for which the wind fields are along parallels and mainly along meridians, respectively. For $\alpha =45\xb0$ the directions of the wind often form an angle of about $45\xb0$ with parallels. The superiority of the modified stereographic grid, in particularly for the order 6, must be caused by the better uniformity. From now on we restrict ourselves to these coordinates, for which the errors in Table 3 are fairly the same for $\alpha =45\xb0$ and $\alpha =90\xb0$. For $\alpha =0\xb0$ the solutions are much simpler than for the other cases, which explains the higher accuracy.
In an additional experiment, not presented in a table, we replaced s= 2 with s= 1 in the third row of Table 3, which implied that a lot of points used in interpolation formulas were interpolation points from another grid, so that the requirement of strict overlapping was violated. The scheme was unstable ($\mathrm{min}\Psi <8\xb7{10}^{3}$ before 10 days), up to about c= 70 and then stabilized. The minimal ${\ell}_{2}$ error, occurring for $c\sim 1000$, was about 900 times larger than for the case s= 2 in Table 3. This shows that strict overlapping should be used for the kind of connections considered here.
In Ullrich et al. (2010), a fourth order finite volume method on a cubed sphere grid, having 25 points in its stencil, is used for $\alpha =45\xb0$, and with N = 160 and N= 320 gridpoints on the equator. We compare it with our 6th order method, having only 13 points in its stencil. The ${\ell}_{2}$ errors in Ullrich et al. (2010) are about 90 and 520 times larger than ours for N= 160 and N= 320, respectively. We could have used the 4th order difference method instead, with only 9 points in its stencil, but the corresponding comparisons would be too unbalanced, since the amount of computational work per gridpoint is also important, not only the order of approximation. Unfortunately, the above is not a clear-cut comparison between grid systems.
We now consider some convergence results, for $\alpha =45\xb0$, and use the same notation as in our first example. From (4, 240, 300, 15.0, 1.3341e − 7), (4, 360, 200, 12.5, 2.5275e − 08) we get the computed order 4.10(4). Similarly from (6, 240, 200, 0, 1.0253e − 10), (6, 360, 200, 0, 9.0908e − 12) we arrive at the computed order 5.98(6).
In Fig. 4, we have visualised the error ${\Psi}_{\mathit{\text{computed}}}-\Psi $ as function of latitude and longitude. The true solution $\Psi $ has two minimum points per period, namely $\theta =45\xb0,\lambda =180\xb0$ and $\theta =-45\xb0,\lambda =0\xb0$ or $360\xb0$. In the vicinities of these points the curvature is much greater than elsewhere, which leads to much larger errors in these vicinities. Everything in the figure can easily be explained by the solution $\Psi $, thus there is no harmful grid imprinting.
We now continue with test problem no. 3 from Williamson et al. (1992). The spherical coordinate systems $(\lambda \prime ,\theta \prime )$ and $(\lambda ,\theta )$ are the same as in Section 4.2.1 and the problem is formulated in the former system. For the velocity components $(u\prime ,v\prime )$, in the system $(\lambda \prime ,\theta \prime )$ we have $v\prime =0$ and $u\prime $ depends only on $\theta \prime $ and has compact support (cf. Williamson et al., 1992). The geopotential $\Psi $ is determined by an ordinary differential equation. Transformation to the rotated system $(\lambda ,\theta )$ leads to the same expression for the Coriolis parameter as in Section 4.3.1 and initial conditions were obtained by using Equations (13) and (14). In the last set of formulas $a\omega (\theta \prime )$ should be replaced by $u\prime /\hspace{0.17em}\text{cos}\hspace{0.17em}\theta \prime $. The results are presented in Fig. 5 and Table 4.
In Browning et al. (1989), two overlapping stereographic grids are used together with a centred difference method of order 6. The resolution compares well to our N= 240. Only the case α = 0 was considered and the ${\ell}_{2}$ error was $6.6\xb7{10}^{-5}$, which is about $4.9\times {10}^{3}$ times larger than our error. The grid system in Browning et al. (1989) is different from ours, with much less uniformity. Nevertheless, the difference in accuracy seems surprisingly high.
The order of approximation will now be estimated experimentally, according to Williamson et al. (1992), for $\alpha =60\xb0$ and for our method of formal order 4. The resolutions are determined by N = 160, $N=1.5\times 160=240$, and $N=1.5\times 240=360$, all with $\Delta t=150$s. With smoothing parameter ${10}^{8}c=10$, in all three cases, we find the ${\ell}_{2}$ errors $2.4309\times {10}^{-5},4.0756\times {10}^{-6}$, and $7.8697\times {10}^{-7}$. By using the first two and the last two errors, we found the approximate orders 4.404 and 4.056, respectively. By optimising each case, we found the ${10}^{8}c$ values 21, 11.5, and 8.5, and the ${\ell}_{2}$ errors $2.1458\times {10}^{-5},4.0715\times {10}^{-6}$, and $7.8630\times {10}^{-7}$, and the corresponding orders 4.099 and 4.056.
The general features of Fig. 6 can be explained by using Fig. 5 and comments similar to those made for Fig. 4. We conclude that there is no harmful grid imprinting in this case either.
In Fig. 7, we give time traces of normalised errors, for the two geostrophic flow problems in Williamson et al. (1992). The ${\ell}_{2}(\Psi )$ errors, the ones that were minimised over c at the endpoint five days, grow linearly in time with no visible deviation.
We will use test problem no. 6 from Williamson et al. (1992) in order to study the global conservation properties for mass and energy. This test example was already used in Phillips (1959) and with a numerical procedure in some aspects related to our E-P based method. The integral formulas for mass and energy can be found in, for example Starius (2014). All dependent variables were interpolated to a full latitude–longitude grid before the numerical integration. The sum of squares of the relative errors for the total mass after each day, $\left|\right|re{l}_{\mathit{\text{mass}}}(1:14)|{|}_{2}$, was minimised over c.
In Table 5, we have used exactly the same numerical method as described earlier, without any corrections or fixers of any kind. The schemes does not explicitly conserve mass, but nevertheless the accuracy in the total mass over time is considerable, probably because of cancellation of errors. This situation is not new to us (cf. Starius, 1977a, 1977b). About seven significant decimal digits in the total mass is obviously sufficient for NWP. The relative changes in the total energy seem to be satisfactory, but are much larger than for the total mass. This is probably because, for non-linear problems, energy can move between different wave numbers.
We have introduced a new grid system for the sphere, called E-P, consisting of one equatorial and two polar grids. The former is an ordinary spherical latitude–longitude grid and the latter ones correspond to suitably modified stereographic coordinates. All grids are orthogonal and almost uniform and are connected by minimal overlapping and centred interpolation formulas, for dependent variables.
The superiority of the modified stereographic coordinates compared to the standard ones is most apparent in Table 3, for the two cases with $\alpha =45\xb0$.
The uniformity is optimised in Section 3 by using so called deviation factors, cf. (9) and (12), for the whole grid system. The comparisons made between different grid systems, in Sections 4.2.2, 4.3.1 and 4.3.2, clearly indicate the usefulness of uniformity, which would probably be more pronounced for more typical wave propagation problems. In a forthcoming paper, we plan to return to this point in connection with a comparison between the well-known overset Yin–Yang and our E-P grid systems. Also, a generalisation of the E-P system will be considered by letting the equatorial grid be a suitable reduced grid (cf. Gates and Riegel, 1963; Kurihara, 1965). This will improve the uniformity and increase the area covered by the equatorial grid considerably, and decrease the total number of gridpoints, and be important for wave propagation problems. Further, the circular polar grids, used in this article, will be compared to the use of square polar grids. They are, at the least, somewhat easier to implement, but there might be other differences as well.
In this short article, all focus is on our overset E-P grid system, which can be looked upon as a way to concentrate the unavoidable non-uniformity to narrow bands, where it is handled by overlapping and accurate centred interpolation formulas.
No potential conflict of interest was reported by the author.
Bao, L., Nair, R. D. and Tufo, H. M. 2014. A mass and momentum flux-form higher-order discontinuous Galerkin shallow water model on the cubed-sphere. J. Comput. Phys. 271, 224–243. doi:https://doi.org/10.1016/j.jcp.2013.11.033.
Browning, B. L., Hack, J. J. and Swarztrauber, P. N. 1989. A comparison of three numerical methods for solving differential equations on the sphere. Mon. Wea. Rev. 117, 1058–1075. doi:https://doi.org/10.1175/1520-0493(1989)117<1058:ACOTNM>2.0.CO;2.
Chen, C., Li, X., Shen, X. and Xiao, F. 2014. Global shallow water models based on multi-moment constrained finite volume method and three quasi-uniform spherical grids. J. Comput. Phys. 271, 191–223. doi:https://doi.org/10.1016/j.jcp.2013.10.026.
Clancy, C. and Pudykiewics, J. A. 2013. On the use of exponential time integration methods in atmospheric models. Tellus A65, 20898. doi:https://doi.org/10.3402/tellusa.v65i0.20898.
Gates, W. L. and Riegel, A. 1963. A study of numerical errors in the integration of barotropic flow on a spherical grid. JGR67, 406–423.
Gaudreault, S. and Pudykiewicz, J. A. 2016. An efficient exponential time integration method for the numerical solution if the shallow water equations on the sphere. J. Comput. Phys. 322, 827–848. doi:https://doi.org/10.1016/j.jcp.2016.07.012.
Giraldo, F. X. 2006. High-order triangle-based discontinuous Galerkin methods for hyperbolic equations on a rotating sphere. J. Comput. Phys. 214, 447–465. doi:https://doi.org/10.1016/j.jcp.2005.09.029.
Henshaw, W. D. and Schwendeman, D. W. 2006. Moving overlapping grids with adaptive mesh refinement for high-speed reactive and non-reactive flow. J. Comput. Phys. 216, 744–779. doi:https://doi.org/10.1016/j.jcp.2006.01.005.
Ii, S. and Xiao, F. 2010. A global shallow water model using high order constrained finite volume method and icosahedral grid. J. Comput. Phys. 229, 1774–1796. doi:https://doi.org/10.1016/j.jcp.2009.11.008.
Kageyama, A. and Sato, T. 2004. The Yin-Yang grid: An overset grid in spherical geometry. Geochem. Geophys. Geosyst. 5, Q09005.
Kurihara, Y. 1965. Integration of the primitive equations on a spherical grid. Mon. Wea. Rev. 93, 399–415. doi:https://doi.org/10.1175/1520-0493(1965)093<0399:NIOTPE>2.3.CO;2.
Lanser, D., Blom, J. G. and Verwer, J. G. 2000. Spatial discretisation of the shallow water equations in spherical geometry using Osher’s scheme. J. Comput. Phys. 165, 542–565. doi:https://doi.org/10.1006/jcph.2000.6632.
Lauritzen, P. H., Nair, R. D. and Ullrich, P. A. 2010. A conservative semi-Lagrangian multi-tracer scheme (CSLAM) on a cubed-sphere grid. J. Comput. Phys. 229, 1401–1424. doi:https://doi.org/10.1016/j.jcp.2009.10.036.
Läuter, M., Giraldo, F. X., Handorf, D. and Dethloff, K. 2008. A discontinuous Galerkin method for the shallow water equations using spherical triangular coordinates. J. Comput. Phys. 227, 10226–10243. doi:https://doi.org/10.1016/j.jcp.2008.08.019.
Nair, R. D. and Jablonowski, C. 2008. Moving vortices on the sphere: A test case for horizontal advection problems. Mon. Wea. Rev. 136, 699–711. doi:https://doi.org/10.1175/2007MWR2105.1.
Nair, R. D., Thomas, S. J. and Loft, R. D. 2005. A discontinuous Galerkin transport scheme on the cubed sphere. Mon. Wea. Rev. 133, 814–828. doi:https://doi.org/10.1175/MWR2890.1.
Peixoto, P. S. and Barros, S. R. M. 2013. Analysis of grid imprinting on geodesic spherical icosahedral grids. J. Comput. Phys. 237, 61–78. doi:https://doi.org/10.1016/j.jcp.2012.11.041.
Phillips, N. A. 1959. Numerical integration of the primitive equations on the hemisphere. Mon. Wea. Rev. 87, 333–345. doi:https://doi.org/10.1175/1520-0493(1959)087<0333:NIOTPE>2.0.CO;2.
Pudykiewics, J. A. 2011. On numerical solution of the shallow water equations with chemical reactions on icoshedral geodesic grid. J. Comput. Phys. 230, 1956–1991. doi:https://doi.org/10.1016/j.jcp.2010.11.045.
Putman, W. M. and Lin, S.-J. 2007. Finite-volume transport on various cubed-sphere grids. J. Comput. Phys. 227, 55–78. doi:https://doi.org/10.1016/j.jcp.2007.07.022.
Qaddouri, A. and Lee, V. 2011. The Canadian global environmental multiscale model on the Yin-Yang grid system. Qjr. Meteorol. Soc. 137, 1913–1926. doi:https://doi.org/10.1002/qj.873.
Robert, A. 1982. A semi-Lagrangian and semi-implicit numerical integration scheme for the primitive meteorological equations. J. Meteorol. Soc. Jpn. 60, 319–325. doi:https://doi.org/10.2151/jmsj1965.60.1_319.
Ronchi, C., Iacono, R. and Paolucci, P. S. 1996. The cubed sphere: A new method for the solution of partial differential equations in spherical geometry. J. Comput. Phys. 124, 93–114. doi:https://doi.org/10.1006/jcph.1996.0047.
Sadourny, R. 1972. Conservative finite-differencing of the primitive equations on quasi-uniform spherical grids. Mon. Wea. Rev. 100, 136–144. doi:https://doi.org/10.1175/1520-0493(1972)100<0136:CFAOTP>2.3.CO;2.
Staniforth, A. and Thuburn, J. 2012. Horizontal grids for global weather and climate prediction models: a review. Qjr. Meteorol. Soc. 138, 1–26. doi:https://doi.org/10.1002/qj.958.
Starius, G. 1977a. Composite mesh difference methods for elliptic boundary value problems. Numer. Math. 28, 243–258. doi:https://doi.org/10.1007/BF01394455.
Starius, G. 1977b. Constructing orthogonal curvilinear meshes by solving initial value problems. Numer. Math. 28, 25–48. doi:https://doi.org/10.1007/BF01403855.
Starius, G. 1980. On composite mesh difference methods for hyperbolic differential equations. Numer. Math. 35, 241–255. doi:https://doi.org/10.1007/BF01396411.
Starius, G. 2014. A solution to the pole problem for the shallow water equations on a sphere. TWMS JPAM. 5, 152–170.
Ullrich, P. A., Jablonowski, C. and van Leer, B. 2010. High-order finite-volume methods for the shallow-water equations on the sphere. J. Comput. Phys. 229, 6104–6134. doi:https://doi.org/10.1016/j.jcp.2010.04.044.
Volkov, E. A. 1968. The method of composite meshes for finite and infinite regions with piecewise smooth boundaries. Proc. Steklov Inst. Math. 96, 145–185.
Williamson, D. L. 1968. Integration of the barotropic vorticity equation on a spherical geodesic grid. Tellus. 20, 642–653.
Williamson, D. L. 2007. The evolution of dynamical cores for global atmospheric model. JMSJ. 85B, 241–269. doi:https://doi.org/10.2151/jmsj.85B.241.
Williamson, D. L., Drake, J. B., Hack, J. J., Jakob, R. and Swarztrauber, P. N. 1992. A standard test set for numerical approximation to the shallow water equations in spherical geometry. J. Comput. Phys. 102, 211–224. doi:https://doi.org/10.1016/S0021-9991(05)80016-6.