2.1. 

Spherical coordinates and the shallow water equations

The transformation between Cartesian and spherical coordinates on the sphere is

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 $\Psi$ is the geopotential, $f=2\Omega \hspace{0.17em}\text{sin}\hspace{0.17em}\theta$ the Coriolis parameter, Ω the rotation rate of the sphere, and u and v the eastward and northward velocity components, respectively.

2.2. 

Standard and modified stereographic coordinates

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

where $q(\theta )=2\hspace{0.17em}\text{tan}\hspace{0.17em}(\phi /2),$ with the colatitude $\phi$ equal to $\frac{1}{2}\pi -\theta$ and $\frac{1}{2}\pi +\theta ,$ on the northern and southern hemisphere, respectively. We also need the map factor $m_f(\theta )=2/(1+\alpha \hspace{0.17em}\text{sin}\hspace{0.17em}\theta ),$ where α = 1 on the northern and −1 on the southern hemisphere.

Since the velocities in spherical coordinates, $u=a \hspace{0.17em}\text{cos}\hspace{0.17em}\theta d\lambda /dt,$$v=a d\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 use of the minus sign in the definition of U, intended to get nicer formulas, namely

on the northern and southern polar regions, respectively. In (4), U and V have been scaled by the factor $1/m_f(\theta ).$ 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 $\lambda =\frac{1}{2}\pi ,\pi ,\frac{3}{2}\pi ,$ 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\prime ,y\prime )$ related to $(x_{st},y_{st})$ by

The modified coordinates $(x\prime ,y\prime )$ 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).

3.1. 

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 ${\sigma }_K=\Delta {\lambda }_1/\Delta \theta ,$ where $\Delta \theta$ will be the same for all discretised segments and determined later. We consider the segment $S_i^d$ and require that the distances between grid points on the parallel ${\theta }_{i-1}$ are equal to $\Delta {\lambda }_1,$ that is $\Delta {\lambda }_i\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{i-1}=\Delta {\lambda }_1,$ where $i=2,\dots ,K.$ Further, the maximal deviation factor for the segment $S_i^d$ is $\max (\Delta {\lambda }_i\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{i-1}/\Delta \theta , \Delta \theta /(\Delta {\lambda }_i\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_i) ),$ which is obviously minimised, with respect to $\Delta {\lambda }_i/\Delta \theta ,$ when the two deviation factors in the $\max$ above are equal. We now have the equalities

and substituting $\Delta {\lambda }_i$ in the second equality by $\Delta {\lambda }_1/\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{i-1}$ and using ${\sigma }_K=\Delta {\lambda }_1/\Delta \theta ,$ gives

The transition from $S_K^d$ to a polar grid gives that ${\sigma }_K=m_f({\theta }_K)=2/(1+\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_K),$ where $m_f(\theta )$ 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 ${\sigma }_K=2/(1+\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_K)$ and by using $\text{cos}\hspace{0.17em}{\theta }_i=\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{i-1}/{\sigma }_K^2={\sigma }_K^{-2i}$ we find θ_i, $i=1,2,\dots ,K$−1. Note that (8) implies that ${\theta }_K\to \frac{1}{2}\pi ,$ and hence ${\sigma }_K\to 1,$ as $K\to \infty .$

Now $\Delta \theta =\Delta {\lambda }_1/{\sigma }_K$ can be evaluated and for practical reasons modified by defining the integer m to be $\frac{1}{2}\pi /(\Delta {\lambda }_1/{\sigma }_K)$ rounded to the nearest integer and then redefining $\Delta \theta =\frac{1}{2}\pi /m.$ With $\Delta \theta$ known we now round the latitudes θ_i, $i=1,2,\dots ,K,$ to their closest multiples of $\Delta \theta .$ Finally we set $n_i=n\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{i-1},$ round it to the nearest integer, and define $\Delta {\lambda }_i=2\pi /n_i,i=2,\dots ,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 $s_{10}=9,$ and n = 2000, K = 10, which gives p = 20.8% and $s_{10}=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 $n_1+2(n_2+\dots +n_{K-1})+n_K,$ 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 ${\theta }_{10}={65.8}^o$) 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 ${\theta }_{K+1},{\theta }_{K+2},\dots$ can be determined by the equation $\text{cos}\hspace{0.17em}{\theta }_i={\sigma }_K^{-2i},i>K,$ and then also $n_i=n\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{i-1} \text{and} \Delta {\lambda }_i=2\pi /n_i,i>K,$ where the n_i have been rounded in a suitable way.

We conclude this subsection with two examples. For initial K = 3, we have ${\sigma }_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 ${67}^o.$ The maximal deviation factor for the whole reduced grid system is still σ₃. In a second example, with initial K = 1 and ${\sigma }_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.

