Articles | Volume 15, issue 6
https://doi.org/10.5194/gmd-15-2561-2022
https://doi.org/10.5194/gmd-15-2561-2022
Development and technical paper
 | 
28 Mar 2022
Development and technical paper |  | 28 Mar 2022

Improved double Fourier series on a sphere and its application to a semi-implicit semi-Lagrangian shallow-water model

Hiromasa Yoshimura
Abstract

One way to reduce the computational cost of a spectral model using spherical harmonics (SH) is to use double Fourier series (DFS) instead of SH. The transform method using SH usually requires O(N3) operations, where N is the truncation wavenumber, and the computational cost significantly increases at high resolution. On the other hand, the method using DFS requires only O(N2log N) operations. This paper proposes a new DFS method that improves the numerical stability of the model compared with the conventional DFS methods by adopting the following two improvements: a new expansion method that employs the least-squares method (or the Galerkin method) to calculate the expansion coefficients in order to minimize the error caused by wavenumber truncation, and new basis functions that satisfy the continuity of both scalar and vector variables at the poles. Partial differential equations such as the Poisson equation and the Helmholtz equation are solved by using the Galerkin method. In the semi-implicit semi-Lagrangian shallow-water model using the new DFS method, the Williamson test cases and the Galewsky test case give stable results without the appearance of high-wavenumber noise near the poles, even without horizontal diffusion and without a zonal Fourier filter. In the Eulerian advection model using the new DFS method, the Williamson test case 1, which simulates a cosine bell advection, also gives stable results without horizontal diffusion but with a zonal Fourier filter. The shallow-water model using the new DFS method is faster than that using SH, especially at high resolutions and gives almost the same results, except that very small oscillations near the truncation wavenumber in the kinetic energy spectrum appear only in the shallow-water model using SH.

1 Introduction

Global spectral atmospheric models using the spectral transform method with spherical harmonics (SH) as basis functions are widely used. They are used in the Japan Meteorological Agency (JMA, 2019) and the Meteorological Research Institute (MRI; Yukimoto et al., 2011, 2019) for a range of applications, including operational weather prediction, operational seasonal prediction, and global warming projection. The spectral model has the advantage that the horizontal derivatives are accurate, and the semi-implicit scheme, which improves numerical stability, can be easily applied because the Helmholtz equation and the Poisson equation are easily solved in spectral space. The application of the semi-implicit semi-Lagrangian scheme allows for time steps longer than the Courant–Friedrichs–Lewy (CFL) condition, which makes the model computationally efficient. In the spectral model using SH, the Legendre transform used in the latitudinal direction significantly increases the computational cost at high resolutions since the Legendre transform usually requires O(N3) operations and O(N3) memory usage (unless using the fast Legendre transform or on-the-fly computation of the associated Legendre functions shown below), where N is the truncation wavenumber. One way to reduce the operation count and the memory usage at high resolutions with large N is to use the fast Legendre transform (Suda, 2005; Tygert, 2008; Wedi et al., 2013; Wedi, 2014), which requires only O(N2(log N)3) operations and also effectively reduces the memory usage. In the fast Legendre transform, the threshold parameter affecting the accuracy–cost balance is chosen so that a loss of accuracy is sufficiently small. Dueben et al. (2020) presented global simulations of the atmosphere at 1.45 km grid spacing in the SH model using the fast Legendre transform. Another approach to improve the Legendre transform is on-the-fly computation of the associated Legendre functions (Schaeffer, 2013; Ishioka, 2018), which still requires O(N3) operations but only O(N2) memory usage. This low memory usage also contributes to speeding up calculations by taking advantage of the cache memory.

Alternatively, we can use double Fourier series (DFS) as basis functions to reduce the operation count and the memory usage in the global spectral model. In the DFS model, the fast Fourier transform (FFT; Cooley and Tukey, 1965; Swarztrauber, 1982) is used in not only the longitudinal (zonal) direction but also the latitudinal (meridional) direction. The FFT requires only O(N2log N) operations and O(N) or O(N2) memory usage, and it is faster than the fast Legendre transform.

In DFS models (and also in SH models), the scalar variable F(λ,θ) is zonally expanded as

(1) F λ , θ m = - M M F m θ e i m λ ,

where λ is longitude, θ is colatitude, and M is the zonal truncation wavenumber. Several methods have been proposed for meridional expansion with DFS. Merilees (1973b), Boer and Steinberg (1975), and Spotz et al. (1998) performed the Fourier transform meridionally along a great circle. Spotz et al. (1998) showed that by using the spherical harmonic filter, the explicit DFS shallow-water model using the pseudo-spectral method can produce results comparable to the SH model in terms of accuracy and stability. However, the spherical harmonic filter consists of the forward SH transform (from grid space to spectral space) followed by the inverse SH transform (from spectral space to grid space), which increases the computational cost.

Orszag (1974) and Boyd (1978) expanded Fm(θ) meridionally as

(2a)Fmθfmθfor evenm,sinθfmθfor oddm,(2b)fmθn=0Nfn,mcosnθ,

where N is the meridional truncation wavenumber. The coefficients fn,m for odd m are calculated by the forward Fourier cosine transform of Fmθ/sinθ. Orszag (1974) imposed the following conditions at the poles:

(3) f m 0 = 0  and  f m π = 0  for  m 2 ,

which can be expressed in terms of the expansion coefficients fn,m as

(4) n = 0 n  is even N f n , m = 0 and n = 1 n  is  odd N f n , m = 0  for  m 2 .

Satisfying the above conditions ensures that the scalar variable F(λ,θ) and its gradient F are continuous at the poles. In Orszag (1974), only fN-1,m and fN,m were modified to satisfy Eq. (4), but this is not the best way to satisfy the same conditions as Eqs. (3) or (4), as will be shown in Sect. 4.

Yee (1981), Akahori et al. (2001), and Layton and Spotz (2003) expanded Fm(θ) as

(5) F m θ = n = 0 N F n , m cos n θ  for even m , n = 1 N F n , m sin n θ for odd m .

In the semi-implicit semi-Lagrangian shallow-water model in Layton and Spotz (2003), the spherical harmonic filter was applied to the prognostic variables for stability and accuracy. Layton and Spotz (2003) explained that the expansion with Eq. (5) permits discontinuity at the poles and non-isotropic waves, which may lead to a prohibitive time step restriction and numerical instability, and these problems can be avoided by applying the spherical harmonic filter.

Cheong (2000a, b) proposed expanding Fm(θ) as

(6) F m θ n = 0 N F n , m cos n θ for m = 0 , n = 1 N F n , m sin n θ for  odd m , n = 1 N F n , m sin θ sin n θ for even m 0 .

The meridional basis functions sin θsin nθ for even m (≠0) are different from Eq. (5). The coefficients Fn,m for  even m0 are calculated by the forward Fourier sine transform of Fmθ/sinθ. The basis functions in Eq. (6) automatically satisfy the same conditions at the poles as Eq. (3) for even m and guarantee the continuity of the scalar variable F at the poles, which is an advantage compared with the basis functions in Eq. (5). On the other hand, Eq. (6) does not automatically satisfy the conditions in Eq. (3) for odd m and does not guarantee the continuity of F at the poles. The shallow-water model and the vorticity equation model using a semi-implicit Eulerian scheme ran stably without the spherical harmonic filter by using high-order horizontal diffusion with O(N2) operations to smooth out the high-wavenumber components (Cheong, 2000b; Cheong et al., 2002; Kwon et al., 2004). The semi-implicit Eulerian hydrostatic atmospheric model also ran stably with high-order horizontal diffusion (Cheong, 2006; Koo and Hong, 2013; Park et al., 2013). However, the computational results of these models appear to be a little different from (i.e., slightly worse than) the models using SH. One reason for this seems to be the appearance of high-wavenumber oscillation resulting from the meridional wavenumber truncation with N=N2J/3 or J/2 for even m (≠0) (see Sect. 4) and the use of high-order horizontal diffusion to smooth out the oscillation, where J is the number of grid points in the latitudinal direction.

Yoshimura and Matsumura (2005) and Yoshimura (2012) stably ran the two-time-level semi-implicit semi-Lagrangian hydrostatic and nonhydrostatic atmospheric models using the DFS basis functions of Cheong in Eq. (6). These models used the same fourth-order horizontal diffusion as the SH models and did not require the spherical harmonics filter or the strong high-order horizontal diffusion for stability. The numerical stability of the models was improved by adopting the following methods.

  1. The semi-Lagrangian scheme is used, which avoids the numerical instability due to the nonlinear advection term.

  2. The meridional truncation with N=J-1 and N=J is used, which enables the accurate reconstruction of the given grid data with the expansion coefficients (Cheong et al., 2004) and avoids the error due to the meridional truncation.

  3. U=usin θ and V=vsin θ instead of u/sinθ and v/sinθ are transformed from grid space to spectral space, where u is the zonal wind and v is the meridional wind.

The results of these models were very similar to those of the SH models. However, we found the following two problems in these models.

  1. High-wavenumber noise appears near the poles.

  2. The meridional truncation wavenumber N needs to be equal to J for even m (≠0) because N<J (e.g., N2J/3) for even m (≠0) causes the high-wavenumber oscillation (see Sect. 4) and numerical instability.

To solve these problems, we propose a new DFS method that adopts the following two improvements.

  1. A new expansion method to calculate DFS expansion coefficients of scalar and vector variables, which adopts the least-squares method (or the Galerkin method) to minimize the error due to the meridional wavenumber truncation, is used.

  2. New DFS basis functions that automatically satisfy the pole conditions in Eq. (3) are introduced, which guarantees continuity of not only scalar variables but also vector variables at the poles.

We also use the Galerkin method to solve partial differential equations such as the Poisson equation, the Helmholtz equation, and the shallow-water equations.

Section 2 describes the arrangement of equally spaced latitudinal grid points used in the new DFS method. Section 3 describes the details of the new DFS method using the new DFS expansion method and the new DFS basis functions and also includes the essential summary of the new DFS method. Section 4 examines the error due to the wavenumber truncation in the DFS method of Orszag (1974), the old DFS method (Cheong, 2000a, b; Yoshimura and Matsumura, 2005), and the new DFS method. Section 5 examines the accuracy of the old and new DFS methods and the SH method for the Laplacian operator and the Helmholtz equation. Section 6 compares the results of the shallow-water test cases between the model using the new DFS method, the model using the old DFS method of Yoshimura and Matsumura (2005), and the model using the SH method. Section 7 presents conclusions and perspectives.

2 Arrangement of equally spaced latitudinal grid points

In DFS models, equally spaced latitudinal grid points are used. We use the following three ways of arranging equally spaced latitudinal grid points in the model using the new DFS method:

(7a)Grid0:J=J0,θj=πj+0.5/J0,j=0,,J0-1,(7b)Grid1:J=J0+1,θj=πj/J0,j=0,,J0,(7c)Grid-1:J=J0-1,θj=πj/J0,j=1,,J0-1,

where θj is the colatitude at each grid point and J0 is the number of latitudinal grid points in Grid [0]. When the grid intervals in Grids [0], [1], and [1] are set as equal, the number of grid points J in Grid [1] is J0+1, and the number of grid points J in Grid [1] is J0−1. Figure 1 shows Grids [0], [1], and [1] when J0=4 and the grid interval Δθ=π/4. Grid [0] has been widely used in DFS models (e.g., Merilees, 1973b; Orszag, 1974; Cheong, 2000a, b; Yoshimura and Matsumura, 2005) and in DFS expansion (e.g., Cheong et al., 2004). Grid [1] was used in DFS expansion (e.g., Yee, 1981; Cheong et al., 2004). Grid [1] was used, for example, in the SH model using Clenshaw–Curtis quadrature (Hotta and Ujiie, 2018). All of Grids [0], [1], and [1] were used in SH expansion (Swarztrauber and Spotz, 2000).

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f01

Figure 1Grid [0], Grid[1], and Grid [1] are three ways of arranging equally spaced latitudinal grid points. Red circles show the positions of the grid points when the grid interval Δθ=π/4.

Download

In the new DFS method, the wind vector components u and v (instead of u/sinθ and v/sinθ or usin θ and vsin θ) are transformed from grid space to spectral space and vice versa, as shown in Sect. 3.5 and 3.6 below. This makes it possible to use Grid [1], which has grid points at the poles.

3 Improved double Fourier series on the sphere

In Sect. 3, we describe the new basis functions for scalar and vector variables and the new method to calculate expansion coefficients which minimizes the error due to wavenumber truncation. We compare the new DFS method with the SH method to see the difference between them. We also describe how to calculate the Laplacian operator, the Poisson equation, the Helmholtz equation, and horizontal diffusion in the new DFS method. The essential summary (cookbook) of the new DFS method is in Sect. 3.10.

3.1 New basis functions for a scalar variable

We propose the following new DFS basis functions that automatically satisfy the continuity conditions at the poles in Eq. (3). The scalar variable T(λ,θ) is expanded zonally as

(8) T λ , θ m = 0 M T m c θ cos m λ + m = 1 M T m s θ sin m λ ,

where Tmcθ and Tmsθ are calculated from T(λ,θ) by the forward Fourier transform as

(9a)Tmcθ=a2π02πcosmλTλ,θdλ,a1 for m=02 for m1,(9b)Tmsθ=1π02πsinmλTλ,θdλ.

The variables Tmcθ and Tmsθ are meridionally expanded as

(10) T m c s θ T m c s , N θ n = N min, m N max, m T n , m c s S n , m θ ,

where

(11)Sn,mθcosnθform=0,sinθcosnθform=1,sinθsinnθfor evenm2,sin2θsinnθfor oddm3,(12)Nmin,m=0form=0,0form=1,1for evenm2,1for oddm3,Nmax,m=Nform=0,N-1form=1,N-1for evenm2,N-2for oddm3.

Here, the superscript c(s) means c or s, and thus, for example, Tmcsθ means Tmcθ or Tmsθ. In Eq. (8), cos mλ and sin mλ are used instead of eimλ as zonal basis functions for convenience in calculating the expansion coefficients using the least-squares method described below in Sect. 3.3 and 3.6. In Eq. (11), the meridional basis functions sin 2θsin nθ for odd m≥3 are especially different from the basis functions of Cheong in Eq. (6). Either sin nθ or sin θcos nθ can be used as the basis functions for m=1 because it can be shown using Eq. (A2a) from Appendix A that sinθcosnθn=0,,N-1 are the linear combination of sinnθn=1,,N and vice versa. Here we use sin θcos nθ for m=1 because it can be more easily divided by sin θ, which is convenient for calculating T.

Using Eqs. (A2a)–(A2c), Eq. (10) can be converted as follows:

(13) T m c s , N θ = n = 0 N T n , m c s cos n θ for  even m , n = 1 N T n , m c s sin n θ for  odd m ,

where

(14a)Tn,mcs=Tn,mc(s)n=0,,Nfor  m=0,(14b)Tn,mcs=Tn-1,mcs-Tn+1,mcs2n=1,,N for m=1except  forT1,mcs=2T0,mcs-T2,mcs2n=1,(14c)Tn,mcs=-Tn-1,mcs+Tn+1,mcs2n=0,,Nfor  evenm2,(14d)Tn,mcs=-Tn-2,mcs+2Tn,mcs-Tn+2,mcs4n=1,,Nfor oddm3except  forT1,mcs=3T1,mcs-T3,mcs4n=1.

The value of Nmax,m in Eq. (12) is determined so that the maximum value of n for each m in Eq. (13) becomes N. In Grid [0] and Grid [1] (see Sect. 2), the upper limit of N is J0−1 for each m. In Grid [1], the upper limit of N is J0−1 for m≥2, but J0-2=J-1 for m= 0 or 1. The reason for this is shown in Appendix C.

When calculating the values of Tmcs,Nθ in grid space from Tn,mcs in spectral space, the coefficients Tn,mcs are calculated from Tn,mcs using Eq. (14), and inverse discrete cosine and sine transforms (see Appendix B) are performed using Eq. (13). The calculation of Tn,mcs in spectral space from Tmcsθ in grid space is described in Sect. 3.3 below.

The truncated variable TN,Mλ,θ is defined as

(15) T N , M λ , θ m = 0 M T m c , N θ cos m λ + m = 1 M T m s , N θ sin m λ .

From Eq. (10), the values of Tmc,Nθ at the poles are finite for m=0, and the values of Tmcs,Nθ at the poles are zero for m≠0. Therefore, TN,Mλ,θ is continuous at the poles.

3.2 Gradient of a scalar variable

The gradient TN,M=TλN,M,TϕN,M is obtained as follows:

(16a)TλN,M1asinθTN,Mλ=m=1MTλ,mc,Nθcosmλ+m=1MTλ,ms,Nθsinmλ,(16b)Tλ,mc,NθmasinθTms,Nθ=n=Nmin,mNmax,  mTn,msmSn,mθasinθ,(16c)Tλ,ms,Nθ-masinθTmc,Nθ=n=Nmin,mNmax,  m-Tn,mcmSn,mθasinθ,

