the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Improved double Fourier series on a sphere and its application to a semi-implicit semi-Lagrangian shallow-water model
Hiromasa Yoshimura
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.
- Article
(5491 KB) - Full-text XML
-
Supplement
(495 KB) - BibTeX
- EndNote
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
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
where N is the meridional truncation wavenumber. The coefficients fn,m for odd m are calculated by the forward Fourier cosine transform of . Orszag (1974) imposed the following conditions at the poles:
which can be expressed in terms of the expansion coefficients fn,m as
Satisfying the above conditions ensures that the scalar variable F(λ,θ) and its gradient ∇F are continuous at the poles. In Orszag (1974), only 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
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
The meridional basis functions sin θsin nθ for even m (≠0) are different from Eq. (5). The coefficients are calculated by the forward Fourier sine transform of . 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 or 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.
-
The semi-Lagrangian scheme is used, which avoids the numerical instability due to the nonlinear advection term.
-
The meridional truncation with and 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.
-
U=usin θ and V=vsin θ instead of and 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.
-
High-wavenumber noise appears near the poles.
-
The meridional truncation wavenumber N′ needs to be equal to J for even m (≠0) because (e.g., 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.
-
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.
-
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.
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:
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 . 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).
In the new DFS method, the wind vector components u and v (instead of and 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.
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
where and are calculated from T(λ,θ) by the forward Fourier transform as
The variables and are meridionally expanded as
where
Here, the superscript c(s) means c or s, and thus, for example, means or . 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 are the linear combination of 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:
where
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 for m= 0 or 1. The reason for this is shown in Appendix C.
When calculating the values of in grid space from in spectral space, the coefficients are calculated from using Eq. (14), and inverse discrete cosine and sine transforms (see Appendix B) are performed using Eq. (13). The calculation of in spectral space from in grid space is described in Sect. 3.3 below.
The truncated variable is defined as
From Eq. (10), the values of at the poles are finite for m=0, and the values of at the poles are zero for m≠0. Therefore, is continuous at the poles.
3.2 Gradient of a scalar variable
The gradient is obtained as follows:
where a is the radius of the Earth and ϕ is the latitude. From Eqs. (16b), (16c), and (A2b), we obtain
where
The equations for are the same as Eq. (19), except that and are replaced with and , respectively. From Eqs. (17b), (13), and (14), we obtain
where
From Eqs. (18)–(21), it can be seen that , , , and at the poles are finite for m=1 and zero for m≠1, and moreover the following relations are satisfied for m=1:
Thus, it is guaranteed that is continuous at the poles.
3.3 New method to calculate expansion coefficients for a scalar variable
One way to calculate the coefficients from in Eq. (10) is to perform a forward cosine transform of for m=1, a sine transform of for even m (≥2), and a sine transform of for odd m (≥3). However, this approach with the meridional wavenumber truncation N<J leads to large high-wavenumber oscillations (see Sect. 4). Dividing by sin 2θ reduces the numerical stability of the model more significantly than dividing 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 by sin θ or sin 2θ before the forward cosine or sine transforms. The coefficients in Eq. (10) are calculated as follows. First, in Eq. (10) are expanded like Eq. (5) as
where the expansion coefficients are calculated by the forward discrete cosine transform for even m and the forward discrete sine transform for odd m from the values of at the grid points (see Appendix B).
Next, and are calculated using the least-squares method to minimize the following error E (the squared L2 norm of the residual):
where the residual R(λ,θ) is
From Eqs. (24), (25), and (15) and the equations and used in the least-squares method to minimize E, we obtain
From Eq. (10), we derive
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
From Eqs. (28) and (D4) in Appendix D, we obtain
By substituting Eqs. (10) and (23) into Eq. (29), the following equations for and are derived, as shown in Appendix E.
For m=0,
For m=1,
with the exception of the following underlined values:
For even m (≥2),
with the exception of the following underlined values:
For odd m(≥3),
with the exception of the following underlined values:
From Eq. (30a), for m=0 is obtained. From Eqs. (30d), the following linear simultaneous equations for m≥3 are derived:
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 are calculated from . A penta-diagonal matrix C can be lower–upper (LU) decomposed as
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 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, and in Eq. (8) are expanded with the associated Legendre functions Pn,m(θ) as
where m≥0. The functions Pn,m(θ) satisfy the following orthogonality relations for each m:
By the modified Robert expansion (Merilees, 1973a; Orszag, 1974), the associated Legendre functions Pn,m(θ) are expressed as
Conversely, the functions can be expressed as the linear combination of . Substituting Eq. (35) into Eq. (33) gives
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 for m=2 and for m=3 in Eq. (36) are the linear combinations of and 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 in the SH model is about 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 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 in Eq. (33) are calculated from by the forward Legendre transform as
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 is replaced with (e.g., Sneeuw and Bun, 1996), although the values of calculated from are different from those calculated from . In the new DFS method, the values of calculated from in Eq. (29) are the same as those calculated from 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):
where the residual RSH(λ,θ) is
From Eqs. (38), (39), and (33) and the equations and used in the least-squares method to minimize ESH, we derive
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 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 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 , which satisfy
This is similar to Eq. (34), but the latitudinal weight is constant. in Eq. (8) are expanded with as
By using the least squares method or the Galerkin method with Eq. (42), we obtain the same equations as Eqs. (24)–(29), except that and Sn,m(θ) are replaced with and , respectively. From Eq. (29) with and Sn,m(θ) replaced by and and Eqs. (41) and (42), we derive
Thus, and in Eqs. (43) and (42) are calculated uniquely. This unique solution is the same as in Eq. (29) obtained by the least-squares method with the non-orthogonal basis functions Sn,m(θ) because is the linear combination of 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
where is the zonal wind and is the meridional wind. The scalar variables χ and ψ are expanded like Eqs. (8) and (10) as
The truncated variables and are defined as
From Eqs. (44)–(47), the equations for the wind vector components and are derived as
The vector in Eqs. (48) and (49) can also be represented as
where we define the new DFS vector basis functions , , , and as
From Eqs. (48), (49), and (16)–(21), we obtain
where
The equations for are the same as Eqs. (53b)–(53d), except that , , and are replaced with , , and , respectively. The equations for are the same as Eqs. (53a)–(53d), except that , , and are replaced with , , and , respectively. The equations for are the same as Eqs. (53b)–(53d), except that , , and are replaced with , , and , respectively.
From Eqs. (52) and (53), it can be seen that , , , and at the poles are finite for m=1 and zero for m≠1. Moreover, the following relations are satisfied for m=1:
Thus, it is guaranteed that the wind vector 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 , , , and in Eqs. (48)–(50) using the least-squares method to minimize the error of and from u(λ,θ) and v(λ,θ) due to the meridional wavenumber truncation. First, the wind vector components u and v are expanded zonally as
where and are calculated from u(λ,θ) and v(λ,θ) by the forward Fourier transform as
The variables and are meridionally expanded as
where and are calculated from and by the forward discrete cosine or sine transform (see Appendix B).
Next, , , , and are calculated to minimize the following error F (the squared L2 norm of the residual vector):
where the residual vector is defined as
From Eqs. (58) and (59) and the equations , , , and used in the least-squares method, we obtain
From Eq. (50), we derive
Equations (60) and (61) show that the residual vector 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
By substituting Eqs. (52) and (57) into Eq. (62a), the following equations for and are derived as shown in Appendix H.
For m=0,
The coefficient 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,
with the exception of the following underlined values:
For even m≥2,
with no exception.
For odd m≥3,
with the exception of the following underlined values:
Similarly, from Eq. (62b) we derive the same equations as Eqs. (63b)–(63d), except that χc, ψs, , and are replaced with χs, −ψc, , and , respectively. From Eq. (62c), we derive the same equations as Eqs. (63a)–(63d), except that χc, ψs, , and are replaced with −ψc, χs, , and , respectively. From Eq. (62d), we derive the same equations as Eqs. (63b)–(63d), except that χc, ψs, , and are replaced with ψs, χc, , and , respectively.
Equation (63a) is easily solved. From Eq. (63d) and from the same equations as Eq. (63d), except that χc, ψs, , and are replaced with ψs, χc, , and , respectively, we derive the following linear simultaneous equations for m≥3:
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 and in Eq. (64) can be calculated efficiently. From the same equations as Eqs. (63b)–(63d), except that χc, ψs, , and are replaced with χs, −ψc, , and , respectively, and the same equations as Eqs. (63b)–(63d), except that χc, ψs, , and are replaced with −ψc, χs, , and , respectively, the simultaneous equations similar to Eq. (64) are also derived. Thus, the expansion coefficients , , , and are calculated from , , , and using Eqs. (63a)–(63d) and the similar equations. The expansion coefficients , , , and are calculated from , , , and using Eq. (53) for and the similar equations for , , and .
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
where ∇2 is the Laplacian operator, the variables f and g are expanded zonally like Eq. (8) as
The variables , , , and are expanded meridionally like Eq. (10) as
We define the truncated variables fN,M(θ) and gN,M(θ) as
From Eqs. (65) and (68), we obtain
Here we use the Galerkin method to calculate the Laplacian operator and the Poisson equation and obtain
where the residual is as follows:
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:
When calculating g by applying the Laplacian operator to a given f, and can also be calculated from and using the least-squares method. The equations and give the equivalent equations to Eq. (70). When calculating f from a given g in the Poisson equation, and can also be calculated from and using the least-squares method. However, the equations derived from and 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
For m=0, we calculate by using
instead of Eq. (73) following Yee (1981) and Cheong (2000a) for ease of calculation. For , the exact solutions of can be obtained from Eq. (74) because the new DFS meridional basis functions for are the linear combination of the associated Legendre functions for 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
except for the following underlined values:
For m=1, by substituting Eq. (67) into Eq. (73) and using Eqs. (A2d), (A4a), and (A5b), we obtain
except for the following underlined values:
For even m≥2, by substituting Eq. (67) into Eq. (73) and using Eqs. (A2c), (A4b), and (A5d), we obtain
except for the following underlined values:
For odd m≥3, by substituting Eq. (67) into Eq. (73) and using Eqs. (A2c), (A2e), (A4b), and (A5d), we obtain
except for the following underlined values:
From Eq. (75), we obtain the following two linear simultaneous equations with tri-diagonal or penta-diagonal matrices:
where and are the vectors whose components are (n is even) and (n is odd), respectively, and and are the vectors whose components are (n is even) and (n is odd), respectively; An_even,m, Bn_even,m, An_odd,m and Bn_odd,m are tri-diagonal or penta-diagonal matrices. and are calculated by
which can be solved efficiently as in Eq. (31). We have verified that all the eigenvalues of the matrices and 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
except when m=0 and n is even. For m=0, disappears in Eq. (75a). The coefficients (even n≥2) are calculated from (even n≥2) by using Eq. (75a). The value is calculated from (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
where f is calculated from given g. From Eq. (76), the Poisson equation in Eq. (65) is represented as
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
From Eq. (81), f is calculated from g by
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
is represented as
From Eq. (84), f is calculated from g by
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
where f is calculated from g. Equation (86) can be converted into
where . The calculation of Eq. (86) is accomplished by successive calculations of the following Helmholtz equations:
which are represented as
From Eqs. (89), we obtain the equation to calculate f from g as
Here, and 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, and , 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.
-
Define DFS expansion for a scalar variable with zonal expansion in Eq. (8) and meridional expansion in Eq. (10).
-
For the inverse transform from spectral space to grid point space,
- a.
calculate the coefficients from by using Eqs. (14),
- b.
calculate in Eq. (13) from by inverse cosine and sine Fourier transforms in Appendix B,
- c.
calculate the grid point values in Eq. (15) from by inverse Fourier transform.
- a.
-
For the forward transform from grid point space to spectral space,
- a.
calculate in Eq. (9) from the grid point values T(λ,θ) by forward Fourier transform,
- b.
calculate the coefficients in Eq. (23) from by forward cosine and sine transforms in Appendix B,
- c.
calculate the coefficients from by using Eqs. (30) and (31), where the coefficients are calculated so that Eq. (29) derived from the Galerkin method (or the least-squares method) is satisfied.
- a.
The essential summary for a vector variable is as follows.
-
Represent DFS expansion for a vector variable by Eq. (50).
-
For the inverse transform from spectral space to grid point space,
- a.
calculate the coefficients and from and by using Eq. (53) and the similar equations,
- b.
calculate and in Eq. (52) from and by inverse cosine and sine transforms in Appendix B,
- c.
calculate the grid point values and in Eqs. (48a) and (49a) from and by inverse Fourier transform.
- a.
-
For the forward transform from grid point space to spectral space,
- a.
calculate and in Eq. (56) from the grid point values u(λ,θ) and v(λ,θ) by forward Fourier transform,
- b.
calculate the coefficients and in Eq. (57) from the grid point values and by forward cosine and sine transforms in Appendix B,
- c.
calculate the coefficients and from and by using Eqs. (63) and (64), where and are calculated so that Eq. (62) derived from the Galerkin method (or the least-squares method) is satisfied.
- a.
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 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).
In the DFS method of Orszag, a very large error occurs, especially for odd |m| (≥3) (Fig. 2c), when 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 for m≠0. Because 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 for even m (≥2) to improve stability. However, there is a problem with the latitudinal derivative in the old DFS method even when 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.
We examine the accuracy of the old and new DFS methods for the Laplacian operator in Eq. (65) and the Helmholtz equation
Here, we give the function f as
where , ϕ is latitude, λ is longitude, a is the radius of the Earth, and r is the distance between (λ,ϕ) and the center . The function f is similar to the cosine bell in Williamson test case 1 (Williamson et al., 1992), but 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, and , where θ is colatitude. At this time, g is calculated as follows:
Equation (94) is satisfied at any position of the center. The function h in Eq. (91) is calculated by
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.
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
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
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
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
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 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 . 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 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 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 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.
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
Figure 5 shows the time series of forecast errors of the height for a 5 d integration in Williamson test case 2 with 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.
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).
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 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.
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 . In our semi-implicit semi-Lagrangian models, we usually use N satisfying , which is called linear truncation. However, here N is determined to satisfy 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 . The results of the new DFS models and the SH model with in Table 8 are very similar to those with in Table 7 when J0 is the same.
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 . 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.
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.
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.
We have developed the new DFS method to improve the numerical stability of the DFS model, which has the following two improvements.
-
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.
-
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, 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 or cubic truncation are usually used in the Eulerian model and are also starting to be used in the semi-Lagrangian model instead of the linear truncation 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.
Here we list the trigonometric identities used in transforming the expressions in this paper. The following identities are satisfied:
From Eq. (A1), the following identities are derived:
From Eq. (A1), the following orthogonal relations in longitude are derived:
Similarly, from Eq. (A1), the following orthogonal relations in latitude are derived:
By using Eqs. (A1) and (A2), the following relations are derived:
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:
When using Grid [0], inverse and forward discrete sine transforms are performed as follows:
When using Grid [1], inverse and forward discrete cosine transforms are performed as follows:
When using Grid [1], inverse and forward discrete sine transforms are performed as follows:
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:
where n is from 0 to because the number of the meridional grid points is 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 as follows:
We can expand as follows:
The expansion coefficients can be obtained from in the same way as in Eq. (B4b) by forward discrete sine transform:
From Eqs. (B6) and (B7), we obtain
By using Eq. (A2a), we obtain
By substituting Eq. (B10) into Eq. (B9) and comparing the left and right sides of the equation, we obtain
We can calculate from g(θj) using Eq. (B6), calculate from using Eq. (B8), and calculate gn from using Eq. (B11).
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.
in Eq. (23) is calculated by the forward Fourier cosine or sine transform as follows:
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:
From Eqs. (D1) and (D2), the following equations are satisfied:
From Eqs. (D3), (11), (12), and (A2a)–(A2c), we derive
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):
From Eqs. (D5), (11), (12), and (A2a)–(A2c), we derive
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).
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:
where . From Eqs. (29), (E1), (E2), and (A4b), Eq. (30d) is derived.
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:
The values of and 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.
The global mean value of in Eq. (15) can be calculated in spectral space by the following equation (Cheong 2000a):
where Eq. (A2a) is used.
The latitudinal area weight at each colatitude θj is calculated as follows.
-
The latitudinal distribution of for each j is given as
where in Grid [0] and Grid [1] and in Grid [−1] (see Sect. 2).
-
From , the meridional expansion coefficients are calculated by forward discrete cosine transform as described in Appendix B.
-
The value of G calculated from 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.
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 , , and are replaced with , , and , respectively, and the same equations as Eq. (53), except that , , and are replaced with , , and , respectively, Eq. (62a) is converted into the following equation for odd m≥3:
When n≥4, by using Eq. (A4a), Eq. (H1) can be converted into
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).
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:
where
is horizontal divergence; Δt is a time step; the superscripts −, 0, and + mean past time (t−Δt), present time (t), and future time (t+Δt), respectively; the superscript (+) means future time (t+Δt) extrapolated in time, for example, ; the subscript D means the departure point; the absence of the subscript D means the arrival point; 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 is used here because 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 (t+Δt). 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 (t+Δt) 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
However, when Δt is longer than 30 min, using extrapolated in time to calculate xD causes numerical instability in our experiments. To avoid instability when Δt is 1 h, here we use
instead of Eq. (I4), where is a provisional future value obtained by discretizing Eq. (98) in an explicit semi-Lagrangian scheme. From Eq. (I5), we obtain
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
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
where
is vorticity, is the vertical unit vector, , and .
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).
-
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 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.
-
RD and Rζ are calculated by
using Eqs. (75) and (77). ζ+ is obtained from Rζ using Eq. (I10).
-
Equations (I8a) and (I12) are substituted into Eq. (I9) and we obtain
D+ is calculated by solving Eq. (I14) using Eqs. (83) and (85).
-
h+ is calculated from D+ and Rh using Eq. (I8).
-
χ+ and ψ+ are calculated from D+ and ζ+ by solving the following Poisson equations
using Eqs. (75) and (78).
-
is calculated from χ+ and ψ+ using Eq. (53) for and the similar equations for , , and .
-
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). is transformed meridionally using Eqs. (18)–(21). can also be calculated from and 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 become unnecessary.
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).
The results of model experiments are available at https://climate.mri-jma.go.jp/pub/archives/Yoshimura_DFS_SW_Testcase/ (Yoshimura, 2022).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd-15-2561-2022-supplement.
The contact author has declared that there are no competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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.
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.
This paper was edited by Ignacio Pisso and reviewed by three anonymous referees.
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.
- Abstract
- Introduction
- Arrangement of equally spaced latitudinal grid points
- Improved double Fourier series on the sphere
- The error due to meridional wavenumber truncation in DFS expansion methods
- Tests of the DFS methods with the Laplacian operator and the Helmholtz equation
- Evaluation of the DFS methods using shallow-water test cases
- Conclusions and perspectives
- Appendix A: Trigonometric identities
- Appendix B: Discrete Fourier cosine and sine transforms in latitude
- Appendix C: The upper limit of the meridional truncation wavenumber N
- Appendix D: Equations for the derivation of Eqs. (29) and (62)
- Appendix E: Derivation of Eq. (30) from Eq. (29)
- Appendix F: Zonal Fourier filter
- Appendix G: Calculation of global mean and latitudinal area weight
- Appendix H: Derivation of Eq. (63) from Eq. (62)
- Appendix I: Two-time-level semi-implicit semi-Lagrangian scheme for time integration
- Code availability
- Data availability
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References
- Supplement
- Abstract
- Introduction
- Arrangement of equally spaced latitudinal grid points
- Improved double Fourier series on the sphere
- The error due to meridional wavenumber truncation in DFS expansion methods
- Tests of the DFS methods with the Laplacian operator and the Helmholtz equation
- Evaluation of the DFS methods using shallow-water test cases
- Conclusions and perspectives
- Appendix A: Trigonometric identities
- Appendix B: Discrete Fourier cosine and sine transforms in latitude
- Appendix C: The upper limit of the meridional truncation wavenumber N
- Appendix D: Equations for the derivation of Eqs. (29) and (62)
- Appendix E: Derivation of Eq. (30) from Eq. (29)
- Appendix F: Zonal Fourier filter
- Appendix G: Calculation of global mean and latitudinal area weight
- Appendix H: Derivation of Eq. (63) from Eq. (62)
- Appendix I: Two-time-level semi-implicit semi-Lagrangian scheme for time integration
- Code availability
- Data availability
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References
- Supplement