3.2. 

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 $S_i^d$ with q−1 intervals in the following segment $S_{i+1}^d.$ As for method I we require that $\Delta {\lambda }_i\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{i-1}=\Delta {\lambda }_1,$ where $i=2,\dots ,K,$ and further that each n_i, with i < K, is divisible by q. This implies that $\Delta {\lambda }_{i+1}/\Delta {\lambda }_i=n_i/n_{i+1}=q/(q-1).$ Further, $\Delta {\lambda }_i\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{i-1}=\Delta {\lambda }_{i+1}\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_i$ together with (7) leads to $\Delta {\lambda }_{i+1}/\Delta {\lambda }_i=\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_{i-1}/\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_i={\sigma }^2,$ where we write σ instead of σ_K, since σ is now independent of K. We choose q and determine σ and ${\vartheta }_0,$ where ${\vartheta }_0$ is the width of the reduced grid, by using

We now determine θ_i by $\text{cos}\hspace{0.17em}{\theta }_i={\sigma }^{-2i},i=1,2,\dots$ until ${\theta }_i>{\vartheta }_0.$ Further, we can choose K from $|{\theta }_K-{\vartheta }_0|={\min }_i|{\theta }_i-{\vartheta }_0|$ and then redefine ${\vartheta }_0={\theta }_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 $S_K^d.$ Since $n_1=n$ is known we can use $n_{i+1}=(q-1)n_i/q,\text{where}i=1,\dots ,K$-1, and then define $\Delta {\lambda }_i=2\pi /n_i,i=2,\dots ,K.$

The above works only if n is divisible by $q^{K-1},$ 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 $n_{i_0},$ with not divisible by q, then $\lambda =2\pi$ will not correspond to grid points in $S_{i_0+1}^d.$ 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 $S_{i_0+1}^d$ onward, we define $n_{i+1}=(q-1)n_i/q,$ for $i=i_0,\dots ,K$ and then round them to integers and define $\Delta {\lambda }_{i+1}=2\pi /n_{i+1},$ for $i=i_0,\dots ,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 $n_{i_0+1}$ is rounded to an integer divisible by q, then the interpolation between $S_{i_0+1}^d$ and $S_{i_0+2}^d$ can be done with the q−1 interpolation formulas. We still use method I for the interpolation between $S_{i_0}^d$ and $S_{i_0+1}^d.$ This can be repeated by choosing $n_{i_0+3}$ divisible by q, using method I between $S_{i_0+2}^d$ and $S_{i_0+3}^d,$ and then the q−1 interpolation formulas between $S_{i_0+3}^d$ and $S_{i_0+4}^d,$ etc.

We conclude this subsection by considering an example with q = 6, which preliminarily gives $\sigma =\sqrt{6/5}=1.0954.$ We have ${\vartheta }_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 ${\vartheta }_0={\theta }_3,$ and obtain a new $\sigma =2/(1+\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_3)=1.1019.$ If n is divisible by 6², then method II works perfectly well for K = 3. However, if we want K > 3 and n is not divisible by 6³, we can use a modification as described above for the segments following $S_3^d.$

3.3. 

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 S^d be a segment on the northern hemisphere bounded by the latitudes and consisting of a lat-lon net with the longitudinal step $\Delta \lambda .$ We define σ as the deviation factor for net rectangles at θ₁, that is $\sigma =\Delta \lambda \hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_1/\Delta \theta ,$ where $\Delta \theta$ is the latitudinal step, not yet determined. The ratio of maximum to minimum grid length on S^d is

from which it follows

The minimum of γ is $\text{cos}\hspace{0.17em}{\theta }_1/\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_2,$ which is attained for all $\sigma \le \hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_1/\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_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 $\sqrt{\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_1/\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_2},$ which is attained if and only if $\Delta \lambda /\Delta \theta =1/\sqrt{\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_1\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_2}.$ We get square net rectangles at the latitude θ for $\text{cos}\hspace{0.17em}\theta \Delta \lambda =\Delta \theta ,$ 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 $\left|\theta \right|\le {\vartheta }_0,$ which will be determined later. Our northern polar grid will go down to ${\vartheta }_0$ and the corresponding map factor is $m_f=2/(1+\hspace{0.17em}\text{sin}\hspace{0.17em}{\vartheta }_0).$ According to (12) in (S-E-P) the greatest and smallest grid lengths in our polar grids are $\Delta \theta$ and $\Delta \theta /m_f=\Delta \theta h({\vartheta }_0),$ respectively, where $h({\vartheta }_0)=(1+\hspace{0.17em}\text{sin}\hspace{0.17em}{\vartheta }_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

where $\sigma =\Delta \lambda /\Delta \theta .$ We now consider the equation $h({\vartheta }_0)/(\sigma \hspace{0.17em}\text{cos}\hspace{0.17em}{\vartheta }_0)=1,$ with a unique root ${\theta }^{*}(\sigma )$ in the interval $(0\pi /2).$ Elementary computations show that

We look at two cases, namely

In both cases $\gamma (\sigma ,{\vartheta }_0)\ge 1/\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }^{*}(1)=1.25,$ with equality if and only if σ = 1 and ${\vartheta }_0={\theta }^{*}(1)\approx {36.87}^o.$ 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.