(17a)TϕN,M1aTN,Mϕ=-1aTN,Mθ=m=0MTϕ,mc,Nθcosmλ+m=1MTϕ,ms,Nθsinmλ,(17b)Tϕ,mcs,Nθ-1aTmcs,Nθθ=n=Nmin,mNmax,  m-Tn,mcs1aSn,mθθ,

where a is the radius of the Earth and ϕ is the latitude. From Eqs. (16b), (16c), and (A2b), we obtain

(18) T λ , m c s , N θ = 0 for m = 0 , n = 0 N - 1 T λ , n , m c s cos n θ for m = 1 , n = 1 N - 1 T λ , n , m c s sin n θ for  even m 2 , n = 0 N - 1 T λ , n , m c s cos n θ = n = 1 N - 2 T λ , n , m c s sin θ sin n θ for odd m 3 ,

where

(19a)Tλ,n,mc=1amTn,msn=0,,N-1 for m=1,(19b)Tλ,n,mc=1amTn,msn=1,,N-1for  evenm2,(19c)Tλ,n,mc=1am-Tn-1,ms+Tn+1,ms2n=0,,N-1for  oddm3.

The equations for Tλ,n,ms are the same as Eq. (19), except that Tλ,n,mc and Tn,ms are replaced with Tλ,n,ms and -Tn,mc, respectively. From Eqs. (17b), (13), and (14), we obtain

(20) T ϕ , m c s , N θ = n = 1 N T ϕ , n , m c s sin n θ  for  m = 0 , n = 0 N T ϕ , n , m c s cos n θ  for  m = 1 , n = 1 N T ϕ , n , m c s sin n θ for  even m 2 , n = 0 N T ϕ , n , m c s cos n θ = n = 1 N - 1 T ϕ , n , m c s sin θ sin n θ for  odd m 3 ,

where

(21a)Tϕ,n,mcs=-1a-nTn,mcsn=1,,N for m=0,(21b)Tϕ,n,mcs=-1anTn-1,mcs-Tn+1,mcs2n=0,,N for m=1,except  forTϕ,1,mcs=-1a2T0,mcs-T2,mcs2n=1,(21c)Tϕ,n,mcs=-1anTn-1,mcs-Tn+1,mcs2n=1,,Nfor  even m2,(21d)Tϕ,n,mcs=-1an-Tn-2,mcs+2Tn,mcs-Tn+2,mcs4n=0,,Nfor  oddm3,except forTϕ,1,mcs=-1a3T1,mcs-T3,mcs4n=1.

From Eqs. (18)–(21), it can be seen that Tλ,mc,Nθ, Tλ,ms,Nθ, Tϕ,mc,Nθ, and Tϕ,ms,Nθ at the poles are finite for m=1 and zero for m≠1, and moreover the following relations are satisfied for m=1:

(22a)Tλ,m=1c,Nθ=-Tϕ,m=1s,Nθ=1an=0N-1Tn,m=1sat  θ=0North Pole,(22b)Tλ,m=1s,Nθ=Tϕ,m=1c,Nθ=-1an=0N-1Tn,m=1cat  θ=0North Pole,(22c)Tλ,m=1c,Nθ=Tϕ,m=1s,Nθ=1an=0N-1-1nTn,m=1satθ=πSouth Pole,(22d)Tλ,m=1s,Nθ=-Tϕ,m=1c,Nθ=-1an=0N-1-1nTn,m=1cat  θ=πSouth  Pole.

Thus, it is guaranteed that TN,M=TλN,M,TϕN,M is continuous at the poles.

3.3 New method to calculate expansion coefficients for a scalar variable

One way to calculate the coefficients Tn,mcs from Tmcsθ in Eq. (10) is to perform a forward cosine transform of Tmcsθ/sinθ for m=1, a sine transform of Tmcsθ/sinθ for even m (≥2), and a sine transform of Tmcsθ/sin2θ for odd m (≥3). However, this approach with the meridional wavenumber truncation N<J leads to large high-wavenumber oscillations (see Sect. 4). Dividing Tmcsθ by sin 2θ reduces the numerical stability of the model more significantly than dividing Tmcsθ by sin θ.

Here we propose a new method to calculate expansion coefficients using the least-squares method to minimize the error due to the meridional wavenumber truncation. This method avoids dividing Tmcsθ by sin θ or sin 2θ before the forward cosine or sine transforms. The coefficients Tn,mcs in Eq. (10) are calculated as follows. First, Tmcsθ in Eq. (10) are expanded like Eq. (5) as

(23) T m c s θ T ̃ m c s , N θ n = 0 N T ̃ n , m c s cos n θ for  even m , n = 1 N T ̃ n , m c s sin n θ for  odd m ,

where the expansion coefficients T̃n,mcs are calculated by the forward discrete cosine transform for even m and the forward discrete sine transform for odd m from the values of Tmcsθ at the grid points (see Appendix B).

Next, Tn,mc and Tn,ms are calculated using the least-squares method to minimize the following error E (the squared L2 norm of the residual):

(24) E 1 2 π 2 0 2 π 0 π R λ , θ 2 d θ d λ ,

where the residual R(λ,θ) is

(25) R λ , θ T N , M λ , θ - T λ , θ .

From Eqs. (24), (25), and (15) and the equations E/Tn,mc=0 and E/Tn,ms=0 used in the least-squares method to minimize E, we obtain

(26a)12π202π0πTmc,NθTn,mccosmλRλ,θdθdλ=0,(26b)12π202π0πTms,NθTn,mssinmλRλ,θdθdλ=0.

From Eq. (10), we derive

(27) T m c , N θ T n , m c = T m s , N θ T n , m s = S n , m θ .

Equations (26) and (27) show that the residual R(λ,θ) is orthogonal to each of the new DFS basis functions Sm,n(θ)cos mλ and Sm,n(θ)sin mλ, which means that Eq. (26) is the same as the equation derived using the Galerkin method.

From Eqs. (26), (27), (25), (15), (9), and (A3), we derive

(28) 0 π S n , m θ T m c s , N θ d θ = 0 π S n , m θ T m c s θ d θ .

From Eqs. (28) and (D4) in Appendix D, we obtain

(29) 0 π S n , m θ T m c s , N θ d θ = 0 π S n , m θ T ̃ m c s , N θ d θ .

By substituting Eqs. (10) and (23) into Eq. (29), the following equations for Tn,mc and Tn,ms are derived, as shown in Appendix E.

For m=0,

(30a) T n , m c s = T ̃ n , m c s 0 n N .

For m=1,

(30b) - T n - 2 , m c s + 2 T n , m c s - T n + 2 , m c s = - 2 T ̃ n - 1 , m c s + 2 T ̃ n + 1 , m c s 0 n N - 1 ,

with the exception of the following underlined values:

1¯T1,mcs-T3,mcs=2T̃2,mcsn=1,-2¯T0,mcs+2T2,mcs-T4,mcs=-2T̃1,mcs+2T̃3,mcsn=2.

For even m (≥2),

(30c) - T n - 2 , m c s + 2 T n , m c s - T n + 2 , m c s = 2 T ̃ n - 1 , m c s - 2 T ̃ n + 1 , m c s 1 n N - 1 ,

with the exception of the following underlined values:

3¯T1,mcs-T3,mcs=4¯T̃0,mcs-2T̃2,mcsn=1.

For odd m(≥3),

(30d) T n - 4 , m c s - 4 T n - 2 , m c s + 6 T n , m c s - 4 T n + 2 , m c s + T n + 4 , m c s = - 4 T ̃ n - 2 , m c s + 8 T ̃ n , m c s - 4 T ̃ n + 2 , m c s 1 n N - 2 ,

with the exception of the following underlined values:

10¯T1,mcs-5¯T3,mcs+T5,mcs=12¯T̃1,mcs-4T̃3,mcsn=1,5¯T2,mcs-4T4,mcs+T6,mcs=8T̃2,mcs-4T̃4,mcs      n=2,-5¯T1,mcs+6T3,mcs-4T5,mcs+T7,mcs=-4T̃1,mcs+8T̃3,mcs-4T̃5,mcsn=3.

From Eq. (30a), Tn,mcs for m=0 is obtained. From Eqs. (30d), the following linear simultaneous equations for m≥3 are derived:

(31) C n _ odd , m T 1 , m c s T 3 , m c s T 5 , m c s T 7 , m c s = D n _ odd , m T ̃ 1 , m c s T ̃ 3 , m c s T ̃ 5 , m c s T ̃ 7 , m c s , C n _ even , m T 2 , m c s T 4 , m c s T 6 , m c s T 8 , m c s = D n _ even , m T ̃ 2 , m c s T ̃ 4 , m c s T ̃ 6 , m c s T ̃ 8 , m c s ,

where the matrices Cn_odd,m and Cn_even,m are penta-diagonal. From Eqs. (30b) and (30c), the equations similar to Eq. (31) for m=1 and even m (≥2) with tri-diagonal matrices Cn_odd,m and Cn_even,m are derived. By using Eq. (31), the expansion coefficients Tn,mcs are calculated from T̃n,mcs. A penta-diagonal matrix C can be lower–upper (LU) decomposed as

(32) C = LU , L 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 , U 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 .

To solve LUx=b, we solve Ly=b with forward substitution first, and then we solve Ux=y with backward substitution. There are also other methods to solve Eq. (31). For example, the method using LU decomposition considering penta-diagonal matrices as 2×2 block tri-diagonal matrices makes SIMD (Single Instruction/Multiple Data) operations more effective. The method using cyclic reduction for block tri-diagonal matrices (e.g., Gander and Golub, 1997) is suitable for vectorization and parallelization. The calculation with these methods for each m requires O(N) operations. The simultaneous equations with tri-diagonal matrices C can be solved in a similar way. Therefore, the calculation of Tn,mcs for all m and n with Eq. (30) requires only O(N2) operations.

3.4 Comparison of new DFS with SH

Here we compare the new DFS method with the SH method to see the difference between them. In the SH method, Tmcθ andTmsθ in Eq. (8) are expanded with the associated Legendre functions Pn,m(θ) as

(33) T m c s θ T m c s ,SH , N θ n = m N T n , m c s , SH P n , m θ ,

where m≥0. The functions Pn,m(θ) satisfy the following orthogonality relations for each m:

(34) 0 π P n , m θ P n , m θ sin θ d θ = 1 ( or 2 ) for n = n , 0 for n n .

By the modified Robert expansion (Merilees, 1973a; Orszag, 1974), the associated Legendre functions Pn,m(θ) are expressed as

(35) P n , m θ = l = 0 when n - m - l is  even n - m a n , m , l sin m θ cos l θ .

Conversely, the functions sinmθcosn-mθnm can be expressed as the linear combination of Pl,mθl=m,,n. Substituting Eq. (35) into Eq. (33) gives

(36) T m c s , SH , N θ = n = 0 N - m T n , m c s , SH sin m θ cos n θ ,

where m≥0. Equation (36) is similar to Eq. (10) in the following sense: the basis functions for m=0 and m=1 in Eq. (36) are the same as Eq. (11). The basis functions sin2θcosnθn=0,,N-2 for m=2 and sin3θcosnθn=0,,N-3 for m=3 in Eq. (36) are the linear combinations of sinθsinnθn=1,,N-1 and sin2θsinnθn=1,,N-2 in Eq. (11), respectively (see Eq. A2a), and vice versa. The basis functions for m≥4 in Eq. (36) are different from Eq. (11). The number of expansion coefficients in Eq. (33) or Eq. (36) in the SH method is smaller than in Eq. (10) in the new DFS method for each m≥4. From Eqs. (8) and (33), the number of expansion coefficients Tn,mc,SH in the SH model is about N2/2 when M=N. This triangular truncation used in the SH method gives a uniform resolution over the sphere. From Eqs. (8) and (10), the number of the expansion coefficients Tn,mc in the DFS method is about N2 when M=N. This rectangular truncation used in the DFS model gives almost the same resolution as the grid spacing of the regular longitude–latitude grids. Therefore, the zonal Fourier filter (see Appendix F) is used in the DFS model to give a more uniform resolution.

We compare the method used to calculate the expansion coefficients in the new DFS method with that in the SH method. The SH expansion coefficients Tn,mcs,SH in Eq. (33) are calculated from Tmcsθ by the forward Legendre transform as

(37) T n , m c s , SH = 0 π P n , m θ T m c s θ sin θ d θ ,

where Gaussian quadrature or Clenshaw–Curtis quadrature (e.g., Hotta and Ujiie, 2018) is usually used for integration. They can also be calculated using the same equations as Eq. (37), except that Tmcsθ is replaced with T̃mcs,Nθ (e.g., Sneeuw and Bun, 1996), although the values of Tn,mcs,SH calculated from T̃mcs,Nθ are different from those calculated from Tmcsθ. In the new DFS method, the values of Tn,mcs calculated from T̃mcs,Nθ in Eq. (29) are the same as those calculated from T̃mcsθ in Eq. (28) (see Eq. D4 in Appendix D).

Equation (37) can be derived using the least-squares method that minimizes the error ESH (the squared L2 norm of the residual):

(38) E SH 1 4 π 0 2 π 0 π R SH λ , θ 2 sin θ d θ d λ ,

where the residual RSH(λ,θ) is

(39) R SH λ , θ m = 0 M T m c , SH , N θ cos m λ + m = 1 M T m s , SH , N θ sin m λ - T λ , θ .

From Eqs. (38), (39), and (33) and the equations ESH/Tn,mc,SH=0 and ESH/Tn,ms,SH=0 used in the least-squares method to minimize ESH, we derive

(40a)02π0πPn,mθcosmλRSHλ,θsinθdθdλ=0,(40b)02π0πPn,mθsinmλRSHλ,θsinθdθdλ=0.

Equation (40) is the same as the equation obtained using the Galerkin method. From Eqs. (40), (33), (34), (9), and (A3), we derive Eq. (37).

In Eqs. (37) and (38), the latitudinal weight sin θ appears, unlike in Eqs. (24) and (28), which is another difference between the SH and the new DFS methods. In the DFS method, the constant latitudinal weight is used in Eq. (24), although the latitudinal area weight described below in Appendix G is usually used as the latitudinal weight at the grid points, for example, for the calculation of the global mean.

When calculating the coefficients Tn,mcs in Eq. (10), we can also consider the least-squares method, not by using E in Eq. (24) but instead by using E with latitudinal weight sin θ like in Eq. (38). However, minimizing E derives the simultaneous equations for calculating Tn,mcs with dense matrices, which leads to O(N3) operations. When using E, the simultaneous equations with penta-diagonal or tri-diagonal matrices require only O(N2) operations. Therefore, we choose to use E instead of E.

The new DFS meridional basis functions Sn,m(θ) for each m are not orthogonal but independent. Therefore, by using Gram–Schmidt orthogonalization, the basis functions can be converted to orthogonalized basis functions Sn,mOθ, which satisfy

(41) 1 π 0 π S n , m O θ S n , m O θ d θ = 1 for n = n , 0 for n n .

This is similar to Eq. (34), but the latitudinal weight is constant. Tmcsθ in Eq. (8) are expanded with Sn,mOθ as

(42) T m c s θ T m c s , O , N θ n = N min, m N max,   m T n , m c s ,O S n , m O θ .

By using the least squares method or the Galerkin method with Eq. (42), we obtain the same equations as Eqs. (24)–(29), except that Tmcs,Nθ and Sn,m(θ) are replaced with Tmcs,O,Nθ and Sn,mOθ, respectively. From Eq. (29) with Tmcs,Nθ and Sn,m(θ) replaced by Tmcs,O,Nθ and Sn,mOθ and Eqs. (41) and (42), we derive

(43) T n , m c s ,O = 1 π 0 π S n , m O θ T ̃ m c s , N θ d θ .

Thus, Tn,mcs,O and Tmcs,O,Nθ in Eqs. (43) and (42) are calculated uniquely. This unique solution Tmcs,O,Nθ is the same as Tmcs,Nθ in Eq. (29) obtained by the least-squares method with the non-orthogonal basis functions Sn,m(θ) because Sn,mOθn=Nmin,m,,Nmax,  m is the linear combination of Sn,mθn=Nmin,m,,Nmax,  m for each m and vice versa.

3.5 New basis functions for a wind vector

The velocity potential χ and the stream function ψ can be converted into the wind vector components u and v using the equations

(44a)u=1acosϕχλ-1aψϕ=1asinθχλ+1aψθ,(44b)v=1acosϕψλ+1aχϕ=1asinθψλ-1aχθ,

where u=acosϕdλ/dt is the zonal wind and v=adϕ/dt is the meridional wind. The scalar variables χ and ψ are expanded like Eqs. (8) and (10) as

(45)χλ,θψλ,θm=0Mχmcθψmcθcosmλ+m=1Mχmsθψmsθsinmλ,(46)χmcsθψmcsθχmcs,Nθψmcs,Nθn=Nmin,mNmax,  mχn,mcsψn,mcsSn,mθ,

The truncated variables ψN,Mλ,θ and χN,Mλ,θ are defined as

(47) χ N , M λ , θ ψ N , M λ , θ m = 0 M χ m c , N θ ψ m c , N θ cos m λ + m = 1 M χ m s , N θ ψ m s , N θ sin m λ ,

From Eqs. (44)–(47), the equations for the wind vector components uN,Mλ,θ and vN,Mλ,θ are derived as

(48a)uN,Mλ,θ1asinθχN,Mλ,θλ+1aψN,Mλ,θθ=m=0Mumc,Nθcosmλ+m=1Mums,Nθsinmλ,(48b)umc,Nθmχms,Nθasinθ+1aψmc,Nθθ=n=Nmin,mNmax,  mχn,msmSn,mθasinθ+ψn,mc1aSn,mθθ,(48c)ums,Nθ-mχmc,Nθasinθ+1aψms,Nθθ=n=Nmin,mNmax,  m-χn,mcmSn,mθasinθ+ψn,ms1aSn,mθθ,

(49a)vN,Mλ,θ1asinθψN,Mλ,θλ-1aχN,Mλ,θθ=m=0Mvmc,Nθcosmλ+m=1Mvms,Nθsinmλ,(49b)vmc,Nθmψms,Nθasinθ-1aχmc,Nθθ=n=Nmin,mNmax,  mψn,msmSn,mθasinθ-χn,mc1aSn,mθθ,(49c)vms,Nθ-mψmc,Nθasinθ-1aχms,Nθθ=n=Nmin,mNmax,  m-ψn,mcmSn,mθasinθ-χn,ms1aSn,mθθ.

The vector uN,M,vN,M in Eqs. (48) and (49) can also be represented as

(50) u N , M λ , θ , v N , M λ , θ = m = 0 M n = N min, m N max,   m χ n , m c V n , m 1 + χ n , m s V n , m 2 + ψ n , m c V n , m 3 + ψ n , m s V n , m 4 ,

where we define the new DFS vector basis functions Vn,m1, Vn,m2, Vn,m3, and Vn,m4 as

(51a)Vn,m1λ,θ-mSn,mθasinθsinmλ,-1aSn,mθθcosmλ,(51b)Vn,m2λ,θmSn,mθasinθcosmλ,-1aSn,mθθsinmλ,(51c)Vn,m3λ,θ1aSn,mθθcosmλ,-mSn,mθasinθsinmλ,(51d)Vn,m4λ,θ1aSn,mθθsinmλ,mSn,mθasinθcosmλ.

From Eqs. (48), (49), and (16)–(21), we obtain

(52) u m c s , N θ v m c s , N θ = n = 1 N u n , m c s v n , m c s sin n θ for   m = 0 , n = 0 N u n , m c s v n , m c s cos n θ for  m = 1 , n = 1 N u n , m c s v n , m c s sin n θ for  even  m 2 , n = 0 N u n , m c s v n , m c s cos n θ = n = 1 N - 1 u n , m c s v n , m c s sin θ sin n θ  for odd   m 3 ,

where

(53a)un,mc=1a-nψn,mcn=1,,Nfor  m=0,(53b)un,mc=1amχn,ms+nψn-1,mc-ψn+1,mc2n=0,,Nfor  m=1,except  for  u1,mc=1amχ1,ms+2ψ0,mc-ψ2,mc2n=1,(53c)un,mc=1amχn,ms+nψn-1,mc-ψn+1,mc2n=1,,Nfor  even  m2,(53d)un,mc=1am-χn-1,ms+χn+1,ms2+n-ψn-2,mc+2ψn,mc-ψn+2,mc4n=0,,Nfor  odd  m3except  for  u1,mc=1amχ2,ms2+3ψ1,mc-ψ3,mc4n=1.

The equations for un,ms are the same as Eqs. (53b)–(53d), except that un,mc, χn,ms, and ψn,mc are replaced with un,ms, -χn,mc, and ψn,ms, respectively. The equations for vn,mc are the same as Eqs. (53a)–(53d), except that un,mc, χn,ms, and ψn,mc are replaced with vn,mc, ψn,ms, and -χn,mc, respectively. The equations for vn,ms are the same as Eqs. (53b)–(53d), except that un,mc, χn,ms, and ψn,mc are replaced with vn,ms, -ψn,mc, and -χn,ms, respectively.

From Eqs. (52) and (53), it can be seen that umc,Nθ, ums,Nθ, vmc,Nθ, and vms,Nθ at the poles are finite for m=1 and zero for m≠1. Moreover, the following relations are satisfied for m=1:

(54a)um=1c,Nθ=-vm=1s,Nθ=1an=0N-1χn,m=1s+ψn,m=1cat  θ=0North  Pole,(54b)um=1s,Nθ=vm=1c,Nθ=1an=0N-1-χn,m=1c+ψn,m=1sat  θ=0North  Pole,(54c)um=1c,Nθ=vm=1s,Nθ=1an=0N-1-1nχn,m=1s-ψn,m=1cat  θ=πSouth  Pole,(54d)um=1s,Nθ=-vm=1c,Nθ=1an=0N-1-1n-χn,m=1c-ψn,m=1sat  θ=πSouth  Pole.

Thus, it is guaranteed that the wind vector uN,M,vN,M in Eqs. (48) and (49) is continuous at the poles.

3.6 New method to calculate expansion coefficients for a wind vector

We propose a new method that calculates the expansion coefficients χn,mc, χn,ms, ψn,mc, and ψn,ms in Eqs. (48)–(50) using the least-squares method to minimize the error of uN,Mλ,θ and vN,Mλ,θ from u(λ,θ) and v(λ,θ) due to the meridional wavenumber truncation. First, the wind vector components u and v are expanded zonally as

(55) u λ , θ v λ , θ m = 0 M u m c θ v m c θ cos m λ + m = 1 M u m s θ v m s θ sin m λ ,

where umcsθ and vmcsθ are calculated from u(λ,θ) and v(λ,θ) by the forward Fourier transform as

(56a)umcθvmcθ=a2π02πcosmλuλ,θvλ,θdλ,a1 form=02 form1,(56b)umsθvmsθ=1π02πsinmλuλ,θvλ,θdλ.

The variables umcsθ and vmcsθ are meridionally expanded as

(57) u m c s θ v m c s θ u ̃ m c s , N θ v ̃ m c s , N θ n = 1 N u ̃ n , m c s v ̃ n , m c s sin n θ for  even   m , n = 0 N u ̃ n , m c s v ̃ n , m c s cos n θ for  odd m ,

where ũn,mcs and ṽn,mcs are calculated from umcsθ and vmcsθ by the forward discrete cosine or sine transform (see Appendix B).

Next, χn,mc, χn,ms, ψn,mc, and ψn,ms are calculated to minimize the following error F (the squared L2 norm of the residual vector):

(58) F 1 2 π 2 0 2 π 0 π R u λ , θ 2 + R v λ , θ 2 d θ d λ ,

where the residual vector Ruλ,θ,Rvλ,θ is defined as

(59a)Ruλ,θuN,Mλ,θ-uλ,θ,(59b)Rvλ,θvN,Mλ,θ-vλ,θ.

From Eqs. (58) and (59) and the equations F/χm,nc=0, F/χn,ms=0, F/ψn,mc=0, and F/ψn,ms=0 used in the least-squares method, we obtain

(60a)12π202π0πuN,Mλ,θχn,mcRuλ,θ+vN,Mλ,θχn,mcRvλ,θdθdλ=0,(60b)12π202π0πuN,Mλ,θχn,msRuλ,θ+vN,Mλ,θχn,msRvλ,θdθdλ=0,(60c)12π202π0πuN,Mλ,θψn,mcRuλ,θ+vN,Mλ,θψn,mcRvλ,θdθdλ=0,(60d)12π202π0πuN,Mλ,θψn,msRuλ,θ+vN,Mλ,θψn,msRvλ,θdθdλ=0.

From Eq. (50), we derive

(61a)uN,Mλ,θχn,mc,vN,Mλ,θχn,mc=Vn,m1λ,θ,(61b)uN,Mλ,θχn,ms,vN,Mλ,θχn,ms=Vn,m2λ,θ,(61c)uN,Mλ,θψn,mc,vN,Mλ,θψn,mc=Vn,m3λ,θ,(61d)uN,Mλ,θψn,ms,vN,Mλ,θψn,ms=Vn,m4λ,θ.

Equations (60) and (61) show that the residual vector Ruλ,θ,Rvλ,θ is orthogonal to each of the vector basis function, which means that Eq. (60) is the same as the equation obtained by the Galerkin method. From Eqs. (60), (61), (51), (48a), (49a), (56), (A3), and (D6), we derive

(62a)1π0π-mSn,mθasinθums,Nθ-ũms,Nθ-1aSn,mθθvmc,Nθ-ṽmc,Nθdθ=0,(62b)1π0πmSn,mθasinθumc,Nθ-ũmc,Nθ-1aSn,mθθvms,Nθ-ṽms,Nθdθ=0,(62c)1π0π1aSn,mθθumc,Nθ-ũmc,Nθ-mSn,mθasinθvms,Nθ-ṽms,Nθdθ=0,(62d)1π0π1aSn,mθθums,Nθ-ũms,Nθ+mSn,mθasinθvmc,Nθ-ṽmc,Nθdθ=0.

By substituting Eqs. (52) and (57) into Eq. (62a), the following equations for χn,mc and ψn,ms are derived as shown in Appendix H.

For m=0,

(63a) 1 a n χ n , m c = v ̃ n , m c 1 n N .

The coefficient χm=0,n=0c is determined so that the global means of χ are zero. See Eq. (G1) for an explanation of the calculation of the global mean.

For m=1,

(63b) 1 a - n - 1 2 χ n - 2 , m c - 2 m ψ n - 1 , m s + 4 m 2 + 2 n 2 + 2 χ n , m c - 2 m ψ n + 1 , m s - n + 1 2 χ n + 2 , m c = 2 n - 1 v ̃ n - 1 , m c - 4 m u ̃ n , m s - 2 n + 1 v ̃ n + 1 , m c 0 n N - 1 ,

with the exception of the following underlined values:

1a8¯m2+4¯χ0,mc-4¯mψ1,ms-2¯χ2,mc=-8¯mũ0,ms-4¯ṽ1,mcn=0,1a-4¯mψ0,ms+4m2+4χ1,mc+=n=1,1a-2¯χ0,mc-2mψ1,ms+=n=2.

For even m≥2,

(63c) 1 a - n - 1 2 χ n - 2 , m c - 2 m ψ n - 1 , m s + 4 m 2 + 2 n 2 + 2 χ n , m c - 2 m ψ n + 1 , m s - n + 1 2 χ n + 2 , m c = 2 n - 1 v ̃ n - 1 , m c - 4 m u ̃ n , m s - 2 n + 1 v ̃ n + 1 , m c 1 n N - 1 ,

with no exception.

For odd m≥3,

(63d) 1 a n - 2 2 χ n - 4 , m c + 2 m ψ n - 3 , m s + - 4 m 2 - 4 n 2 + 8 n - 8 χ n - 2 , m c - 2 m ψ n - 1 , m s + 8 m 2 + 6 n 2 + 8 χ n , m c - 2 m ψ n + 1 , m s + - 4 m 2 - 4 n 2 - 8 n - 8 χ n + 2 , m c + 2 m ψ n + 3 , m s + n + 2 2 χ n + 4 , m c = 4 n - 2 v ̃ n - 2 , m c - 8 m u ̃ n - 1 , m s - 8 n v ̃ n , m c + 8 m u ̃ n + 1 , m s + 4 n + 2 v ̃ n + 2 , m c 1 n N - 2 ,

with the exception of the following underlined values:

1a12¯m2+18¯χ1,mc-4¯mψ2,ms+-4m2-21¯χ3,mc+=-16¯mũ0,ms-12¯ṽ1,mc+n=1,1a-4¯mψ1,ms+8m2+32χ2,mc+=n=2,1a-4m2-21¯χ1,mc-2mψ2,ms+=n=3.

Similarly, from Eq. (62b) we derive the same equations as Eqs. (63b)–(63d), except that χc, ψs, ṽc, and ũs are replaced with χs, ψc, ṽs, and -ũc, respectively. From Eq. (62c), we derive the same equations as Eqs. (63a)–(63d), except that χc, ψs, ṽc, and ũs are replaced with ψc, χs, ũc, and -ṽs, respectively. From Eq. (62d), we derive the same equations as Eqs. (63b)–(63d), except that χc, ψs, ṽc, and ũs are replaced with ψs, χc, -ũs, and -ṽc, respectively.

Equation (63a) is easily solved. From Eq. (63d) and from the same equations as Eq. (63d), except that χc, ψs, ṽc, and ũs are replaced with ψs, χc, -ũs, and -ṽc, respectively, we derive the following linear simultaneous equations for m≥3:

(64) E m χ 1 , m c ψ 2 , m s χ 3 , m c ψ 4 , m s : = F m u ̃ 0 , m s v ̃ 1 , m c u ̃ 2 , m s v ̃ 3 , m c : , E m ψ 1 , m s χ 2 , m c ψ 3 , m s χ 4 , m c : = F m - v ̃ 0 , m c - u ̃ 1 , m s - v ̃ 2 , m c - u ̃ 3 , m s : ,

where the matrices Em are nona-diagonal. From Eqs. (63b) and (63c), we derive equations similar to Eq. (64) for m=1 and even m (≥2) with penta-diagonal matrices Em. The simultaneous equations with nona-diagonal or penta-diagonal matrices Em can be solved in a similar way to Eq. (31), and the expansion coefficients χn,mc and ψn,ms in Eq. (64) can be calculated efficiently. From the same equations as Eqs. (63b)–(63d), except that χc, ψs, ṽc, and ũs are replaced with χs, ψc, ṽs, and -ũc, respectively, and the same equations as Eqs. (63b)–(63d), except that χc, ψs, ṽc, and ũs are replaced with ψc, χs, ũc, and -ṽs, respectively, the simultaneous equations similar to Eq. (64) are also derived. Thus, the expansion coefficients χn,mc, χn,ms, ψn,mc, and ψn,ms are calculated from ũn,mc, ũn,ms, ṽn,mc, and ṽn,ms using Eqs. (63a)–(63d) and the similar equations. The expansion coefficients un,mc, un,ms, vn,mc, and vn,ms are calculated from χn,mc, χn,ms, ψn,mc, and ψn,ms using Eq. (53) for un,mc and the similar equations for un,ms, vn,mc, and vn,ms.

This method to calculate the DFS expansion coefficients of χ and ψ from u and v using the least-squares method (or the Galerkin method with the DFS vector basis functions) is similar to the vector harmonic transform method (Browning et al., 1989; Swarztrauber, 1993), where the SH expansion coefficients of the divergence D=∇2χ and the vorticity ζ=∇2ψ are calculated from the grid-point values of u and v using the Galerkin spectral method with the orthogonal vector SH basis functions.

3.7 Laplacian operator and Poisson equation

The calculation of the Laplacian operator and the Poisson equation in the new DFS method is described here. In the equation

(65) g λ , θ = 2 f λ , θ = 1 a 2 1 sin 2 θ 2 f λ 2 + 1 sin θ θ sin θ f θ ,

where 2 is the Laplacian operator, the variables f and g are expanded zonally like Eq. (8) as

(66) f λ , θ g λ , θ m = 0 M f m c θ g m c θ cos m λ + m = 1 M f m s θ g m s θ sin m λ .

The variables fmcθ, fmsθ, gmcθ, and gmsθ are expanded meridionally like Eq. (10) as

(67) f m c s θ g m c s θ f m c s , N θ g m c s , N θ n = N min, m N max,   m f n , m c s g n , m c s S n , m θ .

We define the truncated variables fN,M(θ) and gN,M(θ) as

(68) f N , M λ , θ g N , M λ , θ m = 0 M f m c , N θ g m c , N θ cos m λ + m = 1 M f m s , N θ g m s , N θ sin m λ .

From Eqs. (65) and (68), we obtain

(69) 2 f N , M λ , θ = m = 0 M 1 a 2 - m 2 sin 2 θ f m c, N θ + 1 sin θ θ sin θ f m c, N θ θ cos m λ + m = 1 M 1 a 2 - m 2 sin 2 θ f m s , N θ + 1 sin θ θ sin θ f m s, N θ θ sin m λ .

Here we use the Galerkin method to calculate the Laplacian operator and the Poisson equation and obtain

(70) 1 2 π 2 0 2 π 0 π S n , m θ cos m λ sin m λ R g λ , θ d θ d λ = 0 ,

where the residual is as follows:

(71) R g λ , θ g N , M λ , θ - 2 f N , M λ , θ ,

and it is orthogonal to each of the new DFS basis functions Sm,n(θ)cos mλ and Sm,n(θ)sin mλ.

We can also use the least-squares method instead of the Galerkin method so that the following error H (the squared L2 norm of the residual) is minimized:

(72) H 1 2 π 2 0 2 π 0 π R g λ , θ 2 d θ d λ .