3.4. 

Square net rectangles at mid-latitudes

Grids for NWP are generally constructed to be uniform (having square net rectangles) at mid-latitudes, $\theta =\pm \pi /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 $S_{i_0}^d$ with and then use (9) with $\theta =\pi /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 $\text{cos}\hspace{0.17em}{\overline{\theta }}_{i_0}={\overline{\sigma }}^{-2i_0},$ which together with (12) gives ${\overline{\sigma }}^{-4i_0+2}=1/2.$ We then determine $\overline{\sigma }$ and ${\overline{\theta }}_i,$$i=1,\dots ,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 $\Delta \lambda ={0.1}^o$ 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 $\theta =\pm {45}^o.$ If the Yang grid is congruent with the Yin grid the total number of grid points will be $(2·\sqrt{2}·3/16)n^2,$ which for n = 3600 is about $6.873·{10}^6.$ 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·{10}^6.$

4.1. 

Smoothing

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 L_h be the five-point operator associated with minus the Laplacian $-{\nabla }^2,$ and defined by

To the first, second, and third equation in (2), we add the terms $-c{L}_h^{p}u,-c{L}_h^{p}v,$ and $-c{L}_h^p\Psi ,$ 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 L_h for the whole sphere. We observe that $L_h^p\Psi =O(\Delta {\theta }^{2p}),$ for any smooth function $\Psi .$ In the temporal discretisations the smoothing terms will be multiplied by $\Delta t.$ The smoothing constant $c\ge 0$ will often be determined by minimising ${\mathcal{\ell }}_2(H,T),$${\overline{\mathcal{\ell }}}_2(H,T),{\overline{r}}_m(T)$ or ${\overline{r}}_e(T).$ We have used the MATLAB^® routine fminbnd, for the numerical minimisation.

In an earlier work we have used L_h 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.

4.2. 

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 $u_0=2\pi 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 ${\mathcal{\ell }}_2(H,12),$${\overline{r}}_m(12)$ and ${\overline{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 ${\mathcal{\ell }}_2(H,12)$ leads to slightly more accurate values of H, for K = 1 than for K = 3, but the values of ${\overline{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 ${\overline{r}}_m(12)$ leads to much greater values of c than for the other alternatives, which results in greater values for ${\mathcal{\ell }}_2(H,12)$ and ${\overline{r}}_e(12),$ see for example the case K = 1. For this reason we reject the use of ${\overline{r}}_m(12)$ for the determination of c and note that ${\overline{r}}_e(12)$ works better for this example if we focus on the other errors ${\mathcal{\ell }}_2(H,12)$ and ${\overline{r}}_m(12)$ and not on the values of c themselves, which are much smaller compared to the ones corresponding to ${\mathcal{\ell }}_2(H,12).$

4.3. 

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 $u_0=2\pi a/12$ m/day, $k_1=133681$ and $k_2=10m^2/s^2.$ We will compute ${\Psi }^{*}=g{H}^{*},$ 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 $\Delta 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 $\Delta t=300$ than for $\Delta t=200.$ For $2p=6$ the accuracy obtained for K = 1,2,3,4 is very much the same, particularly for $\Delta 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.1010e^{-9},4.1408e^{-9},6.2201e^{-9}),$ 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, ${\overline{r}}_e(15)$ is minimised in order to investigate whether it may be used instead of ${\overline{\mathcal{\ell }}}_2({\Psi }^{*},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 ${\overline{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.

4.4. 

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 $({\lambda }_c,{\theta }_c)$ and the initial conditions are

where $u_0=20$ m/s and $h_0=5960$ m. The geopotential corresponding to the conical mountain is and $r=\min (\sqrt{{(\lambda -{\lambda }_c)}^2+{(\theta -{\theta }_c)}^2},R).$ We have computed $\Psi ={\Psi }^{*}+{\Psi }_s,$ which occurs in the two momentum equations, whereas ${\Psi }^{*}=\Psi -{\Psi }_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.

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

where c_min guarantees that there is enough smoothing for stability and c_max 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 c_max can be infinite. For the present example we can choose $c_{\mathit{min}}=0.1\times {10}^{-6}$ and $c_{\mathit{max}}=18\times {10}^{-6}$ for $2p=6,$ and $c_{\mathit{min}}=0.5\times {10}^{-6}$ and $c_{\mathit{max}}=140\times {10}^{-6}$ 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.