When calculating g by applying the Laplacian operator to a given f, gn,mc and gn,ms can also be calculated from H/gn,mc and H/gn,ms using the least-squares method. The equations H/gn,mc and H/gn,ms give the equivalent equations to Eq. (70). When calculating f from a given g in the Poisson equation, fn,mc and fn,ms can also be calculated from H/fn,mc and H/fn,ms using the least-squares method. However, the equations derived from H/fn,mc and H/fn,ms are different from Eq. (70). If we use different equations for calculating g from f and f from g, the original values are changed when calculating g from f followed by calculating f from g, which may be not good for numerical stability. Therefore, we use Eq. (70) obtained with the Galerkin method for calculating both g from f and f from g. Generally, it cannot be said that the least-squares method is superior to the Galerkin method or vice versa, and here we choose to use the Galerkin method because of the reason described above.

From Eqs. (68)–(71) and (A3) we derive

(73) 0 π S n , m θ g m c s , N θ - 1 a 2 - m 2 sin 2 θ f m c s , N θ + 1 sin θ θ sin θ f m c s , N θ θ d θ = 0 .

For m=0, we calculate gn,mcs by using

(74) g m c s , N θ = 1 a 2 - m 2 sin 2 θ f m c s , N θ + 1 sin θ θ sin θ f m c s , N θ θ ,

instead of Eq. (73) following Yee (1981) and Cheong (2000a) for ease of calculation. For 0m3, the exact solutions of gn,mcs can be obtained from Eq. (74) because the new DFS meridional basis functions for 0m3 are the linear combination of the associated Legendre functions for 0m3 and vice versa as described in Sect. 3.4.

For m=0, by substituting Eq. (67) into Eq. (74) multiplied by sin 2θ, transforming using Eqs. (A2d) and (A5b), and comparing both sides of the equation, we obtain

(75a) - g n - 2 , m c s + 2 g n , m c s - g n + 2 , m c s = 1 a 2 n - 1 n - 2 f n - 2 , m c s - 2 n 2 f n , m c s + n + 1 n + 2 f n + 2 , m c s 0 n N ,

except for the following underlined values:

1¯g1,mcs-g3,mcs=n=1,-2¯g0,mcs+2g2,mcs-g4,mcs=n=2.

For m=1, by substituting Eq. (67) into Eq. (73) and using Eqs. (A2d), (A4a), and (A5b), we obtain

(75b) - g n - 2 , m c s + 2 g n , m c s - g n + 2 , m c s = 1 a 2 n - 1 n f n - 2 , m c s - 2 n 2 + 4 m 2 f n , m c s + n + 1 n f n + 2 , m c s 0 n N - 1 ,

except for the following underlined values:

1¯g1,mcs-g3,mcs=n=1,-2¯g0,mcs+2g2,mcs-g4,mcs=1a24¯f0,mcs+n=2.

For even m≥2, by substituting Eq. (67) into Eq. (73) and using Eqs. (A2c), (A4b), and (A5d), we obtain

(75c) - g n - 2 , m c s + 2 g n , m c s - g n + 2 , m c s = 1 a 2 n - 1 n f n - 2 , m c s - 2 n 2 + 4 m 2 f n , m c s + n + 1 n f n + 2 , m c s 1 n N - 1 ,

except for the following underlined values:

3¯g1,mcs-g3,mcs=n=1.

For odd m≥3, by substituting Eq. (67) into Eq. (73) and using Eqs. (A2c), (A2e), (A4b), and (A5d), we obtain

(75d) g n - 4 , m c s - 4 g n - 2 , m c s + 6 g n , m c s - 4 g n + 2 , m c s + g n + 4 , m c s = 1 a 2 - n - 2 n - 1 f n - 4 , m c s + 4 n 2 - 6 n + 4 + 4 m 2 f n - 2 , m c s - 6 n 2 + 4 + 8 m 2 f n , m c s + 4 n 2 + 6 n + 4 + 4 m 2 f n + 2 , m c s - n + 2 n + 1 f n + 4 , m c s 1 n N - 2 ,

except for the following underlined values:

10¯g1,mcs-5¯g3,mcs+g5,mcs=1a2-12¯+12¯m2f1,mcs+n=1,5¯g2,mcs-4g4,mcs+g6,mcs=n=2,-5¯g1,mcs+6g3,mcs-4g5,mcs+g7,mcs=1a224¯+4m2f1,mcs+n=3.

From Eq. (75), we obtain the following two linear simultaneous equations with tri-diagonal or penta-diagonal matrices:

(76) A n _ even , m g n _ even , m c s = B n _ even , m f n _ even , m c s , A n _ odd , m g n _ odd , m c s = B n _ odd , m f n _ odd , m c s

where gn_even,mcs and gn_odd,mcs are the vectors whose components are gn,mcs (n is even) and gn,mcs (n is odd), respectively, and fn_even,mcs and fn_odd,mcs are the vectors whose components are fn,mcs (n is even) and fn,mcs (n is odd), respectively; An_even,m, Bn_even,m, An_odd,m and Bn_odd,m are tri-diagonal or penta-diagonal matrices. gn_even,mcs and gn_odd,mcs are calculated by

(77) g n _ even , m c s = A n _ even , m - 1 B n _ even , m f n _ even , m c s , g n _ odd , m c s = A n _ odd , m - 1 B n _ odd , m f n _ odd , m c s ,

which can be solved efficiently as in Eq. (31). We have verified that all the eigenvalues of the matrices An_even,m-1Bn_even,m and An_odd,m-1Bn_odd,m are negative real numbers for several truncation wavenumbers M and N, but we have not yet proven that this is true for all truncation wavenumbers.

In the Poisson equation, f is calculated from given g in Eq. (65). We calculate f from g by the reverse calculation of g from f in Eq. (77). That is, we calculate f from g by

(78) f n _ even , m c s = B n _ even , m - 1 A n _ even , m g n _ even , m c s , f n _ odd , m c s = B n _ odd , m - 1 A n _ odd , m g n _ odd , m c s ,

except when m=0 and n is even. For m=0, fn=0,m=0c disappears in Eq. (75a). The coefficients fn,m=0c (even n≥2) are calculated from gn,m=0c (even n≥2) by using Eq. (75a). The value fn=0,m=0c is calculated from fn,m=0c (even n≥2) so that the global mean of f is zero using Eq. (G1).

In Eq. (65), the global mean of g must be zero because the global mean of the right-hand side of Eq. (65) is zero. Before calculating f from a given g in the Poisson equation, we should subtract the global mean from g (Cheong, 2000b). See Eq. (G1) for an explanation of the calculation of the global mean.

3.8 The Helmholtz equation

The Helmholtz equation is

(79) f - ε 2 f = 1 - ε 1 a 2 1 sin 2 θ 2 λ 2 + 1 sin θ θ sin θ θ f = g ,

where f is calculated from given g. From Eq. (76), the Poisson equation in Eq. (65) is represented as

(80) A g = B f ,

where the subscripts n_even, n_odd, and m and the superscripts c and s are omitted. Similarly, by using the Galerkin method, Eq. (79) is represented as

(81) A f - ε B f = A g .

From Eq. (81), f is calculated from g by

(82) f = A - ε B - 1 A g .

Since AεB is a penta-diagonal or tri-diagonal matrix, Eq. (82) can be efficiently solved as in Eq. (31).

Similarly, the Helmholtz-like equation

(83) f - ε 2 f = 2 g ,

is represented as

(84) A f - ε B f = B g .

From Eq. (84), f is calculated from g by

(85) f = A - ε B - 1 B g .

3.9 Horizontal diffusion

The horizontal diffusion is calculated in the similar way as in Cheong et al. (2004). Here we describe how to calculate fourth-order diffusion. Higher-order diffusion can be calculated similarly.

The equation for fourth-order hyper-diffusion is

(86) f + ε 4 f = g ,

where f is calculated from g. Equation (86) can be converted into

(87) 1 + i ε 2 1 - i ε 2 f = g ,

where i=-1. The calculation of Eq. (86) is accomplished by successive calculations of the following Helmholtz equations:

(88a)1+iε2f=g,(88b)1-iε2f=f,

which are represented as

(89a)A+iεBf=Ag,(89b)A-iεBf=Af.

From Eqs. (89), we obtain the equation to calculate f from g as

(90) f = A - i ε B - 1 A A + i ε B - 1 A g .

Here, A-iεB and A+iεB are complex matrices and f and g are real column vectors. For efficient computation, two real column vectors can be combined into one complex column vector (Cheong et al., 2004); for example, f=fc-ifs and g=gc-igs, where the superscript c indicates the zonal cosine component, and the superscript s indicates the zonal sine component.

3.10 Essential summary (cookbook) of the new DFS method

The essential summary for a scalar variable is as follows.

  1. Define DFS expansion for a scalar variable with zonal expansion in Eq. (8) and meridional expansion in Eq. (10).

  2. For the inverse transform from spectral space to grid point space,

    • a.

      calculate the coefficients Tn,mcs from Tn,mc(s) by using Eqs. (14),

    • b.

      calculate Tmcs,Nθ in Eq. (13) from Tn,mcs by inverse cosine and sine Fourier transforms in Appendix B,

    • c.

      calculate the grid point values TN,Mλ,θ in Eq. (15) from Tmcs,Nθ by inverse Fourier transform.

  3. For the forward transform from grid point space to spectral space,

    • a.

      calculate Tmcsθ in Eq. (9) from the grid point values T(λ,θ) by forward Fourier transform,

    • b.

      calculate the coefficients T̃n,mcs in Eq. (23) from Tmcsθ by forward cosine and sine transforms in Appendix B,

    • c.

      calculate the coefficients Tn,mcs from T̃n,mcs by using Eqs. (30) and (31), where the coefficients Tn,mcs are calculated so that Eq. (29) derived from the Galerkin method (or the least-squares method) is satisfied.

The essential summary for a vector variable is as follows.

  1. Represent DFS expansion for a vector variable by Eq. (50).

  2. For the inverse transform from spectral space to grid point space,

    • a.

      calculate the coefficients un,mcs and vn,mcs from χn,mcs and ψn,mcs by using Eq. (53) and the similar equations,

    • b.

      calculate umcs,Nθ and vmcs,Nθ in Eq. (52) from un,mcs and vn,mcs by inverse cosine and sine transforms in Appendix B,

    • c.

      calculate the grid point values uN,Mλ,θ and vN,Mλ,θ in Eqs. (48a) and (49a) from umcs,Nθ and vmcs,Nθ by inverse Fourier transform.

  3. For the forward transform from grid point space to spectral space,

    • a.

      calculate umcsθ and vmcsθ in Eq. (56) from the grid point values u(λ,θ) and v(λ,θ) by forward Fourier transform,

    • b.

      calculate the coefficients ũn,mcs and ṽn,mcs in Eq. (57) from the grid point values umcsθ and vmcsθ by forward cosine and sine transforms in Appendix B,

    • c.

      calculate the coefficients χn,mcs and ψn,mcs from ũn,mcs and ṽn,mcs by using Eqs. (63) and (64), where χn,mcs and ψn,mcs are calculated so that Eq. (62) derived from the Galerkin method (or the least-squares method) is satisfied.

4 The error due to meridional wavenumber truncation in DFS expansion methods

Here we examine the error due to the meridional wavenumber truncation when the same continuity conditions at the poles as Eq. (3) are satisfied. In the DFS method of Orszag (1974) using Eq. (2), only fN-1,m and fN,m are modified to satisfy Eq. (4) equivalent to Eq. (3). In the old DFS method using Eq. (6), which is proposed in Cheong (2000a, b) and used in Yoshimura and Matsumura (2005), the DFS meridional basis functions automatically satisfy the pole conditions in Eq. (3) for even m but not for odd m. In the new DFS method using Eqs. (10)–(12), the DFS meridional basis functions automatically satisfy the condition in Eq. (3) for both even and odd m. We examine the error due to the wavenumber truncation in these DFS methods while comparing it with the SH method.

Figure 2 shows the error due to the wavenumber truncation. The number of latitudinal grid points is J=64. The initial values of Fm(θj) are set to one at the grid points north of 30 N (except for the North Pole), and zero at the grid points south of 30 N. Grid [0] is used in the DFS methods, and the Gaussian grid is used in the SH method. There are no grid points at the poles. Since the values at the poles are zero due to the pole conditions in Eq. (3), the initial values abruptly change around the North Pole. The initial values are meridionally transformed from grid space to spectral space (forward transform), truncated with N=42, and then transformed back from spectral space to grid space (inverse transform) to obtain the truncated reconstruction of Fm(θj).

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f02

Figure 2Change in values at the grid points due to the meridional wavenumber truncation. We use Grid [0] with the number of latitudinal grid points J=64. Initial values (black) are meridionally transformed from grid space to spectral space, truncated with N=42, and transformed back from spectral space to grid space. (a) Values for even |m|≥2 when using the DFS method of Orszag (blue), the old DFS method (green), and the new DFS method (red) with Grid [0]. (b) Values for m=2 (orange), m=14 (deep sky blue), and m=30 (lime) when using the SH expansion method with the Gaussian grid. Panel (c) is the same as (a) except for the values for odd |m|≥3. Panel (d) is the same as (b) except that m=3 (orange), m=15 (deep sky blue), and m=31 (lime).

Download

In the DFS method of Orszag, a very large error occurs, especially for odd |m| (≥3) (Fig. 2c), when fN-1,m and fN,m are modified to satisfy the pole conditions in Eq. (4). Dividing Fm(θj) by sin θ before the forward Fourier cosine transform for odd m also contributes to the large error.

In the old DFS method, large high-wavenumber oscillations appear for even m (≠0) in Fig. 2a. Although the basis functions for even m (≠0) in the old DFS method are the same as those in the new method, the expansion coefficients are calculated differently in the two methods. In the old DFS method, the simple meridional truncation with N<J after the forward Fourier sine transform of a variable divided by sin θ causes the large high-wavenumber oscillations. The large oscillations appear especially when the initial values abruptly change around the poles. In the case shown in Fig. 2, the initial values at the grid points near the North Pole are one, but the value at the North Pole abruptly becomes zero due to the pole conditions of Eq. (3). The result in the old DFS method for odd |m| (≥3) is not shown in Fig. 2c because the method does not satisfy the condition of Eq. (3) for odd m

In the new DFS method, the usual small oscillations from the Gibbs phenomenon appear in Fig. 2. The error is small because the expansion coefficients are calculated using the least-squares method (or the Galerkin method) to minimize the error. Because of this, the truncation with arbitrary N<J does not cause large oscillations in the new DFS method. The values for even m (≥2) and odd m (≥3) in the new DFS method are similar to those for m=2 and m=3 in the SH method, respectively. In the SH method, when m is large, the values become close to zero at high latitudes.

When using the basis functions of Orszag in Eq. (2), we can also obtain results equivalent to the new DFS method by calculating the expansion coefficients using the least-squares method with Lagrange multipliers in order to minimize the error while satisfying the pole conditions in Eq. (4).

Figure 3a is the same as Fig. 2a except that N=63. In the old DFS method using Eq. (6), we set N=63 for m=0 and N=64 for m≠0. Because N=J for even m (≥2), the forward transform followed by the inverse transform does not change the initial values at the grid points, and the oscillations do not appear in the old DFS method. For this reason, Yoshimura and Matsumura (2005) and Yoshimura (2012) set N=J for even m (≥2) to improve stability. However, there is a problem with the latitudinal derivative in the old DFS method even when N=J for even m (≥2). Fig. 3b is the same as Fig. 3a except that it also shows the values between grid points calculated from the expansion coefficients by using Eq. (6) or Eq. (10). The large oscillations appear in the old DFS method with Grid [0], and it makes the latitudinal derivative at the grid points unrealistically large. In the new DFS method with the least-squares method, the large oscillations do not appear.

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f03

Figure 3(a) The same as Fig. 2a except for N=63. Panel (b) is the same as (a) except that the values between grid points calculated from the expansion coefficients are also shown.

Download

Table 1Normalized L2 errors of Laplacian operator calculation (2f). We use the old DFS method with Grid [0]; the new DFS methods with Grid [0], Grid [1], and Grid [-1]; and the SH method. J0 is the number of latitudinal grid points in Grid [0]. The truncation wavenumber N2J0/3.

Download Print Version | Download XLSX

Table 2The same as Table 1 except that the global mean values of calculated 2f are shown.

Download Print Version | Download XLSX

5 Tests of the DFS methods with the Laplacian operator and the Helmholtz equation

We examine the accuracy of the old and new DFS methods for the Laplacian operator in Eq. (65) and the Helmholtz equation

(91) 1 - ε 2 f = h .

Here, we give the function f as

(92) f = H 4 1 + cos π r R 2 if  r < R , 0 if  r R ,

(93) r = a cos - 1 sin ϕ c sin ϕ + cos ϕ c cos ϕ cos λ - λ c ,

where H=1000,R=a/3, ϕ is latitude, λ is longitude, a is the radius of the Earth, and r is the distance between (λ,ϕ) and the center λc,ϕc=3π/2,π/2-0.05. The function f is similar to the cosine bell in Williamson test case 1 (Williamson et al., 1992), but 1+cosπr/R is squared so that the second derivative of f is continuous. To easily calculate the exact values of 2f, the center is temporarily set to the North Pole, that is, λc,ϕc=0,π/2 and r=acos-1sinϕ=aθ, where θ is colatitude. At this time, g is calculated as follows:

(94) g = 2 f = 1 a 2 1 sin 2 θ 2 f λ 2 + 1 sin θ θ sin θ f θ = - cos θ sin θ H 2 a 2 π a R 1 + cos π r R sin π r R + H 2 a 2 π a R 2 sin 2 π r R - 1 + cos π r R cos π r R .

Equation (94) is satisfied at any position of the center. The function h in Eq. (91) is calculated by

(95) h = 1 - ε 2 f = f - ε g ,

where ε=0.01a2 and f and g are given by Eqs. (92) and (94).

To examine the accuracy for the Laplacian operator, f is given by Eq. (92), and 2f is calculated from f with the old DFS method (Cheong, 2000a), the new DFS method (see Sect. 3.7), and the SH method. The calculated values are compared with the exact values of 2f in Eq. (94). Here, the exact values of 2f are truncated by the forward transform followed by the inverse transform in order to see the error that does not include the error due to inability to resolve at the resolution. Table 1 shows the normalized L2 error between the calculated values and the exact values, which is normalized by the L2 norm of the exact values. The differences in error between the methods are small, but the results of the SH method are a little better than the old and new DFS methods. Table 2 shows the global mean of calculated 2f. The exact value of the global mean of 2f is zero. In Table 2, the global means calculated with each method are very close to zero. The global means of 2f in the DFS methods using Grid [1] and Grid [1] are not as close to zero as those in the DFS methods using Grid [0] and the SH method. This seems to be because the accuracy of the meridional discrete cosine and sine transforms in the DFS methods using Grid [1] and Grid [1] is not as good as that in the DFS methods using Grid [0].

To examine the accuracy of the solution of the Helmholtz equation, h is given in Eq. (95), and the Helmholtz equation in Eq. (91) is solved with the old DFS method (Cheong, 2000a), the new DFS method (see Sect. 3.8), and the SH method. The calculated values are compared with the exact solution f in Eq. (92). The exact values of f are also truncated as described above. Table 3 shows the normalized L2 error between the calculated values and the exact values. The differences in error between the methods are small, and which method is better depends on the resolution and the arrangement of the grid points.

Table 3The same as Table 1 except that L2 errors of the solution of the Helmholtz equation are shown.

Download Print Version | Download XLSX

6 Evaluation of the DFS methods using shallow-water test cases

We ran the Williamson test cases 1, 2, 5 ,and 6 (Williamson et al., 1992), and the Galewsky test case (Galewsky et al., 2004) in the model using the new DFS method described in Sect. 3, the model using the old DFS method of Yoshimura and Matsumura (2005), and the model using the SH method. By comparing the results of these model, we evaluated the old and new DFS methods.

6.1 Shallow-water equations on a sphere

The prognostic equations of the shallow-water model on a sphere are

(96)dvdt=-2Ω×vH-gh,(97)dh-hsdt=-h-hsv,

where t is time, v is the horizontal wind vector, h is the height, hs is the surface height, g is the acceleration due to gravity, Ω is the three-dimensional angular velocity of the Earth's rotation, and the subscript H indicates the horizontal component. Equation (96) is converted for the advective treatment of the Coriolis term (Temperton, 1997) into

(98) d v + 2 Ω × r d t = - g h ,

where r is the three-dimensional position vector from the Earth's center. Equation (97) is converted for the spatially averaged Eulerian treatment of mountains (Ritchie and Tanguay, 1996) into

(99) d h d t = - h - h s v + v h s .

Equations (98) and (99) are integrated in time using a two-time-level semi-implicit semi-Lagrangian scheme (see Appendix I).

6.2 Models

We ran the shallow-water test cases in the semi-implicit semi-Lagrangian shallow-water model or the Eulerian advection model (see Sect. 6.3) using the new DFS method (hereafter, the new DFS model). We also ran the same test cases in the model using the old DFS method of Yoshimura and Matsumura (2005) with the basis functions of Cheong (2000a, b) (hereafter, the old DFS model), and in the model using the SH method (hereafter, the SH model) for comparison. The new DFS model was run for Grid [0], [1], and [1]. In the old DFS model, Grid [0] is used. In the SH model, the Gaussian grid is used. We use a regular longitude–latitude grid and not a reduced grid. We use the time step Δt=3600 s at about 300 km resolution with J0=64, Δt=1800 s at about 120 km resolution with J0=160, Δt=1200 s at about 60 km resolution with J0=320, Δt=600 s at about 20 km resolution with J0=960, and Δt=90 s at about 1.3 km resolution with J0=15360, where J0 is the number of latitudinal grid points in Grid [0]. The number of latitudinal grid points J is J0 in Grid [0] (and in the Gaussian grid), J0+1 in Grid [1], and J0−1 in Grid [1] (see Sect. 2). The number of longitudinal grid points I is 2J0. The meridional truncation wavenumber N and the zonal wavenumber M are set to be equal. In the Eulerian advection model, shorter time steps are used as shown in Sect. 6.3. Horizontal diffusion is not used in all test cases. The zonal Fourier filter described in Appendix F is used in the DFS models. We have confirmed that numerical instability occurs in some test cases in the old DFS shallow-water model without the zonal Fourier filter, but stable integration is possible in all test cases shown here in the new DFS semi-Lagrangian shallow-water model, even without the zonal Fourier filter. In the new DFS Eulerian advection model, the zonal Fourier filter is necessary (see Sect. 6.3).

The zonal Fourier transforms in all the models and the meridional Fourier cosine and sine transforms in the DFS models are calculated using the Netlib BIHAR library, which includes a double-precision version of the Netlib FFTPACK library (Swarztrauber, 1982). The meridional Legendre transform in the SH model is calculated using the ISPACK library (Ishioka, 2018), which adopts on-the-fly computation of the associated Legendre functions. We use the ISPACK library's optimization option for Intel AVX512, which is highly optimized by using assembly language together with Fortran.

6.3 Williamson test case 1

Williamson test case 1 simulates a cosine bell advection. In the semi-Lagrangian models, the advection is calculated in the semi-Lagrangian scheme, and the horizontal derivatives calculated from the expansion coefficients are not used for the advection calculation. Therefore, we also use the Eulerian scheme here to simulate the advection in the DFS and SH models to test the expansion methods. The advection equation is

(100) d h d t = h t + v h .

In the Eulerian models, the advection term v⋅∇h is evaluated using the spectral transform method. The advection equation is integrated by the leapfrog scheme with the Robert–Asselin time filter (Robert, 1966; Asselin, 1972) with a coefficient of 0.05. The horizontal diffusion is not used, but the zonal Fourier filter is used in the old and new DFS methods. In Eq. (F1), the value M0=20 is used in the DFS semi-Lagrangian shallow-water models. However, the larger the value M0 is, the higher the longitudinal resolution around the pole becomes. Because of this, when the Eulerian scheme is used and M0 is large, a time step must be very short due to the CFL condition. Therefore M0 should be as small as possible. We have tested M0=0, but this degrades the result of Williamson test case 1. We have also tested M0=1, and this result is good. Therefore, we use M0=1 in the Eulerian models.

Figure 4 shows the predicted height after a 12 d integration in Williamson test case 1 when using the Eulerian advection models at the resolution J0=64. The meridional truncation wavenumber N and the zonal truncation wavenumber M are set as N=M=422J0/3 because the two-thirds rule (Orszag, 1971) is used in order to avoid aliasing in the nonlinear advection term. The time step is 30 min. The angle between the solid body rotation and the polar axis α is  π/2-0.05. The results for DFS [0], DFS [1], DFS [1], and SH are very similar. Instability occurs in the old DFS model without horizontal diffusion. This is probably because of the appearance of high-wavenumber oscillations due to the wavenumber truncation with N2J0/3 for even m (≠0) in the old DFS method, as shown in Sect. 4. Table 4 shows the normalized L2 errors of the predicted height after a 12 d integration when using the Eulerian advection models. The time steps are 30, 15, 7.5, and 2.5 min at the resolution J0=64, 160, 320, and 960 (N=42, 106, 213 and 639), respectively. The errors are very close between the models at each resolution. At the resolution N=639, the new DFS model without horizontal diffusion is unstable when the time step is 200 s. The SH model without horizontal diffusion is stable when the time step is 240 s and unstable when the time step is 300 s. One reason for this difference in time step is probably that the longitudinal resolution near the poles is higher in the new DFS model with M0=1 than in the SH model. When the fourth-order horizontal diffusion in Eq. (86) with ε=a4Δt/7.2×3600×1072N2 is used, the both new DFS and SH models are stable when the time step is 240 s and are unstable when the time step is 300 s. The old DFS model is unstable even when the same fourth-order horizontal diffusion is used. Higher-order horizontal diffusion, which effectively smooths out the high-wavenumber components, is necessary to stabilize the Eulerian old DFS model (Cheong, 2000b; Cheong et al., 2002).

Table 4Normalized L2 errors of the predicted height after a 12 d integration in Williamson test case 1 when using the Eulerian advection models. The truncation wavenumber N2J0/3.

Download Print Version | Download XLSX

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f04

Figure 4Predicted height (m) in the Eulerian models after a 12 d integration in Williamson test case 1: (a) the new DFS model with Grid [0], (b) the new DFS model with Grid [1], (c) the new DFS model with Grid [1], and (d) the SH model. The number of longitudinal (I) and latitudinal (J) grid points is shown in the form I×J. In the top row, the black contour shows the predicted height, and the red contour shows the reference solution. In the bottom row, color shading shows the difference between the predicted height and the reference solution.

Download

Table 5 shows the same information as Table 4 but using the semi-Lagrangian scheme. In the semi-Lagrangian models, the forward transform followed by the inverse transform are executed at every time step, but the expansion coefficients are not used for the advection calculation. The time steps are the same as described in Sect. 6.2. The errors are very close between the models. At the resolution J0= 64, the errors in the semi-Lagrangian models are larger than those in the Eulerian models, but at the resolutions J0=160, 320, and 960, the errors in the semi-Lagrangian models are smaller than those in the Eulerian models.

Table 5The same as Table 4 but using the semi-Lagrangian models and the truncation wavenumber NJ0-1.

Download Print Version | Download XLSX

The conservation of mass in Williamson test case 1 was also examined, and the results are shown in Sect. S2 in the Supplement.

6.4 Williamson test case 2

Williamson test case 2 simulates a steady-state nonlinear zonal geostrophic flow. In this test case, the angle between the solid body rotation and the polar axis α is given, and the zonal and meridional components of 2Ω×r become

(101) 2 Ω × r = 2 Ω a cos θ cos α + cos λ sin θ sin α , - 2 Ω a sin λ sin α .

Figure 5 shows the time series of forecast errors of the height for a 5 d integration in Williamson test case 2 with α=π/2-0.05 in the models at the resolution J0=64 and N=63 (DFS) or N=62 (SH) using no horizontal diffusion. The normalized L1, L2, and L errors are almost the same between the new DFS models using Grids [0], [1], and [1]; the old DFS model; and the SH model. Table 6 shows the normalized L2 errors of the predicted height after a 5 d integration. The errors are almost the same between the old DFS, new DFS, and SH models at each resolution.

Table 6The same as Table 5 but for the errors after a 5 d integration in Williamson test case 2.

Download Print Version | Download XLSX

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f05

Figure 5Time series of prediction error of height (m) for 5 d (120 h) integration in Williamson test case 2 α=π/2-0.05. The number of longitudinal grid points is I=128. The number of latitudinal grid points in Grid [0] is J0=64. The truncation wavenumber is N=63. Solid, dashed, and dotted lines represent normalized L1, L2, and L errors, respectively. The colors blue, green, red, purple, and orange represent the models using SH, old DFS with Grid [0], new DFS with Grid [0], new DFS with Grid [1], and new DFS with Grid [1], respectively.

Download

The conservation of mass, energy, and vorticity in the Williamson test cases 2, 5, and 6 was also examined, and the results are shown in Sect. S2 in the Supplement.

6.5 Williamson test case 5

Williamson test case 5 simulates zonal flow over an isolated mountain. Figure 6 shows the predicted height after a 15 d integration in Williamson test case 5 with h0=5960 m. The result of the high-resolution SH model at the resolution J=960 and N=958 is regarded as the reference solution. Horizontal diffusion is not used. The errors with respect to the reference solution are almost the same for the new DFS models, the old DFS model, and the SH model at the resolution J0=64. Table 7 shows the normalized L2 errors of the predicted height after a 15 d integration. The errors are almost the same between the old DFS, new DFS, and SH models at each resolution. The errors do not decrease when the resolution increases, which is different from the results in the other test cases. This may be because the mountain topography is not a differentiable function, and the mountain is added impulsively on to an initially balanced flow (Galewsky et al., 2004).

Table 7The same as Table 5 but for the errors after a 15 d integration in Williamson test case 5. The result of the high-resolution SH model with J0=960 and N=958 is regarded as the reference solution.

Download Print Version | Download XLSX

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f06

Figure 6Predicted height (m) after a 15 d integration in Williamson test case 5: (a) the new DFS model with Grid [0], (b) the new DFS model with Grid [1], (c) the new DFS model with Grid [1], (d) the old DFS model with Grid [0], (e) the SH model, and (f) the SH model at high resolution (which is regarded as the reference solution). The number of longitudinal (I) and latitudinal (J) grid points is shown in the form I×J. N is the truncation wavenumber. Color shading shows the error with respect to the reference solution.

Download

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f07

Figure 7Longitudinal distributions of meridional wind (m s−1) at the grid points near the South Pole after a 15 d integration in Williamson test case 5. Results of the models using Grid [0] with (a) I=128, J0=64, and N=63 and (b) I=1920, J0=960, and N=959. Green (red) lines represent the old (new) DFS models.

Download

Figure 7 shows the longitudinal distributions of meridional wind at the grid points near the South Pole after a 15 d integration in the old and new DFS models using Grid [0] at the resolutions J0=64 and J0=960. While the zonal wavenumber 1 component is dominant in the new DFS model at the resolution J0=64 and N=63, high zonal wavenumber noise appears in the old DFS model at the same resolution. One possible reason is that the latitudinal derivative at the grid points can be unrealistically large in the old DFS method, even when N=J0 for even m (≥2), as described in Sect. 4 (Fig. 3b). The new DFS expansion method with the least-squares method does not have this problem. By using the new expansion method with the least-squares method, the high zonal wavenumber noise does not appear even in the model that uses the same DFS basis functions as in Eq. (11) except that the basis function for odd m (≥3) is sin θcos nθ instead of sin 2θsin nθ. In the old DFS model at high resolution with J0=960 and N=959, the high-wavenumber noise is not seen in Fig. 7. The higher the resolution, the smaller the high-wavenumber noise becomes.

Figure 8 shows the kinetic energy spectra of the horizontal winds (Lambert, 1984) after a 15 d integration in Williamson test case 5. The kinetic energy spectra in the DFS models are calculated from the SH expansion coefficients, which are obtained by first calculating the Gaussian grid point values from the DFS coefficients using Eq. (10) for the new DFS method and Eq. (6) for the old DFS method and then calculating the SH expansion coefficients from the Gaussian grid point values by using a forward Legendre transform. In the old DFS model with J0=64 and N=63, the high-wavenumber components are larger than in the other models, which is related to the high-wavenumber noise near the South Pole in Fig. 7. In the old DFS model with J0=960, the high-wavenumber components are a little larger than in the other models, but the differences are slight.

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f08

Figure 8Kinetic energy spectrum of horizontal winds (m2 s−2) after a 15 d integration in Williamson test case 5. Results of the models with (a) I=128, J0=64, and N=63 (DFS) or N=62 (SH) and (b) I=1920, J0=960, and N=959 (DFS) or 958 (SH). The colors blue, green, red, purple, and orange represent the models using SH, old DFS with Grid [0], new DFS with Grid [0], new DFS with Grid [1], and new DFS with Grid [1], respectively.

Download

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f09

Figure 9The same as Fig. 6 but with the truncation wavenumber N.

Download

Figure 9 shows the predicted height after a 15 d integration in Williamson test case 5, which is the same as Fig. 6 but using the truncation wavenumber N=422J0/3. In our semi-implicit semi-Lagrangian models, we usually useN satisfying NJ0-1, which is called linear truncation. However, here N is determined to satisfy N2J0/3 to eliminate aliasing errors with quadratic nonlinearity, which is called quadratic truncation. When using the quadratic truncation, the new DFS models with Grids [0], [1], and [1] are stable without horizontal diffusion, but the old DFS model without strong high-order horizontal diffusion is unstable. The numerical instability in the old DFS model probably occurs because of the high-wavenumber oscillations due to the quadratic wavenumber truncation for even m (≠0) (see Sect. 4), as is also the case in Williamson test case 1 with the Eulerian model. The results of the new DFS models are almost the same as those of the SH model. Table 8 is the same as Table 7 except for N2J0/3. The results of the new DFS models and the SH model with N2J0/3 in Table 8 are very similar to those with NJ0-1 in Table 7 when J0 is the same.

Table 8The same as Table 7 but with N2J0/3.

Download Print Version | Download XLSX

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f10

Figure 10The same as Fig. 8 but with the truncation wavenumber N.

Download

Figure 10 shows the kinetic energy spectrum of the horizontal winds after a 15 d integration in Williamson test case 5, which is the same as Fig. 8 except for the truncation wavenumber N2J0/3. At the resolution J0=64 and N=42, the high-wavenumber components are a little larger in the SH model than in the new DFS model. At the resolution J0=960 and N=639, very small oscillations appear in the high-wavenumber region in the SH model but not in the new DFS models. In the SH model, the wind components u and v divided by sin θ are transformed from grid space to spectral space (Ritchie, 1988; Temperton, 1991), which seems to reduce the accuracy and cause the small oscillations in the high-wavenumber region. Another way to transform u and v from grid space to spectral space in the SH model is to use the vector harmonic transform (see Sect. 3.6). Both ways are algebraically equivalent (Temperton, 1991). However, the latter way avoids dividing u and v by sin θ and provides remarkable stability and accuracy (Swarztrauber, 2004). The latter way is similar to the new DFS expansion method for u and v using the least-squares method described in Sect. 3.6 and probably eliminates the small oscillations in the SH model. Alternatively, using D and ζ instead of u and v as prognostic variables may eliminate the small oscillations.

6.6 Williamson test case 6

Figure 11 shows the predicted height after a 14 d integration in Williamson test case 6. The error is similar between the old and new DFS models using Grid [0] and the SH model. The error in the new DFS model using Grid [1] is the smallest. This is probably because Grid [1] has grid points at the poles, where the minimum height exists, and on the Equator, where the maximum height exists. The error in the new DFS model using Grid [1] is the second smallest. This is probably because Grid [1] has grid points on the Equator. Table 9 shows the normalized L2 errors of the predicted height after a 14 d integration. The error in the new DFS model using Grid [1] is the smallest, and the error in the new DFS model using Grid [1] is the second smallest at each resolution. The errors in the old and new DFS models using Grid [0] and in the SH model are very close.

Table 9The same as Table 7 but the errors after a 14 d integration in Williamson test case 6.

Download Print Version | Download XLSX

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f11

Figure 11The same as Fig. 6 but for predicted height (m) after a 14 d integration in Williamson test case 6.

Download

6.7 Galewsky test case

The Galewsky test case simulates a barotropically unstable midlatitude jet. Figure 12 shows the predicted vorticity after a 6 d integration in the Galewsky test case for the models at 1.3 km resolution with J0=15360 and the quadratic truncation N=10 239, without horizontal diffusion. The result in the new DFS model using Grid [0] is almost the same as in the SH model. The old DFS model is unstable for the same reason as that shown in Sect. 6.5 (Fig. 9). Figure 13 shows the kinetic energy spectrum of horizontal winds after a 6 d integration in the Galewsky test case. The results are almost the same for the DFS models using Grid [0], [1], and [1] and the SH model, but very small oscillations appear near the truncation wavenumber in the SH model. This is probably for the same reason as in Williamson test case 5 in Fig. 10.

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f12

Figure 12Predicted vorticity (s−1) after a 6 d integration in the Galewsky test case. (a) The new DFS model with Grid [0]. (b) The SH model at 1.3 km resolution with I=30 720, J0=15 360, and N=10 239.

Download

The results of the Galewsky-like test case using the north–south symmetric initial conditions are shown in Sect. S3 in the Supplement.

6.8 Elapsed time

Figure 14 shows the elapsed time for the 15 d integration in Williamson test case 5 in the SH model and the new DFS model using Grid [0] at 20 km resolution with J0=960 and N=958 (SH) or N=959 (DFS) and that for the 6 d integration in the Galewsky test case at 1.3 km resolution with J0=15 360 and N=10 239. We use one node (with two Intel Xeon Gold 6248 CPUs with 20 cores per CPU) of the FUJITSU Server PRIMERGY CX2550 M5 in the MRI. The source code written in Fortran is compiled with the Intel compiler. OpenMP (Open Multi-Processing) parallelization is used, but MPI (Message Passing Interface) parallelization is not used. The elapsed time in the SH model is larger than in the DFS model, although the Legendre transform used in the SH model is highly optimized for Intel AVX512. The higher the resolution, the larger the difference of the elapsed time between the models becomes. This is because the Legendre transform used in the SH model requires O(N3) operations, whereas the Fourier cosine and sine transforms used in the DFS model require only O(N2log N) operations. If the fast Legendre transform, which requires only N2(log N)3 operation, is used instead of the usual Legendre transform in the SH model, the difference in the elapsed time between the models will be reduced at high resolutions. We have not tested the fast Legendre transform yet because we do not have subroutines for the fast Legendre transform.

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f13

Figure 13Kinetic energy spectrum of horizontal winds (m2 s−2) after a 6 d integration in the Galewsky test case. (a) Results of the models with I=30 720, J0=15 360, and N=10 239. The colors blue, red, purple, and orange represent the models using SH, DFS with Grid [0], DFS with Grid [1], and DFS with Grid [1], respectively. Panel (b) is the same as (a) but shows the high-wavenumber region.

Download

https://gmd.copernicus.org/articles/15/2561/2022/gmd-15-2561-2022-f14

Figure 14Elapsed time (s) for (a) 15 d integration in Williamson test case 5 in the SH model and the new DFS model at 20 km resolution with I=1920, J0=960, and N=959. (b) The 6 d integration in the Galewsky test case at 1.3 km resolution with I=30 720, J0=15 36, and N=10 239. There is no monitoring output during elapsed time measurement.

Download

7 Conclusions and perspectives

We have developed the new DFS method to improve the numerical stability of the DFS model, which has the following two improvements.

  1. A new expansion method with the least-squares method is used to calculate the expansion coefficients so that the error due to the meridional wavenumber truncation is minimized. The method also avoids dividing by sin θ before taking the forward Fourier cosine or sine transform.

  2. New DFS basis functions are used, which guarantees that not only scalar variables but also vector variables and the gradient of scalar variables are continuous at the poles.

The equations obtained with the least-squares method are equivalent to those obtained with the Galerkin method. We also use the Galerkin method to solve partial differential equations such as the Poisson equation, the Helmholtz equation, and the shallow-water equations.

To test the new DFS method, we conducted experiments for Williamson test cases 2, 5, and 6 and the Galewsky test case in the semi-implicit semi-Lagrangian shallow-water models using the new DFS method with three types of equally spaced latitudinal grids with or without the poles. We also ran Williamson test case 1, which simulates a cosine bell advection, in the Eulerian and semi-Lagrangian advection models. We compared the results between the new DFS models using the new DFS method, the old DFS model using the method of Yoshimura and Matsumura (2005) with the basis functions of Cheong (2000a, b), and the SH model.

The high zonal wavenumber noise of the meridional wind appears near the poles in the old DFS model, but not in the new DFS models in Williamson test case 5. One possible reason is that the latitudinal derivative at the grid points can be unrealistically large in the old DFS method even when the truncation wavenumber N for even m (≠0) is equal to the number of latitudinal grid points J, while the new DFS expansion method with the least-squares method does not have this problem. In the old DFS model, N<J for even m (≠0) causes numerical instability. In the new DFS model, an arbitrary meridional wavenumber truncation N (<J) can be used without the stability problem because the error due to meridional wavenumber truncation is small when using the new DFS expansion method with the least-squares method. This is one of the merits of the new DFS method because quadratic truncation N2J/3 or cubic truncation NJ/2 are usually used in the Eulerian model and are also starting to be used in the semi-Lagrangian model instead of the linear truncation NJ-1 for stability and efficiency at high resolutions (Wedi, 2014; Hotta and Ujiie, 2018; Dueben et al., 2020). We have also confirmed that in the new DFS model, stable integration is possible in all test cases shown here even without using the zonal Fourier filter (unlike in the old DFS model). Thus, the numerical stability of the semi-implicit semi-Lagrangian model using the new DFS method is very good. In Williamson test case 1, the Eulerian advection model using the new DFS method also gives stable results without horizontal diffusion but with a zonal Fourier filter. The Eulerian advection model using the old DFS method is unstable without horizontal diffusion or with weak fourth-order horizontal diffusion. In the old DFS model, the use of the semi-Lagrangian scheme is important for numerical stability. On the other hand, the advection model using the new DFS method is stable, even when the Eulerian scheme is used instead of the semi-Lagrangian scheme.

The results of the new DFS model are almost the same as the SH model. However, in the SH shallow-water model without horizontal diffusion, very small oscillations appear in the high-wavenumber region of the kinetic energy spectrum in some cases, unlike in the new DFS model. This seems to be because the wind components u and v divided by sin θ are transformed from grid space to spectral space in the SH model. The small oscillations with the SH model can probably be eliminated by using the vector harmonic transform, which is similar to the new DFS expansion method for u and v using the least-squares method and avoids dividing u and v by sin θ. Alternatively, using divergence and vorticity instead of u and v as prognostic variables may eliminate the small oscillations.

The elapsed time in the new DFS model is shorter than in the SH model (especially at high resolution) because the Fourier transform requires only O(N2log N) operations, while the Legendre transform in the SH model requires O(N3) operations. We have executed our shallow-water models on Intel CPUs. The execution on GPUs is one important topic, but we have not tested our models in this way because execution on GPUs is not an easy task. MPI parallelization is another important topic. However, in our shallow-water models, we use only OpenMP parallelization and not MPI parallelization due to the simplicity of the source code.

We developed hydrostatic and nonhydrostatic global atmospheric models using the old DFS method (Yoshimura and Matsumura, 2005; Yoshimura, 2012) and conducted typhoon prediction experiments in the nonhydrostatic global atmospheric model using the old DFS method in the Global 7 km mesh nonhydrostatic Model Intercomparison Project for improving TYphoon forecast (TYMIP-G7; Nakano et al., 2017). We have already developed a nonhydrostatic (or hydrostatic) atmospheric model using the new DFS method where both OpenMP and MPI parallelization are used. We will describe the nonhydrostatic DFS model and the MPI parallelization in another paper after improving the nonhydrostatic dynamical core as needed.

Appendix A: Trigonometric identities

Here we list the trigonometric identities used in transforming the expressions in this paper. The following identities are satisfied:

(A1a)sinnθcosnθ=12sinn+nθ+sinn-nθ,(A1b)cosnθsinnθ=12sinn+nθ-sinn-nθ,(A1c)cosnθcosnθ=12cosn+nθ+cosn-nθ,(A1d)sinnθsinnθ=12-cosn+nθ+cosn-nθ.

From Eq. (A1), the following identities are derived:

(A2a)sinθcosnθ=12sinn+1θ-sinn-1θ,(A2b)sinθsinnθ=12-cosn+1θ+cosn-1θ,(A2c)sin2θsinnθ=14-sinn-2θ+2sinnθ-sinn+2θ,(A2d)sin2θcosnθ=14-cosn-2θ+2cosnθ-cosn+2θ,(A2e)sin4θsinnθ=116sinn-4θ-4sinn-2θ+6sinnθ-4sinn+2θ+sinn+4θ.

From Eq. (A1), the following orthogonal relations in longitude are derived:

(A3a)02πcosmλcosmλdλ=2π form=m=0,π form=m0,0 formm,(A3b)02πcosmλsinmλdλ=0,(A3c)02πsinmλsinmλdλ=π for m=m0,0 for mm.

Similarly, from Eq. (A1), the following orthogonal relations in latitude are derived:

(A4a)0πcosnθcosnθdθ=π forn=n=0,12π forn=n0,0 fornn,(A4b)0πsinnθsinnθdθ=12π for n=n0,0 for nn.

By using Eqs. (A1) and (A2), the following relations are derived:

(A5a)θsinlθcosnθ=n+l2sinl-1θcosn+1θ-n-l2sinl-1θcosn-1θ,(A5b)sinθθsinθθsinlθcosnθ=n+ln+l+14sinlθcosn+2θ-2n2-2l2+2l4sinlθcosnθ+n-ln-l-14sinlθcosn-2θ,(A5c)θsinlθsinnθ=n+l2sinl-1θsinn+1θ-n-l2sinl-1θsinn-1θ,(A5d)sinθθsinθθsinlθsinnθ=n+ln+l+14sinlθsinn+2θ-2n2-2l2+2l4sinlθsinnθ+n-ln-l-14sinlθsinn-2θ.
Appendix B: Discrete Fourier cosine and sine transforms in latitude

Forward discrete Fourier cosine and sine transforms are performed in Eqs. (23) and (57), and inverse discrete Fourier cosine and sine transforms are performed in Eqs. (13) and (52) in the latitudinal direction. The calculation of the discrete cosine and sine transforms in Grids [0], [1], and [1] is shown below. Here, g(θj) and h(θj) are grid point values, θj is the colatitude defined in Eq. (7), and gn and hn are expansion coefficients.

When using Grid [0], inverse and forward discrete cosine transforms are performed as follows:

(B1a)gθj=n=0J0-1gncosnθj,(B1b)gn=bJ0j=0J0-1gθjcosnθj,b1 for n=02 for 1nJ0-1.

When using Grid [0], inverse and forward discrete sine transforms are performed as follows:

(B2a)hθj=n=1J0hnsinnθj,(B2b)hn=bJ0j=0J0-1hθjsinnθj,b1 forn=J02 for1nJ0-1.

When using Grid [1], inverse and forward discrete cosine transforms are performed as follows:

(B3a)gθj=n=0J0gncosnθj,(B3b)gn=bJ0j=0J0cgθjcosnθj,b1 forn=0,J02 for1nJ0-1,c1/2 forj=0,J01 for1jJ0-1.

When using Grid [1], inverse and forward discrete sine transforms are performed as follows:

(B4a)hθj=n=1J0-1hnsinnθj,hθ0=hθJ0=0,(B4b)hn=2J0j=1J0-1hθjsinnθj1nJ0-1.

Grid [1] is the same as Grid [1], except that there are no grid points at the North and South poles. The zonal wavenumber components of scalar variables at the poles are zero except for m=0 (see Eqs. 10 and 11), and those of vector variables at the poles are zero except for m=1 (see Eq. 52). When we use Grid [1] and the values at the poles are known to be zero, forward and inverse discrete cosine transforms can be performed using Eq. (B3), and forward and inverse discrete sine transforms can be performed using Eq. (B4) in the same way as for Grid [1]. When we use Grid [1] and the values at the poles are unknown (i.e., the zonal wavenumber components of scalar variables for m=0, and those of vector variables for m=1), the inverse discrete cosine transform can be performed like Eq. (B3a) as follows:

(B5) g θ j = n = 0 J 0 - 2 g n cos n θ j ,

where n is from 0 to J0-2=J-1 because the number of the meridional grid points is J0-1=J in Grid [1]. However, the forward discrete cosine transform cannot be performed like Eq. (B3b). We can calculate the expansion coefficients gn from g(θj) in the following way. Equation (B5) is multiplied by sin θj, and we define g^θj as follows:

(B6) g ^ θ j g θ j sin θ j = n = 0 J 0 - 2 g n sin θ j cos n θ j .

We can expand g^θj as follows:

(B7) g ^ θ j = n = 1 J 0 - 1 g ^ n sin n θ j .

The expansion coefficients g^n can be obtained from g^θj in the same way as in Eq. (B4b) by forward discrete sine transform:

(B8) g ^ n = 2 J 0 j = 1 J 0 - 1 g ^ θ j sin n θ j .

From Eqs. (B6) and (B7), we obtain

(B9) n = 0 J 0 - 2 g n sin θ cos n θ = n = 1 J 0 - 1 g ^ n sin n θ .

By using Eq. (A2a), we obtain

(B10) n = 0 J 0 - 2 g n sin θ cos n θ = g 0 - g 2 2 sin θ + n = 2 J 0 - 3 g n - 1 2 - g n + 1 2 sin n θ + g J 0 - 3 2 sin J 0 - 2 θ + g J 0 - 2 2 sin J 0 - 1 θ .

By substituting Eq. (B10) into Eq. (B9) and comparing the left and right sides of the equation, we obtain

(B11) g ^ n = g 0 - g 2 2  for n = 1 , g n - 1 2 - g n + 1 2  for 2 n J 0 - 3 , g J 0 - 3 2  for n = J 0 - 2 , g J 0 - 2 2  for n = J 0 - 1 .

We can calculate g^θj from g(θj) using Eq. (B6), calculate g^n from g^θj using Eq. (B8), and calculate gn from g^n using Eq. (B11).

Appendix C: The upper limit of the meridional truncation wavenumber N

In the new DFS method, the meridional truncation wavenumber N is used for the new DFS meridional basis functions in Eq. (12) and for the discrete cosine or sine transform of a scalar variable (Eqs. 13 and 23), derivatives of a scalar variable (Eqs. 18 and 20), and a wind vector (Eqs. 52 and 57). In Grid [0], the upper limit of N is J0−1 for each m because the discrete cosine transform in Eq. (B1), where the maximum value of n is J0−1, is used for a scalar variable when m is even, and for vector components when m is odd. In Grid [1], the upper limit of N is J0−1 for each m because the discrete sine transform in Eq. (B4), where the maximum value of n is J0−1, is used for a scalar variable when m is odd, and for vector components when m is even. In Grid [1], the upper limit of N is J0−1 for m≥2 for the same reason as in Grid [1]. However, for m= 0 or 1 in Grid [1], the upper limit of N is J0−2 because the discrete cosine transform in Eq. (B5), where the maximum value of n is J0−2, is used for a scalar variable when m=0 and for vector components when m=1. Thus, the upper limit of N is J0−1, except that the upper limit of N for m= 0 or 1 in Grid [1] is J0−2. For example, in the model using the new DFS method with Grid [1] at the resolution J0=64 and N=63, we set N=63 for m≥2 but N=62 for m= 0 or 1.

Appendix D: Equations for the derivation of Eqs. (29) and (62)

T̃n,mcs in Eq. (23) is calculated by the forward Fourier cosine or sine transform as follows:

(D1) T ̃ n , m c s = b π 0 π cos n θ T m c s θ d θ , b 1 for n = 0 2 for n 0 , for  even   m , 2 π 0 π sin n θ T m c s θ d θ for  odd m .

The equations for the forward discrete Fourier cosine or sine transform are described in Appendix B. From Eq. (23) and (A4), we derive the following equations:

(D2a)bπ0πcosnθT̃mcs,Nθdθ=T̃n,mcsn=0,,Nfor  evenm(D2b)2π0πsinnθT̃mcs,Nθdθ=T̃n,mcsn=1,,Nfor  oddm.

From Eqs. (D1) and (D2), the following equations are satisfied:

(D3a)0πcosnθTmcsθdθ=0πcosnθT̃mcs,Nθdθn=0,,Nfor  evenm,(D3b)0πsinnθTmcsθdθ=0πsinnθT̃mcs,Nθdθn=1,,Nfor  oddm.

From Eqs. (D3), (11), (12), and (A2a)–(A2c), we derive

(D4) 0 π S n , m θ T m c s θ d θ = 0 π S n , m θ T ̃ m c s , N θ d θ n = N min,   m , , N max, m .

From Eqs. (28) and (D4), we derive Eq. (29).

We can also derive the following equations from Eq. (57) in the similar way to the derivation of (D3):

(D5a)0πsinnθumcsθdθ=0πsinnθũmcs,Nθdθn=1,,Nfor  evenm,(D5b)0πcosnθumcsθdθ=0πcosnθũmcs,Nθdθn=0,,Nfor  oddm

From Eqs. (D5), (11), (12), and (A2a)–(A2c), we derive

(D6a)0πmSn,mθsinθumcsθdθ=0πmSn,mθsinθũmcs,Nθdθn=Nmin,  m,,Nmax,m,(D6b)0πSn,mθθumcsθdθ=0πSn,mθθũmcs,Nθdθn=Nmin,  m,,Nmax,m.

We can also derive the same equations as Eq. (D6) except that u is replaced with v. Equation (D6) are used to derive Eq. (62).

Appendix E: Derivation of Eq. (30) from Eq. (29)

Here we derive Eq. (30d) for odd (m≥3) from Eq. (29). Equations (30b) and (30c) can also be derived similarly. By using Eqs. (10), (11), (23), (A2c), and (A2e), Eq. (29) is converted as follows:

(E1)l.h.s  of  Eq.  29  for  odd  m3=0πSn,mθTmc,Nθdθ=0πsin2θsinnθn=1N-2Tn,mcsin2θsinnθdθ=0πsinnθn=1N-2Tn,mc16sinn-4θ-4sinn-2θ+6sinnθ-4sinn+2θ+sinn+4θdθ=0πsinnθ10T1,mc-5T3,mc+T5,mc16sinθ+5T2,mc-4T4,mc+T6,mc16sin2θ+-5T1,mc+6T3,mc-4T5,mc+T7,mc16sin3θ+n=4N+2Tn-4,mc-4Tn-2,mc+6Tn,mc-4Tn+2,mc+Tn+4,mc16sinnθdθ(E2)r.h.s  of  Eq.  (29)  for  odd  m3=0πSn,mθT̃mcs,Nθdθ=0πsin2θsinnθn=1NT̃n,mcsinnθdθ,=0πsinnθn=1NT̃n,mc4-sinn-2θ+2sinnθ-sinn+2θdθ=0πsinnθ3T̃1,mc-T̃3,mc4sinθ+2T̃2,mc-T̃4,mc4sin2θ+n=3N+2-T̃n-2,mc+2T̃n,mc-T̃n+2,mc4sinnθdθ,

where 1nN-2. From Eqs. (29), (E1), (E2), and (A4b), Eq. (30d) is derived.

Appendix F: Zonal Fourier filter

In a regular longitude–latitude grid, the longitudinal grid spacing becomes narrow at high latitudes. In DFS methods, the zonal Fourier filter (Merilees, 1974; Boer and Steinberg, 1975; Cheong, 2000a), which filters out the high zonal wavenumber components at high latitudes, is usually used to obtain a more uniform resolution. The use of a reduced grid (Hortal and Simmons, 1991; Juang, 2004; Miyamoto, 2006; Malardel, 2016) has a similar effect to the zonal Fourier filter. In our atmospheric model using the old DFS method (Yoshimura, 2012), we use the reduced grid of Miyamoto (2006).

In this study, we use the regular longitude–latitude grid with the zonal Fourier filter, not the reduced grid, for the simplicity of its source code. We set the largest zonal wavenumber Mf at each colatitude θj as follows:

(F1) M f θ j = min M , M 0 + M sin θ j .

The values of Tmcθj and Tmsθj in Eq. (8) are set to zero for m>Mf(θj) during the spectral transform. We use the value M0=20 in the DFS shallow-water model to make the resolution similar to that in the reduced grid of Miyamoto (2006). In the DFS Eulerian advection model, we use the value M0=1 as described in Sect. 6.3.

Appendix G: Calculation of global mean and latitudinal area weight

The global mean value of TN,Mλ,θ in Eq. (15) can be calculated in spectral space by the following equation (Cheong 2000a):

(G1) G = 1 4 π 0 2 π 0 π m = 0 M T m c , N θ cos m λ + m = 1 M T m s , N θ sin m λ sin θ d θ d λ = 1 2 0 π n = 0 N T n , m = 0 c cos n θ sin θ d θ = n = 0 when n   is  even N T n , m = 0 c 1 - n 2 ,

where Eq. (A2a) is used.

The latitudinal area weight at each colatitude θj is calculated as follows.

  1. The latitudinal distribution of Tm=0cjθk for each j is given as

    (G2) T m = 0 c j θ k = 1   for k = j 0   for k j ,

    where 0kJ-1 in Grid [0] and Grid [1] and 1kJ in Grid [1] (see Sect. 2).

  2. From Tm=0cjθk, the meridional expansion coefficients Tn,m=0cj0nN are calculated by forward discrete cosine transform as described in Appendix B.

  3. The value of G calculated from Tn,m=0cj using Eq. (G1) is considered the latitudinal area weight w(θj) at colatitude θj.

    In Grid [0] and Grid[1], the distribution of w(θj) is smooth. However, in Grid [1], the distribution of w(θj) is not smooth because of the irregularity with Grid [1] (see Eqs. B5–B11 in Appendix B).

The latitudinal area weight w(θj) is used, for example, to calculate the global mean in the grid space.

Appendix H: Derivation of Eq. (63) from Eq. (62)

Here we describe the derivation of Eq. (63d) for odd m (≥3) from Eq. (62a). Equations (63b) and (63c) can also be derived similarly. By using Eqs. (52), (57), (11), (A2b), (A2c), and the same equations as Eq. (53), except that un,mc, χn,ms, and ψn,mc are replaced with un,ms, -χn,mc, and ψn,ms, respectively, and the same equations as Eq. (53), except that un,mc, χn,ms, and ψn,mc are replaced with vn,mc, ψn,ms, and -χn,mc, respectively, Eq. (62a) is converted into the following equation for odd m≥3:

(H1) 0 π - m - cos n + 1 θ + cos n - 1 θ 2 - m χ 1 , m c 2 a - u ̃ 0 , m s + - m χ 2 , m c 2 a + 3 ψ 1 , m s - ψ 3 , m s 4 a - u ̃ 1 , m s cos θ + n = 2 N m χ n - 1 , m c - χ n + 1 , m c 2 a + n - ψ n - 2 , m s + 2 ψ n , m s - ψ n + 2 , m s 4 a - u ̃ n , m s cos n θ - - n - 2 4 cos n - 2 θ + 2 n 4 cos n θ - n + 2 4 cos n + 2 θ m ψ 1 , m s 2 a - v ̃ 0 , m c + m ψ 2 , m s 2 a + - 3 χ 1 , m c + χ 3 , m c 4 a - v ̃ 1 , m c cos θ + n = 2 N m - ψ n - 1 , m s + ψ n + 1 , m s 2 a + n χ n - 2 , m c - 2 χ n , m c + χ n + 2 , m c 4 a - v ̃ n , m c cos n θ d θ = 0 .

When n≥4, by using Eq. (A4a), Eq. (H1) can be converted into

(H2) 0 π m 2 cos n + 1 θ m χ n , m c - χ n + 2 , m c 2 a + n + 1 - ψ n - 1 , m s + 2 ψ n + 1 , m s - ψ n + 3 , m s 4 a - u ̃ n + 1 , m s cos n + 1 θ - m 2 cos n - 1 θ m χ n - 2 , m c - χ n , m c 2 a + n - 1 - ψ n - 3 , m s + 2 ψ n - 1 , m s - ψ n + 1 , m s 4 a - u ̃ n - 1 , m s cos n - 1 θ + n - 2 4 cos n - 2 θ m - ψ n - 3 , m s + ψ n - 1 , m s 2 a + n - 2 χ n - 4 , m c - 2 χ n - 2 , m c + χ n , m c 4 a - v ̃ n - 2 , m c cos n - 2 θ - 2 n 4 cos n θ m - ψ n - 1 , m s + ψ n + 1 , m s 2 a + n χ n - 2 , m c - 2 χ n , m c + χ n + 2 , m c 4 a - v ̃ n , m c cos n θ + n + 2 4 cos n + 2 θ m - ψ n + 1 , m s + ψ n + 3 , m s 2 a + n + 2 χ n , m c - 2 χ n + 2 , m c + χ n + 4 , m c 4 a - v ̃ n + 2 , m c cos n + 2 θ d θ = 0 .

From Eqs. (H2) and (A4a), Eq. (63d) for n≥4 is derived. Equation (63d) for n≤3 can also be derived from Eqs. (H1) and (A4a).

Appendix I: Two-time-level semi-implicit semi-Lagrangian scheme for time integration

A two-time-level semi-implicit semi-Lagrangian scheme (e.g., Temperton et al., 2001) and the Stable Extrapolation Two-Time-Level Scheme (SETTLS; Hortal, 2002) are adopted to discretize the shallow-water equations in Eqs. (98) and (99) in time as follows:

(I1)v+2Ω×r-v+2Ω×rD0Δt=-ghD++h02+βvghD++h02-βvghD0+h+2,(I2)h+-hD0Δt=-h-hsDD++h-hsD02+vhsD++vhs02+βhh¯DD++h¯D02-βhh¯DD0+h¯D+2,

where

(I3) D v = 1 a 1 cos ϕ u λ + 1 cos ϕ v cos ϕ ϕ ,

is horizontal divergence; Δt is a time step; the superscripts , 0, and + mean past time (t−Δt), present time (t), and future time (tt), respectively; the superscript (+) means future time (tt) extrapolated in time, for example, h+=2h0-h-; the subscript D means the departure point; the absence of the subscript D means the arrival point; h¯ is a constant value of height for semi-implicit linear terms; βv and βh are second-order decentering parameters (Yukimoto et al., 2011). Using βv and βh larger than 1.0 (e.g., 1.2) increases the effect of the semi-implicit scheme improving computational stability, but βv=βh=1.0 is used here because h¯ larger than h is enough for stable calculations in the shallow-water model. The departure point xD is the upstream horizontal position from the arrival point x along the wind vector between present time (t) and future time (tt). Here, the arrival point x is on a grid point, and the departure point xD is not generally on a grid point. Since the right-hand sides of Eqs. (I1) and (I2) are the time average between present time (t) and future time (tt) and the spatial average between the departure point and the arrival point, these equations have second-order precision in time and space. In SETTLS, xD is calculated using

(I4) x D = x - v D + + v 0 2 Δ t .

However, when Δt is longer than 30 min, using vD+ extrapolated in time to calculate xD causes numerical instability in our experiments. To avoid instability when Δt is 1 h, here we use

(I5a)xD=x-vD0+v+2Δt,(I5b)v+vD0+2Ω×rD-2Ω×r-ghD++h02Δt,

instead of Eq. (I4), where v+ is a provisional future value obtained by discretizing Eq. (98) in an explicit semi-Lagrangian scheme. From Eq. (I5), we obtain

(I6) x D = x - Δ t v 0 + Ω × r - g Δ t h + 4 D - Ω × r - g Δ t h 0 4 .

This method using a provisional future value to calculate xD is similar to the method in Gospodinov et al. (2001). Since the value with the subscript D depends on xD, xD is calculated iteratively from Eq. (I6) (e.g., Ritchie et al., 1995; Temperton et al., 2001). Since xD is not generally on the grid point, the value at xD is calculated by spatial interpolation from nearby grid points. In the right-hand side of Eq. (I6), the value at xD with the subscript D is calculated by third-order Lagrange interpolation.

Equations (I1) and (I2) are converted into

(I7a)v++βvΔt2gh+=Rv,(I7b)Rvv0+2Ω×r-Δt2gh+-βvh++βvh0D-2Ω×r-Δt2gh0-βvh0,

(I8a)h++βhΔt2h¯D+=Rh,(I8b)Rhh0+Δt2-h-hsD++vhs++βhh¯D+-βhh¯D0D+Δt2-h-hsD0+vhs0+βhh¯D0.

In Eqs. (I7b) and (I8b), the values at xD with the subscript D are calculated by fifth-order and third-order Lagrange interpolations, respectively, since high-order interpolation of wind vector components increases the accuracy of the model's results in our experiments. From Eq. (I7), we obtain

(I9)D++βvΔt2g2h+=RD,(I10)ζ+=Rζ,

where

(I11) ζ k × v = 1 a 1 cos ϕ v λ - 1 cos ϕ u cos ϕ ϕ ,

is vorticity, kr/r is the vertical unit vector, RDRv, and Rζk×Rv.

We calculate h+ and v+ using the spectral transform method and the Galerkin method with the new DFS method as follows (see Sect. 3.10 for the spectral transform with the new DFS method).

  1. The scalar variable Rh is transformed from grid space to spectral space using Eqs. (9), (23), (30), and (31). The components of the vector variable Rv=Ru,Rv in grid space are transformed to Rχ and Rψ in spectral space using Eqs. (56), (57), (63), and (64), where Rχ and Rψ are the velocity potential and the stream function of Rv, respectively.

  2. RD and Rζ are calculated by

    (I12)RD=2Rχ,(I13)Rζ=2Rψ,

    using Eqs. (75) and (77). ζ+ is obtained from Rζ using Eq. (I10).

  3. Equations (I8a) and (I12) are substituted into Eq. (I9) and we obtain

    (I14) D + - Δ t 2 2 β v β h g h ¯ 2 D + = 2 R χ - Δ t 2 β v g R h .

    D+ is calculated by solving Eq. (I14) using Eqs. (83) and (85).

  4. h+ is calculated from D+ and Rh using Eq. (I8).

  5. χ+ and ψ+ are calculated from D+ and ζ+ by solving the following Poisson equations

    (I15)2χ+=D+,(I16)2ψ+=ζ+,

    using Eqs. (75) and (78).

  6. v+=u+,v+ is calculated from χ+ and ψ+ using Eq. (53) for un,mc and the similar equations for un,ms, vn,mc, and vn,ms.

  7. u+, v+, h+, D+, and h+ in spectral space are transformed to grid space. h+ and D+ are transformed meridionally using Eqs. (14) and (13). u+ and v+ are transformed meridionally using Eq. (52). h+=hλ+,hθ+ is transformed meridionally using Eqs. (18)–(21). hλ+ can also be calculated from hm+c,Nθj and hm+s,Nθj at the latitudinal grid points using Eq. (16) and additionally using Eq. (22) at the poles when using Grid [1], which is more efficient than using Eqs. (18) and (19) because the meridional inverse discrete cosine and sine transforms of hλ+ become unnecessary.

Code availability

The source code of the DFS and SH shallow-water models is available in the Supplement and is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0) license. These models utilize the Netlib BIHAR library and the ISPACK library. The Netlib BIHAR library is available at https://www.netlib.org/bihar/ (Bihar, 2022) and is also included in the Supplement. The ISPACK library is available at https://www.gfd-dennou.org/arch/ispack/ (Ishioka, 2022).

Data availability

The results of model experiments are available at https://climate.mri-jma.go.jp/pub/archives/Yoshimura_DFS_SW_Testcase/ (Yoshimura, 2022).

Supplement

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-15-2561-2022-supplement.

Competing interests

The contact author has declared that there are no competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We are grateful to Keiichi Ishioka (Kyoto University), Tadashi Tsuyuki (MRI), Daisuke Hotta (MRI), Masashi Ujiie (JMA), and other members of the MRI and JMA model development teams for their useful comments. We are also grateful to the referees and the editor for their comments that helped to improve the quality of this paper.

Financial support

This research has been supported by the Integrated Research Program for Advancing Climate Models (TOUGOU) (grant no. JPMXD0717935561) from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan.

Review statement

This paper was edited by Ignacio Pisso and reviewed by three anonymous referees.

References

Akahori, K., Ishiguro, T., Hattori, K., Suda, R., and Sugihara, M.: A global model based on a double Fourier series, Extended abstract of the 3rd International Workshop on Next Generation Climate Models for Advanced High Performance Computing Facilities, Japan, March 2001. 

Asselin, R. A.: Frequency filter for time integrations, Mon. Weather Rev., 100, 487–490, https://doi.org/10.1175/1520-0493(1972)100<0487:FFFTI>2.3.CO;2, 1972. 

Bihar: Netlib BIHAR library, Netlib Repository [code], https://www.netlib.org/bihar/, last access: 20 March 2022. 

Boer, G. J. and Steinberg, L.: Fourier series on spheres, Atmosphere, 13, 180–191, https://doi.org/10.1080/00046973.1975.9648396, 1975. 

Boyd, J. P.: The choice of spectral functions on a sphere for boundary and eigenvalue problems: A comparison of Chebyshev, Fourier and associated Legendre expansions, Mon. Weather Rev., 106, 1184–1191, https://doi.org/10.1175/1520-0493(1978)106<1184:TCOSFO>2.0.CO;2, 1978. 

Browning, G. L., Hack, J. J., and Swarztrauber, P. N.: A comparison of three numerical methods for solving differential equations on the sphere, Mon. Weather Rev., 117, 1058–1075, https://doi.org/10.1175/1520-0493(1989)117<1058:ACOTNM>2.0.CO;2, 1989. 

Cheong, H.-B.: Double Fourier series on a sphere: applications to elliptic and vorticity equations, J. Comput. Phys., 157, 327–349, https://doi.org/10.1006/jcph.1999.6385, 2000a. 

Cheong, H.-B.: Application of double Fourier series to shallow water equations on a sphere, J. Comput. Phys., 165, 261–287, https://doi.org/10.1006/jcph.2000.6615, 2000b. 

Cheong, H.-B.: A dynamical core with double Fourier series: Comparison with the spherical harmonics method, Mon. Weather Rev., 134, 1299–1315, https://doi.org/10.1175/MWR3121.1, 2006. 

Cheong, H.-B., Kwon, I.-H., and Goo, T.-Y.: Further study on the high-order double-Fourier-series spectral filtering on a sphere, J. Comput. Phys., 193, 180–197, https://doi.org/10.1016/j.jcp.2003.07.029, 2004. 

Cheong, H.-B., Kwon, I.-H., Goo, T.-Y., and M.-J. Lee: High-order harmonic spectral filter with the double Fourier series on a sphere, J. Comput. Phys., 177, 313–335, https://doi.org/10.1006/jcph.2002.6997, 2002. 

Cooley, J. W. and Tukey, J. W.: An algorithm for the machine calculation of complex Fourier series, Math. Comput., 19, 297–301, https://doi.org/10.2307/2003354, 1965. 

Dueben, P. D., Wedi, N., Saarinen, S., and Zeman, C.: Global simulations of the atmosphere at 1.45 km grid-spacing with the Integrated Forecasting System, J. Meteorol. Soc. Jpn., 98, 551–572, https://doi.org/10.2151/jmsj.2020-016, 2020. 

Galewsky, J., Scott, R. K., and Polvani, L. M.: An initial-value problem for testing numerical models of the global shallow-water equations, Tellus A, 56, 429–440, https://doi.org/10.3402/tellusa.v56i5.14436, 2004. 

Gander, W. and Golub, G. H.: Cyclic reduction – History and applications, Proceedings of the Workshop on Scientific Computing, Hong Kong, 10–12 March, 1997. 

Gospodinov, I. G., Spiridonov, V. G., and Geleyn, J.-F.: Second-order accuracy of two-time-level semi-Lagrangian schemes, Q. J. Roy. Meteor. Soc., 127, 1017–1033, https://doi.org/10.1002/qj.49712757317, 2001. 

Hortal, M.: The development and testing of a new two-time-level semi-Lagrangian scheme (SETTLS) in the ECMWF forecast model, Q. J. Roy. Meteor. Soc., 128, 1671–1687, https://doi.org/10.1002/qj.200212858314, 2002. 

Hortal, M. and Simmons, A. J.: Use of reduced Gaussian grids in spectral models, Mon. Weather Rev., 119, 1057–1074, https://doi.org/10.1175/1520-0493(1991)119<1057:UORGGI>2.0.CO;2, 1991. 

Hotta, D. and Ujiie, M.: A nestable, multigrid-friendly grid on a sphere for global spectral models based on Clenshaw–Curtis quadrature, Q. J. Roy. Meteor. Soc., 144, 1382–1397, https://doi.org/10.1002/qj.3282, 2018. 

Ishioka, K.: A New Recurrence Formula for Efficient Computation, J. Meteorol. Soc. Jpn., 96, 241–249, https://doi.org/10.2151/jmsj.2018-019, 2018. 

Ishioka, K.: ISPACK library, GFD-DENNOU Club [code], https://www.gfd-dennou.org/arch/ispack/, last access: 20 March 2022. 

Japan Meteorological Agency (JMA): Outline of the operational numerical weather prediction at the Japan Meteorological Agency, Appendix to WMO technical progress report on the global data-processing and forecasting system and numerical weather prediction, 229 pp., http://www.jma.go.jp/jma/jma-eng/jma-center/nwp/outline2019-nwp/index.htm (last access: 20 March 2022), 2019. 

Juang, H.-M. H.: A reduced spectral transform for the NCEP seasonal forecast global spectral atmospheric model, Mon. Weather Rev., 132, 1019–1035, https://doi.org/10.1175/1520-0493(2004)132<1019:ARSTFT>2.0.CO;2, 2004. 

Koo, M.-S. and Hong, S.-Y.: Double Fourier series dynamical core with hybrid sigma-pressure vertical coordinate, Tellus A, 65, 19851, https://doi.org/10.3402/tellusa.v65i0.19851, 2013. 

Kwon, I.-H., Cheong, H.-B., Joh, M., Chung, I.-U., Cho, C.-H., and Lee, W.-J.: Application of double-Fourier-series spectral method to a large size problem: Two dimensional simulations of the shear instability on the sphere, J. Meteorol. Soc. Jpn., 82, 1301–1314, https://doi.org/10.2151/jmsj.2004.1301, 2004. 

Lambert, S. J.: A global available potential energy-kinetic energy budget in terms of the two-dimensional wavenumber for the FGGE year, Atmos. Ocean, 22, 265–282, https://doi.org/10.1080/07055900.1984.9649199, 1984. 

Layton, A. T. and Spotz, W. F.: A semi-Lagrangian double Fourier method for the shallow water equations on the sphere, J. Comput. Phys., 189, 180–196, https://doi.org/10.1016/S0021-9991(03)00207-9, 2003. 

Malardel, S., Wedi, N., Deconinck, W., Diamantakis, M., Kühnlein, C., Mozdzynski, G., Hamrud, M., and Smolarkiewicz, P.: A new grid for the IFS, ECMWF Newsletter, 146, 23–28, https://doi.org/10.21957/zwdu9u5i, 2016. 

Merilees, P. E.: An alternative scheme for the simulation of a series of spherical harmonics, J. Appl. Meteorol. Clim., 12, 224–227, https://doi.org/10.1175/1520-0450(1973)012<0224:AASFTS>2.0.CO;2, 1973a. 

Merilees, P. E.: The pseudospectral approximation applied to the shallow water equations on a sphere, Atmosphere, 11, 13–20, https://doi.org/10.1080/00046973.1973.9648342, 1973b. 

Merilees, P. E.: Numerical experiments with the pseudospectral method in spherical coordinates, Atmosphere, 12, 77–96, https://doi.org/10.1080/00046973.1974.9648374, 1974. 

Miyamoto, K.: Introduction of the Reduced Gaussian Grid into the Operational Global NWP model at JMA, CAS/JSC WGNE Research Activities in Atmospheric and Ocean Modelling, 36, 6.9–6.10, 2006. 

Nakano, M., Wada, A., Sawada, M., Yoshimura, H., Onishi, R., Kawahara, S., Sasaki, W., Nasuno, T., Yamaguchi, M., Iriguchi, T., Sugi, M., and Takeuchi, Y.: Global 7 km mesh nonhydrostatic Model Intercomparison Project for improving TYphoon forecast (TYMIP-G7): experimental design and preliminary results, Geosci. Model Dev., 10, 1363–1381, https://doi.org/10.5194/gmd-10-1363-2017, 2017. 

Orszag, S. A.: On the elimination of aliasing in finite-difference schemes by filtering high-wavenumber components, J. Atmos. Sci., 28, 1074–1074, https://doi.org/10.1175/1520-0469(1971)028<1074:OTEOAI>2.0.CO;2, 1971. 

Orszag, S. A.: Fourier series on spheres, Mon. Weather Rev., 102, 56–75, https://doi.org/10.1175/1520-0493(1974)102<0056:FSOS>2.0.CO;2, 1974. 

Park, H., Hong, S.-Y., Cheong, H.-B., and Koo, M.-S.: A Double Fourier Series (DFS) Dynamical Core in a Global Atmospheric Model with Full Physics. Mon. Weather Rev., 141, 3052–3061, https://doi.org/10.1175/MWR-D-12-00270.1, 2013. 

Ritchie, H.: Application of the semi-Lagrangian method to a spectral model of the shallow water equations, Mon. Weather Rev., 116, 1587–1598, https://doi.org/10.1175/1520-0493(1988)116<1587:AOTSLM>2.0.CO;2, 1988. 

Ritchie, H. and Tanguay, M.: A comparison of spatially averaged Eulerian and semi-Lagrangian treatments of mountains, Mon. Weather Rev., 124, 167–181, https://doi.org/10.1175/1520-0493(1996)124<0167:ACOSAE>2.0.CO;2, 1996. 

Ritchie, H., Temperton, C., and Simmons, A.: Implementation of the semi-Lagrangian method in a high-resolution version of the ECMWF forecast model, Mon. Weather Rev., 123, 489–514, https://doi.org/10.1175/1520-0493(1995)123<0489:IOTSLM>2.0.CO;2, 1995. 

Robert, A. J.: The integration of a low order spectal form of the primitive meteorological equations, J. Meteorol. Soc. Jpn., 44, 237–245, https://doi.org/10.2151/jmsj1965.44.5_237, 1966. 

Schaeffer, N.: Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations, Geochem. Geophys. Geosy., 14, 751–758, https://doi.org/10.1002/ggge.20071, 2013. 

Sneeuw, N. and Bun, R.: Global spherical harmonic computation by two-dimensional Fourier methods, J. Geodesy, 70, 224–232, https://doi.org/10.1007/BF00873703, 1996. 

Spotz, W. F., Taylor, M. A., and Swarztrauber, P. N.: Fast shallow-water equations solvers in latitude-longitude coordinates, J. Comput. Phys., 145, 432–444, https://doi.org/10.1006/jcph.1998.6026, 1998. 

Suda, R.: Fast spherical harmonic transform routine FLTSS applied to the shallow water test set, Mon. Weather Rev., 133, 634–648, https://doi.org/10.1175/MWR-2871.1, 2005. 

Swarztrauber, P. N.: Vectorizing the FFTs, edited by: Rodrigue, G., Parallel Comput., Academic Press, 51–83, https://doi.org/10.1016/B978-0-12-592101-5.50007-5, 1982. 

Swarztrauber, P. N.: The vector harmonic transform method for solving partial differential equations in spherical geometry, Mon. Weather Rev., 121, 3415–3437, https://doi.org/10.1175/1520-0493(1993)121<3415:TVHTMF>2.0.CO;2, 1993. 

Swarztrauber, P. N.: Shallow water flow on the sphere, Mon. Weather Rev., 132, 3010–3018, https://doi.org/10.1175/MWR2829.1, 2004. 

Swarztrauber, P. N. and Spotz, W. F.: Generalized discrete spherical harmonic transforms, J. Comput. Phys., 159, 213–230, https://doi.org/10.1006/jcph.2000.6431, 2000. 

Temperton, C.: On scalar and vector transform methods for global spectral models, Mon. Weather Rev., 119, 1303–1307, https://doi.org/10.1175/1520-0493-119-5-1303.1, 1991. 

Temperton, C.: Treatment of the Coriolis terms in semi-Lagrangian spectral models, Atmos. Ocean, 35, 293–302, https://doi.org/10.1080/07055900.1997.9687353, 1997. 

Temperton, C., Hortal, M., and Simmons, A.: A two-time-level semi-Lagrangian global spectral model, Q. J. Roy. Meteor. Soc., 127, 111–127, https://doi.org/10.1002/qj.49712757107, 2001. 

Tygert, M.: Fast algorithms for spherical harmonic expansions, II, J. Comput. Phys., 227, 4260–4279, https://doi.org/10.1016/j.jcp.2007.12.019, 2008. 

Wedi, N. P., Hamrud, M., and Mozdzynski, G.: A fast spherical harmonics transform for global NWP and climate models, Mon. Weather Rev., 141, 3450–3461, https://doi.org/10.1175/MWR-D-13-00016.1, 2013. 

Wedi, N. P.: Increasing horizontal resolution in numerical weather prediction and climate simulations: illusion or panacea?, Philos. T. R. Soc. A, 372, 20130289, https://doi.org/10.1098/rsta.2013.0289, 2014. 

Williamson, D. L., Drake, J. B., Hack, J. J., Jacob, R., and Swarztrauber, P. N.: A standard test set for numerical approximations to the shallow water equations in spherical geometry, J. Comput. Phys., 102, 211–224, https://doi.org/10.1016/S0021-9991(05)80016-6, 1992. 

Yee, S. Y. K.: Solution of Poisson's equation on a sphere by truncated double Fourier series, Mon. Weather Rev., 109, 501–505, https://doi.org/10.1175/1520-0493(1981)109<0501:SOPEOA>2.0.CO;2, 1981.  

Yoshimura, H.: Development of a nonhydrostatic global spectral atmospheric model using double Fourier series, CAS/JSC WGNE Research Activities in Atmospheric and Ocean Modelling, 42, 3.05–3.06, 2012. 

Yoshimura, H.: DFS shallow-water test cases, Meteorological Research Institute [data set], https://climate.mri-jma.go.jp/pub/archives/Yoshimura_DFS_SW_Testcase/, last access: 20 March 2022. 

Yoshimura, H. and Matsumura, T.: A two-time-level vertically conservative semi-Lagrangian semi-implicit double Fourier series AGCM, CAS/JSC WGNE Research Activities in Atmospheric and Ocean Modelling, 35, 3.25–3.26, 2005. 

Yukimoto, S., Yoshimura, H., Hosaka, M., Sakami, T., Tsujino, H., Hirabara, M., Tanaka, T. Y., Deushi, M., Obata, A., Nakano, H., Adachi, Y., Shindo, E., Yabu, S., Ose, T., and Kitoh, A.: Meteorological Research Institute-Earth System Model Version 1 (MRI-ESM1) – Model Description –, Technical Reports of the Meteorological Research Institute, No. 64, https://doi.org/10.11483/mritechrepo.64, 2011. 

Yukimoto, S., Kawai, H., Koshiro, T., Oshima, N., Yoshida, K., Urakawa, S., Tsujino, H., Deushi, M., Tanaka, T., Hosaka, M., Yabu, S., Yoshimura, H., Shindo, E., Mizuta, R., Obata, A., Adachi, Y., and Ishii, M.: The Meteorological Research Institute Earth System Model Version 2.0, MRI-ESM2.0: Description and Basic Evaluation of the Physical Component, J. Meteorol. Soc. Jpn., 97, 931–965, https://doi.org/10.2151/jmsj.2019-051, 2019. 

Download
Short summary
This paper proposes a new double Fourier series (DFS) method on a sphere that improves the numerical stability of a model compared with conventional DFS methods. The shallow-water model and the advection model using the new DFS method give stable results without the appearance of high-wavenumber noise near the poles. The model using the new DFS method is faster than the model using spherical harmonics (especially at high resolutions) and gives almost the same results.