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

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(N³) 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(N²log N) operations. This paper proposes a new DFS method that improves the numerical stability of the model compared with the conventional DFS methods by adopting the following two improvements: a new expansion method that employs the least-squares method (or the Galerkin method) to calculate the expansion coefficients in order to minimize the error caused by wavenumber truncation, and new basis functions that satisfy the continuity of both scalar and vector variables at the poles. Partial differential equations such as the Poisson equation and the Helmholtz equation are solved by using the Galerkin method. In the semi-implicit semi-Lagrangian shallow-water model using the new DFS method, the Williamson test cases and the Galewsky test case give stable results without the appearance of high-wavenumber noise near the poles, even without horizontal diffusion and without a zonal Fourier filter. In the Eulerian advection model using the new DFS method, the Williamson test case 1, which simulates a cosine bell advection, also gives stable results without horizontal diffusion but with a zonal Fourier filter. The shallow-water model using the new DFS method is faster than that using SH, especially at high resolutions and gives almost the same results, except that very small oscillations near the truncation wavenumber in the kinetic energy spectrum appear only in the shallow-water model using SH.

1 Introduction

Global spectral atmospheric models using the spectral transform method with spherical harmonics (SH) as basis functions are widely used. They are used in the Japan Meteorological Agency (JMA, 2019) and the Meteorological Research Institute (MRI; Yukimoto et al., 2011, 2019) for a range of applications, including operational weather prediction, operational seasonal prediction, and global warming projection. The spectral model has the advantage that the horizontal derivatives are accurate, and the semi-implicit scheme, which improves numerical stability, can be easily applied because the Helmholtz equation and the Poisson equation are easily solved in spectral space. The application of the semi-implicit semi-Lagrangian scheme allows for time steps longer than the Courant–Friedrichs–Lewy (CFL) condition, which makes the model computationally efficient. In the spectral model using SH, the Legendre transform used in the latitudinal direction significantly increases the computational cost at high resolutions since the Legendre transform usually requires O(N³) operations and O(N³) 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(N²(log N)³) 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(N³) operations but only O(N²) 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(N²log N) operations and O(N) or O(N²) 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

$$\begin{array}{}\text{(1)}& F\left(\mathit{\lambda },\mathit{\theta }\right)\cong \sum _{m=-M}^M{F}_m\left(\mathit{\theta }\right)e^{im\mathit{\lambda }},\end{array}$$

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 F_m(θ) meridionally as

$$\begin{array}{}\text{(2a)}& {\displaystyle }{F}_m\left(\mathit{\theta }\right)\cong \left\{\begin{array}{ll}{f}_m\left(\mathit{\theta }\right)& \text{for even}m,\\ \sin \mathit{\theta }{f}_m\left(\mathit{\theta }\right)& \text{for odd}m,\end{array}\right.\text{(2b)}& {\displaystyle }{f}_m\left(\mathit{\theta }\right)\equiv \sum _{n=\mathrm{0}}^N{f}_{n,m}\cos n\mathit{\theta },\end{array}$$

where N is the meridional truncation wavenumber. The coefficients f_n,m for odd m are calculated by the forward Fourier cosine transform of $F_m\left(\mathit{\theta }\right)/\sin \mathit{\theta }$. Orszag (1974) imposed the following conditions at the poles:

$$\begin{array}{}\text{(3)}& f_m\left(\mathrm{0}\right)=\mathrm{0}\text{ and }{f}_m\left(\mathit{\pi }\right)=\mathrm{0}\text{ for }\left|m\right|\ge \mathrm{2},\end{array}$$

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

$$\begin{array}{}\text{(4)}& \sum _{\begin{array}{l}n=\mathrm{0}\\ n\text{ is even}\end{array}}^N{f}_{n,m}=\mathrm{0}\text{and}\sum _{\begin{array}{l}n=\mathrm{1}\\ n\text{ is\hspace{0.17em} odd}\end{array}}^N{f}_{n,m}=\mathrm{0}\text{ for }\left|m\right|\ge \mathrm{2}.\end{array}$$

Satisfying the above conditions ensures that the scalar variable F(λ,θ) and its gradient ∇F are continuous at the poles. In Orszag (1974), only $f_{N-\mathrm{1},m}$ and f_N,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 F_m(θ) as

$$\begin{array}{}\text{(5)}& F_m\left(\mathit{\theta }\right)=\left\{\begin{array}{ll}\sum _{n=\mathrm{0}}^N{F}_{n,m}\cos n\mathit{\theta }& \text{ for even}m,\\ \sum _{n=\mathrm{1}}^N{F}_{n,m}\sin n\mathit{\theta }& \text{for odd}m.\end{array}\right.\end{array}$$

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 F_m(θ) as

$$\begin{array}{}\text{(6)}& F_m\left(\mathit{\theta }\right)\cong \left\{\begin{array}{ll}\sum _{n=\mathrm{0}}^N{F}_{n,m}\cos n\mathit{\theta }& \text{for}m=\mathrm{0},\\ \sum _{n=\mathrm{1}}^{N^{\prime }}{F}_{n,m}\sin n\mathit{\theta }& \text{for\hspace{0.17em} odd}m,\\ \sum _{n=\mathrm{1}}^{N^{\prime }}{F}_{n,m}^{\prime }\sin \mathit{\theta }\sin n\mathit{\theta }& \text{for even}m\left(\ne \mathrm{0}\right).\end{array}\right.\end{array}$$

The meridional basis functions sin θsin nθ for even m (≠0) are different from Eq. (5). The coefficients $F_{n,m}^{\prime }\text{ for\hspace{0.17em} even }m\left(\ne \mathrm{0}\right)$ are calculated by the forward Fourier sine transform of $F_m\left(\mathit{\theta }\right)/\sin \mathit{\theta }$. 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(N²) operations to smooth out the high-wavenumber components (Cheong, 2000b; Cheong et al., 2002; Kwon et al., 2004). The semi-implicit Eulerian hydrostatic atmospheric model also ran stably with high-order horizontal diffusion (Cheong, 2006; Koo and Hong, 2013; Park et al., 2013). However, the computational results of these models appear to be a little different from (i.e., slightly worse than) the models using SH. One reason for this seems to be the appearance of high-wavenumber oscillation resulting from the meridional wavenumber truncation with $N=N^{\prime }\cong \mathrm{2}J/\mathrm{3}$ or $J/\mathrm{2}$ for even m (≠0) (see Sect. 4) and the use of high-order horizontal diffusion to smooth out the oscillation, where J is the number of grid points in the latitudinal direction.

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

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

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

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

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

1. High-wavenumber noise appears near the poles.

2. The meridional truncation wavenumber N^′ needs to be equal to J for even m (≠0) because $N^{\prime }<J$ (e.g., $N^{\prime }\cong \mathrm{2}J/\mathrm{3})$ for even m (≠0) causes the high-wavenumber oscillation (see Sect. 4) and numerical instability.

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

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

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

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

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

2 Arrangement of equally spaced latitudinal grid points

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

$$\begin{array}{}\text{(7a)}& {\displaystyle }\begin{array}{rl}\text{Grid}\left[\mathrm{0}\right]& :J=J^{\mathrm{0}},{\mathit{\theta }}_j=\mathit{\pi }\left(j+\mathrm{0.5}\right)/J^{\mathrm{0}},\\ & j=\mathrm{0},\mathrm{\dots },J^{\mathrm{0}}-\mathrm{1},\end{array}\text{(7b)}& {\displaystyle }\begin{array}{rl}\text{Grid}\left[\mathrm{1}\right]& :J=J^{\mathrm{0}}+\mathrm{1},{\mathit{\theta }}_j=\mathit{\pi }j/J^{\mathrm{0}},\\ & j=\mathrm{0},\mathrm{\dots },J^{\mathrm{0}},\end{array}\text{(7c)}& {\displaystyle }\begin{array}{rl}\text{Grid}\left[-\mathrm{1}\right]& :J=J^{\mathrm{0}}-\mathrm{1},{\mathit{\theta }}_j=\mathit{\pi }j/J^{\mathrm{0}},\\ & j=\mathrm{1},\mathrm{\dots },J^{\mathrm{0}}-\mathrm{1},\end{array}\end{array}$$

where θ_j is the colatitude at each grid point and J⁰ 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 J⁰+1, and the number of grid points J in Grid [−1] is J⁰−1. Figure 1 shows Grids [0], [1], and [−1] when J⁰=4 and the grid interval $\mathrm{\Delta }\mathit{\theta }=\mathit{\pi }/\mathrm{4}$. Grid [0] has been widely used in DFS models (e.g., Merilees, 1973b; Orszag, 1974; Cheong, 2000a, b; Yoshimura and Matsumura, 2005) and in DFS expansion (e.g., Cheong et al., 2004). Grid [1] was used in DFS expansion (e.g., Yee, 1981; Cheong et al., 2004). Grid [−1] was used, for example, in the SH model using Clenshaw–Curtis quadrature (Hotta and Ujiie, 2018). All of Grids [0], [1], and [−1] were used in SH expansion (Swarztrauber and Spotz, 2000).

Figure 1 Grid [0], Grid[1], and Grid [−1] are three ways of arranging equally spaced latitudinal grid points. Red circles show the positions of the grid points when the grid interval $\mathrm{\Delta }\mathit{\theta }=\mathit{\pi }/\mathrm{4}$.

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

3 Improved double Fourier series on the sphere

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

3.1 New basis functions for a scalar variable

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

$$\begin{array}{}\text{(8)}& T\left(\mathit{\lambda },\mathit{\theta }\right)\cong \sum _{m=\mathrm{0}}^M{T}_m^{\mathrm{c}}\left(\mathit{\theta }\right)\cos m\mathit{\lambda }+\sum _{m=\mathrm{1}}^M{T}_m^{\mathrm{s}}\left(\mathit{\theta }\right)\sin m\mathit{\lambda },\end{array}$$

where $T_m^{\mathrm{c}}\left(\mathit{\theta }\right)$ and $T_m^{\mathrm{s}}\left(\mathit{\theta }\right)$ are calculated from T(λ,θ) by the forward Fourier transform as

$$\begin{array}{}\text{(9a)}& {\displaystyle }{T}_m^{\mathrm{c}}\left(\mathit{\theta }\right)={\displaystyle \frac{a}{\mathrm{2}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\cos m\mathit{\lambda }T\left(\mathit{\lambda },\mathit{\theta }\right)\mathrm{d}\mathit{\lambda },a\equiv \left\{\begin{array}{l}\mathrm{1}\text{ for }m=\mathrm{0}\\ \mathrm{2}\text{ for }m\ge \mathrm{1},\end{array}\right.\text{(9b)}& {\displaystyle }{T}_m^{\mathrm{s}}\left(\mathit{\theta }\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\sin m\mathit{\lambda }T\left(\mathit{\lambda },\mathit{\theta }\right)\mathrm{d}\mathit{\lambda }.\end{array}$$

The variables $T_m^{\mathrm{c}}\left(\mathit{\theta }\right)$ and $T_m^{\mathrm{s}}\left(\mathit{\theta }\right)$ are meridionally expanded as

$$\begin{array}{}\text{(10)}& T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\cong T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\equiv \sum _{n=N_{\text{min,}m}}^{N_{\text{max,}m}}{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}{S}_{n,m}\left(\mathit{\theta }\right),\end{array}$$

where

$$\begin{array}{}\text{(11)}& {\displaystyle }{S}_{n,m}\left(\mathit{\theta }\right)\equiv \left\{\begin{array}{ll}\cos n\mathit{\theta }& \text{for}m=\mathrm{0},\\ \sin \mathit{\theta }\cos n\mathit{\theta }& \text{for}m=\mathrm{1},\\ \sin \mathit{\theta }\sin n\mathit{\theta }& \text{for even}m\ge \mathrm{2},\\ {\sin }^{\mathrm{2}}\mathit{\theta }\sin n\mathit{\theta }& \text{for odd}m\ge \mathrm{3},\end{array}\right.\text{(12)}& {\displaystyle }\begin{array}{rl}{N}_{\text{min,}m}& =\left\{\begin{array}{ll}\mathrm{0}& \text{for}m=\mathrm{0},\\ \mathrm{0}& \text{for}m=\mathrm{1},\\ \mathrm{1}& \text{for even}m\ge \mathrm{2},\\ \mathrm{1}& \text{for odd}m\ge \mathrm{3},\end{array}\right.\\ N_{\text{max,}m}& =\left\{\begin{array}{ll}N& \text{for}m=\mathrm{0},\\ N-\mathrm{1}& \text{for}m=\mathrm{1},\\ N-\mathrm{1}& \text{for even}m\ge \mathrm{2},\\ N-\mathrm{2}& \text{for odd}m\ge \mathrm{3}.\end{array}\right.\end{array}\end{array}$$

Here, the superscript c(s) means c or s, and thus, for example, $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ means $T_m^{\mathrm{c}}\left(\mathit{\theta }\right)$ or $T_m^{\mathrm{s}}\left(\mathit{\theta }\right)$. In Eq. (8), cos mλ and sin mλ are used instead of e^imλ 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 ²θsin nθ for odd m≥3 are especially different from the basis functions of Cheong in Eq. (6). Either sin nθ or sin θcos nθ can be used as the basis functions for m=1 because it can be shown using Eq. (A2a) from Appendix A that $\sin \mathit{\theta }\cos n\mathit{\theta }\left(n=\mathrm{0},\mathrm{\dots },N-\mathrm{1}\right)$ are the linear combination of $\sin n\mathit{\theta }\left(n=\mathrm{1},\mathrm{\dots },N\right)$ 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:

$$\begin{array}{}\text{(13)}& T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)=\left\{\begin{array}{ll}\sum _{n=\mathrm{0}}^N{T}_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}\cos n\mathit{\theta }& \text{for\hspace{0.17em} even}m,\\ \sum _{n=\mathrm{1}}^N{T}_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}\sin n\mathit{\theta }& \text{for\hspace{0.17em} odd}m,\end{array}\right.\end{array}$$

where

$$\begin{array}{}\text{(14a)}& {\displaystyle }{T}_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}=T_{n,m}^{\text{c(s)}}\left(n=\mathrm{0},\mathrm{\dots },N\right)\text{for\hspace{0.17em} }m=\mathrm{0},\text{(14b)}& {\displaystyle }\begin{array}{r}{T}_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}={\displaystyle \frac{T_{n-\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{n+\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}}{\mathrm{2}}}\left(n=\mathrm{1},\mathrm{\dots },N\right)\text{ for }m=\mathrm{1}\\ \text{except\hspace{0.17em} for}{T}_{\mathrm{1},m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}={\displaystyle \frac{\mathrm{2}{T}_{\mathrm{0},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}}{\mathrm{2}}}\left(n=\mathrm{1}\right),\end{array}\text{(14c)}& {\displaystyle }\begin{array}{rl}{T}_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}& ={\displaystyle \frac{-T_{n-\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}+T_{n+\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}}{\mathrm{2}}}\left(n=\mathrm{0},\mathrm{\dots },N\right)\\ & \text{for\hspace{0.17em} even}m\ge \mathrm{2},\end{array}\text{(14d)}& {\displaystyle }\begin{array}{rl}{T}_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}& ={\displaystyle \frac{-T_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{2}{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}}{\mathrm{4}}}\left(n=\mathrm{1},\mathrm{\dots },N\right)\\ & \text{for odd}m\ge \mathrm{3}\\ & \text{except\hspace{0.17em} for}{T}_{\mathrm{1},m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}={\displaystyle \frac{\mathrm{3}{T}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}}{\mathrm{4}}}\left(n=\mathrm{1}\right).\end{array}\end{array}$$

The value of N_max,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 J⁰−1 for each m. In Grid [−1], the upper limit of N is J⁰−1 for m≥2, but $J^{\mathrm{0}}-\mathrm{2}\left(=J-\mathrm{1}\right)$ for m= 0 or 1. The reason for this is shown in Appendix C.

When calculating the values of $T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ in grid space from $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ in spectral space, the coefficients $T_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}$ are calculated from $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ using Eq. (14), and inverse discrete cosine and sine transforms (see Appendix B) are performed using Eq. (13). The calculation of $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ in spectral space from $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ in grid space is described in Sect. 3.3 below.

The truncated variable $T^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ is defined as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{T}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)& \equiv \sum _{m=\mathrm{0}}^M{T}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\cos m\mathit{\lambda }\\ & +\sum _{m=\mathrm{1}}^M{T}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\sin m\mathit{\lambda }.\end{array}\end{array}$$

From Eq. (10), the values of $T_m^{\mathrm{c},N}\left(\mathit{\theta }\right)$ at the poles are finite for m=0, and the values of $T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ at the poles are zero for m≠0. Therefore, $T^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ is continuous at the poles.

3.2 Gradient of a scalar variable

The gradient $\mathrm{\nabla }{T}^{N,M}=\left(T_{\mathit{\lambda }}^{N,M},T_{\mathit{\varphi }}^{N,M}\right)$ is obtained as follows:

$$\begin{array}{}\text{(16a)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda }}^{N,M}& \equiv {\displaystyle \frac{\mathrm{1}}{a\sin \mathit{\theta }}}{\displaystyle \frac{\partial T^{N,M}}{\partial \mathit{\lambda }}}=\sum _{m=\mathrm{1}}^M{T}_{\mathit{\lambda },m}^{\mathrm{c},N}\left(\mathit{\theta }\right)\cos m\mathit{\lambda }\\ & +\sum _{m=\mathrm{1}}^M{T}_{\mathit{\lambda },m}^{\mathrm{s},N}\left(\mathit{\theta }\right)\sin m\mathit{\lambda },\end{array}\text{(16b)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda },m}^{\mathrm{c},N}\left(\mathit{\theta }\right)& \equiv {\displaystyle \frac{m}{a\sin \mathit{\theta }}}{T}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\\ & =\sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\left(T_{n,m}^{\mathrm{s}}{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\right),\end{array}\text{(16c)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda },m}^{\mathrm{s},N}\left(\mathit{\theta }\right)& \equiv -{\displaystyle \frac{m}{a\sin \mathit{\theta }}}{T}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\\ & =\sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\left(-T_{n,m}^{\mathrm{c}}{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(17a)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\varphi }}^{N,M}& \equiv {\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial T^{N,M}}{\partial \mathit{\varphi }}}=-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial T^{N,M}}{\partial \mathit{\theta }}}\\ & =\sum _{m=\mathrm{0}}^M{T}_{\mathit{\varphi },m}^{\mathrm{c},N}\left(\mathit{\theta }\right)\cos m\mathit{\lambda }+\sum _{m=\mathrm{1}}^M{T}_{\mathit{\varphi },m}^{\mathrm{s},N}\left(\mathit{\theta }\right)\sin m\mathit{\lambda },\end{array}\text{(17b)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\varphi },m}^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)& \equiv -{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\\ & =\sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\left(-T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\right),\end{array}\end{array}$$

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

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& T_{\mathit{\lambda },m}^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)=\\ & \left\{\begin{array}{l}\mathrm{0}\text{for}m=\mathrm{0},\\ \sum _{n=\mathrm{0}}^{N-\mathrm{1}}{T}_{\mathit{\lambda },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\cos n\mathit{\theta }\text{for}m=\mathrm{1},\\ \sum _{n=\mathrm{1}}^{N-\mathrm{1}}{T}_{\mathit{\lambda },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\sin n\mathit{\theta }\text{for\hspace{0.17em} even}m\ge \mathrm{2},\\ \sum _{n=\mathrm{0}}^{N-\mathrm{1}}{T}_{\mathit{\lambda },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\cos n\mathit{\theta }\left(=\sum _{n=\mathrm{1}}^{N-\mathrm{2}}{T}_{\mathit{\lambda },n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}\sin \mathit{\theta }\sin n\mathit{\theta }\right)\\ \text{for odd}m\ge \mathrm{3},\end{array}\right.\end{array}\end{array}$$

where

$$\begin{array}{}\text{(19a)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda },n,m}^{\mathrm{c}}& ={\displaystyle \frac{\mathrm{1}}{a}}m{T}_{n,m}^{\mathrm{s}}\left(n=\mathrm{0},\mathrm{\dots },N-\mathrm{1}\right)\\ & \text{ for }m=\mathrm{1},\end{array}\text{(19b)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda },n,m}^{\mathrm{c}}& ={\displaystyle \frac{\mathrm{1}}{a}}m{T}_{n,m}^{\mathrm{s}}\left(n=\mathrm{1},\mathrm{\dots },N-\mathrm{1}\right)\\ & \text{for\hspace{0.17em} even}m\ge \mathrm{2},\end{array}\text{(19c)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda },n,m}^{\mathrm{c}}& ={\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{m\left(-T_{n-\mathrm{1},m}^{\mathrm{s}}+T_{n+\mathrm{1},m}^{\mathrm{s}}\right)}{\mathrm{2}}}\left(n=\mathrm{0},\mathrm{\dots },N-\mathrm{1}\right)\\ & \text{for\hspace{0.17em} odd}m\ge \mathrm{3}.\end{array}\end{array}$$

The equations for $T_{\mathit{\lambda },n,m}^{\mathrm{s}}$ are the same as Eq. (19), except that $T_{\mathit{\lambda },n,m}^{\mathrm{c}}$ and $T_{n,m}^{\mathrm{s}}$ are replaced with $T_{\mathit{\lambda },n,m}^{\mathrm{s}}$ and $-T_{n,m}^{\mathrm{c}}$, respectively. From Eqs. (17b), (13), and (14), we obtain

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& T_{\mathit{\varphi },m}^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)=\\ & \left\{\begin{array}{l}\sum _{n=\mathrm{1}}^N{T}_{\mathit{\varphi },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\sin n\mathit{\theta }\text{ for }m=\mathrm{0},\\ \sum _{n=\mathrm{0}}^N{T}_{\mathit{\varphi },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\cos n\mathit{\theta }\text{ for }m=\mathrm{1},\\ \sum _{n=\mathrm{1}}^N{T}_{\mathit{\varphi },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\sin n\mathit{\theta }\text{for\hspace{0.17em} even}m\ge \mathrm{2},\\ \sum _{n=\mathrm{0}}^N{T}_{\mathit{\varphi },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\cos n\mathit{\theta }\left(=\sum _{n=\mathrm{1}}^{N-\mathrm{1}}{T}_{\mathit{\varphi },n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}\sin \mathit{\theta }\sin n\mathit{\theta }\right)\\ \text{for\hspace{0.17em} odd}m\ge \mathrm{3},\end{array}\right.\end{array}\end{array}$$

where

$$\begin{array}{}\text{(21a)}& {\displaystyle }{T}_{\mathit{\varphi },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}=-{\displaystyle \frac{\mathrm{1}}{a}}\left(-n{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\right)\left(n=\mathrm{1},\mathrm{\dots },N\right)\text{ for }m=\mathrm{0},\text{(21b)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\varphi },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}& =-{\displaystyle \frac{\mathrm{1}}{a}}\left[{\displaystyle \frac{n\left(T_{n-\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{n+\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right)}{\mathrm{2}}}\right]\\ & \left(n=\mathrm{0},\mathrm{\dots },N\right)\text{ for }m=\mathrm{1},\\ & \text{except\hspace{0.17em} for}{T}_{\mathit{\varphi },\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}=-{\displaystyle \frac{\mathrm{1}}{a}}\left[{\displaystyle \frac{\mathrm{2}{T}_{\mathrm{0},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}}{\mathrm{2}}}\right]\\ & \left(n=\mathrm{1}\right),\end{array}\text{(21c)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\varphi },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}& =-{\displaystyle \frac{\mathrm{1}}{a}}\left[{\displaystyle \frac{n\left(T_{n-\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{n+\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right)}{\mathrm{2}}}\right]\\ & \left(n=\mathrm{1},\mathrm{\dots },N\right)\text{for\hspace{0.17em} even }m\ge \mathrm{2},\end{array}\text{(21d)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\varphi },n,m}^{\mathrm{c}\left(\mathrm{s}\right)}& =-{\displaystyle \frac{\mathrm{1}}{a}}\left[{\displaystyle \frac{n\left(-T_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{2}{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right)}{\mathrm{4}}}\right]\\ & \left(n=\mathrm{0},\mathrm{\dots },N\right)\text{for\hspace{0.17em} odd}m\ge \mathrm{3},\\ & \text{except for}{T}_{\mathit{\varphi },\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}=-{\displaystyle \frac{\mathrm{1}}{a}}\left[{\displaystyle \frac{\left(\mathrm{3}{T}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right)}{\mathrm{4}}}\right]\\ & \left(n=\mathrm{1}\right).\end{array}\end{array}$$

From Eqs. (18)–(21), it can be seen that $T_{\mathit{\lambda },m}^{\mathrm{c},N}\left(\mathit{\theta }\right)$, $T_{\mathit{\lambda },m}^{\mathrm{s},N}\left(\mathit{\theta }\right)$, $T_{\mathit{\varphi },m}^{\mathrm{c},N}\left(\mathit{\theta }\right)$, and $T_{\mathit{\varphi },m}^{\mathrm{s},N}\left(\mathit{\theta }\right)$ at the poles are finite for m=1 and zero for m≠1, and moreover the following relations are satisfied for m=1:

$$\begin{array}{}\text{(22a)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda },m=\mathrm{1}}^{\mathrm{c},N}\left(\mathit{\theta }\right)& =-T_{\mathit{\varphi },m=\mathrm{1}}^{\mathrm{s},N}\left(\mathit{\theta }\right)\left(={\displaystyle \frac{\mathrm{1}}{a}}\sum _{n=\mathrm{0}}^{N-\mathrm{1}}{T}_{n,m=\mathrm{1}}^{\mathrm{s}}\right)\\ & \text{at\hspace{0.17em} }\mathit{\theta }=\mathrm{0}\left(\text{North Pole}\right),\end{array}\text{(22b)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda },m=\mathrm{1}}^{\mathrm{s},N}\left(\mathit{\theta }\right)& =T_{\mathit{\varphi },m=\mathrm{1}}^{\mathrm{c},N}\left(\mathit{\theta }\right)\left(=-{\displaystyle \frac{\mathrm{1}}{a}}\sum _{n=\mathrm{0}}^{N-\mathrm{1}}{T}_{n,m=\mathrm{1}}^{\mathrm{c}}\right)\\ & \text{at\hspace{0.17em} }\mathit{\theta }=\mathrm{0}\left(\text{North Pole}\right),\end{array}\text{(22c)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda },m=\mathrm{1}}^{\mathrm{c},N}\left(\mathit{\theta }\right)& =T_{\mathit{\varphi },m=\mathrm{1}}^{\mathrm{s},N}\left(\mathit{\theta }\right)\left(={\displaystyle \frac{\mathrm{1}}{a}}\sum _{n=\mathrm{0}}^{N-\mathrm{1}}{\left(-\mathrm{1}\right)}^n{T}_{n,m=\mathrm{1}}^{\mathrm{s}}\right)\\ & \text{at}\mathit{\theta }=\mathit{\pi }\left(\text{South Pole}\right),\end{array}\text{(22d)}& {\displaystyle }\begin{array}{rl}{T}_{\mathit{\lambda },m=\mathrm{1}}^{\mathrm{s},N}\left(\mathit{\theta }\right)& =-T_{\mathit{\varphi },m=\mathrm{1}}^{\mathrm{c},N}\left(\mathit{\theta }\right)\left(=-{\displaystyle \frac{\mathrm{1}}{a}}\sum _{n=\mathrm{0}}^{N-\mathrm{1}}{\left(-\mathrm{1}\right)}^n{T}_{n,m=\mathrm{1}}^{\mathrm{c}}\right)\\ & \text{at\hspace{0.17em} }\mathit{\theta }=\mathit{\pi }\left(\text{South\hspace{0.17em} Pole}\right).\end{array}\end{array}$$

Thus, it is guaranteed that $\mathrm{\nabla }{T}^{N,M}=\left(T_{\mathit{\lambda }}^{N,M},T_{\mathit{\varphi }}^{N,M}\right)$ is continuous at the poles.

3.3 New method to calculate expansion coefficients for a scalar variable

One way to calculate the coefficients $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ from $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ in Eq. (10) is to perform a forward cosine transform of $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)/\sin \mathit{\theta }$ for m=1, a sine transform of $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)/\sin \mathit{\theta }$ for even m (≥2), and a sine transform of $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)/{\sin }^{\mathrm{2}}\mathit{\theta }$ for odd m (≥3). However, this approach with the meridional wavenumber truncation N<J leads to large high-wavenumber oscillations (see Sect. 4). Dividing $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ by sin ²θ reduces the numerical stability of the model more significantly than dividing $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ 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 $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ by sin θ or sin ²θ before the forward cosine or sine transforms. The coefficients $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ in Eq. (10) are calculated as follows. First, $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ in Eq. (10) are expanded like Eq. (5) as

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{T}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)& \cong {\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\\ & \equiv \left\{\begin{array}{ll}\sum _{n=\mathrm{0}}^N{\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\cos n\mathit{\theta }& \text{for\hspace{0.17em} even}m,\\ \sum _{n=\mathrm{1}}^N{\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\sin n\mathit{\theta }& \text{for\hspace{0.17em} odd}m,\end{array}\right.\end{array}\end{array}$$

where the expansion coefficients ${\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ are calculated by the forward discrete cosine transform for even m and the forward discrete sine transform for odd m from the values of $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ at the grid points (see Appendix B).

Next, $T_{n,m}^{\mathrm{c}}$ and $T_{n,m}^{\mathrm{s}}$ are calculated using the least-squares method to minimize the following error E (the squared L₂ norm of the residual):

$$\begin{array}{}\text{(24)}& E\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}R{\left(\mathit{\lambda },\mathit{\theta }\right)}^{\mathrm{2}}\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda },\end{array}$$

where the residual R(λ,θ) is

$$\begin{array}{}\text{(25)}& R\left(\mathit{\lambda },\mathit{\theta }\right)\equiv T^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)-T\left(\mathit{\lambda },\mathit{\theta }\right).\end{array}$$

From Eqs. (24), (25), and (15) and the equations $\partial E/\partial T_{n,m}^{\mathrm{c}}=\mathrm{0}$ and $\partial E/\partial T_{n,m}^{\mathrm{s}}=\mathrm{0}$ used in the least-squares method to minimize E, we obtain

$$\begin{array}{}\text{(26a)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{\displaystyle \frac{\partial T_m^{\mathrm{c},N}\left(\mathit{\theta }\right)}{\partial T_{n,m}^{\mathrm{c}}}}\cos m\mathit{\lambda }R\left(\mathit{\lambda },\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }=\mathrm{0},\text{(26b)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{\displaystyle \frac{\partial T_m^{\mathrm{s},N}\left(\mathit{\theta }\right)}{\partial T_{n,m}^{\mathrm{s}}}}\sin m\mathit{\lambda }R\left(\mathit{\lambda },\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }=\mathrm{0}.\end{array}$$

From Eq. (10), we derive

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\partial T_m^{\mathrm{c},N}\left(\mathit{\theta }\right)}{\partial T_{n,m}^{\mathrm{c}}}}={\displaystyle \frac{\partial T_m^{\mathrm{s},N}\left(\mathit{\theta }\right)}{\partial T_{n,m}^{\mathrm{s}}}}=S_{n,m}\left(\mathit{\theta }\right).\end{array}$$

Equations (26) and (27) show that the residual R(λ,θ) is orthogonal to each of the new DFS basis functions S_m,n(θ)cos mλ and S_m,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

$$\begin{array}{}\text{(28)}& \underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}\left(\mathit{\theta }\right)T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}\left(\mathit{\theta }\right)T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }.\end{array}$$

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

$$\begin{array}{}\text{(29)}& \underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}\left(\mathit{\theta }\right)T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}\left(\mathit{\theta }\right){\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }.\end{array}$$

By substituting Eqs. (10) and (23) into Eq. (29), the following equations for $T_{n,m}^{\mathrm{c}}$ and $T_{n,m}^{\mathrm{s}}$ are derived, as shown in Appendix E.

For m=0,

$$\begin{array}{}\text{(30a)}& T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}={\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathrm{0}\le n\le N\right).\end{array}$$

For m=1,

$$\begin{array}{}\text{(30b)}& \begin{array}{rl}-T_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{2}{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}& =-\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{n-\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{n+\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & \left(\mathrm{0}\le n\le N-\mathrm{1}\right),\end{array}\end{array}$$

with the exception of the following underlined values:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{¯}}{\mathrm{1}}{T}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}=\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\left(n=\mathrm{1}\right),\\ {\displaystyle }& {\displaystyle }-\underset{\mathrm{¯}}{\mathrm{2}}{T}_{\mathrm{0},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{2}{T}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}=-\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}\left(n=\mathrm{2}\right).\end{array}$$

For even m (≥2),

$$\begin{array}{}\text{(30c)}& \begin{array}{rl}-T_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{2}{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}& =\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{n-\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{n+\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & \left(\mathrm{1}\le n\le N-\mathrm{1}\right),\end{array}\end{array}$$

with the exception of the following underlined values:

$$\underset{\mathrm{¯}}{\mathrm{3}}{T}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-T_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}=\underset{\mathrm{¯}}{\mathrm{4}}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{0},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\left(n=\mathrm{1}\right).$$

For odd m(≥3),

$$\begin{array}{}\text{(30d)}& \begin{array}{rl}{T}_{n-\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}& -\mathrm{4}{T}_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{6}{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{T}_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}+T_{n+\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & =-\mathrm{4}{\stackrel{\mathrm{̃}}{T}}_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{8}{\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{\stackrel{\mathrm{̃}}{T}}_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & \left(\mathrm{1}\le n\le N-\mathrm{2}\right),\end{array}\end{array}$$

with the exception of the following underlined values:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{¯}}{\mathrm{10}}{T}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\underset{\mathrm{¯}}{\mathrm{5}}{T}_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}+T_{\mathrm{5},m}^{\mathrm{c}\left(\mathrm{s}\right)}=\underset{\mathrm{¯}}{\mathrm{12}}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}\left(n=\mathrm{1}\right),\\ {\displaystyle }& {\displaystyle }\underset{\mathrm{¯}}{\mathrm{5}}{T}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{T}_{\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}+T_{\mathrm{6},m}^{\mathrm{c}\left(\mathrm{s}\right)}=\mathrm{8}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}\left(n=\mathrm{2}\right),\\ {\displaystyle }& {\displaystyle }-\underset{\mathrm{¯}}{\mathrm{5}}{T}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{6}{T}_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{T}_{\mathrm{5},m}^{\mathrm{c}\left(\mathrm{s}\right)}+T_{\mathrm{7},m}^{\mathrm{c}\left(\mathrm{s}\right)}=-\mathrm{4}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{8}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\displaystyle }& {\displaystyle }-\mathrm{4}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{5},m}^{\mathrm{c}\left(\mathrm{s}\right)}\left(n=\mathrm{3}\right).\end{array}$$

From Eq. (30a), $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ for m=0 is obtained. From Eqs. (30d), the following linear simultaneous equations for m≥3 are derived:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& {\mathbf{C}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}\left[\begin{array}{l}{T}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ T_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ T_{\mathrm{5},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ T_{\mathrm{7},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ \mathrm{⋮}\end{array}\right]={\mathbf{D}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}\left[\begin{array}{l}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\stackrel{\mathrm{̃}}{T}}_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\stackrel{\mathrm{̃}}{T}}_{\mathrm{5},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\stackrel{\mathrm{̃}}{T}}_{\mathrm{7},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ \mathrm{⋮}\end{array}\right],\\ & {\mathbf{C}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}\left[\begin{array}{l}{T}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ T_{\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ T_{\mathrm{6},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ T_{\mathrm{8},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ \mathrm{⋮}\end{array}\right]={\mathbf{D}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}\left[\begin{array}{l}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\stackrel{\mathrm{̃}}{T}}_{\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\stackrel{\mathrm{̃}}{T}}_{\mathrm{6},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\stackrel{\mathrm{̃}}{T}}_{\mathrm{8},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ \mathrm{⋮}\end{array}\right],\end{array}\end{array}$$

where the matrices C_{n_odd,m} and C_{n_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 C_{n_odd,m} and C_{n_even,m} are derived. By using Eq. (31), the expansion coefficients $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ are calculated from ${\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$. A penta-diagonal matrix C can be lower–upper (LU) decomposed as

$$\begin{array}{}\text{(32)}& \begin{array}{rl}\mathbf{C}=\mathbf{LU},\mathbf{L}& \equiv \left[\begin{array}{cccccccc}\ast & \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \ast & \ast & \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \ast & \ast & \ast & \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \ast & \ast & \ast & \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ & & & & \mathrm{⋮}& & & \\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}& \ast & \ast & \ast & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}& \mathrm{0}& \ast & \ast & \ast \end{array}\right],\\ \mathbf{U}& \equiv \left[\begin{array}{cccccccc}\mathrm{1}& \ast & \ast & \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \ast & \ast & \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \ast & \ast & \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \ast & \ast & \mathrm{\cdots }& \mathrm{0}\\ & & & & \mathrm{⋮}& & & \\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \ast \\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right].\end{array}\end{array}$$

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 $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ for all m and n with Eq. (30) requires only O(N²) 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, $T_m^{\mathrm{c}}\left(\mathit{\theta }\right)$ and$T_m^{\mathrm{s}}\left(\mathit{\theta }\right)$ in Eq. (8) are expanded with the associated Legendre functions P_n,m(θ) as

$$\begin{array}{}\text{(33)}& T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\cong T_m^{\mathrm{c}\left(\mathrm{s}\right)\text{,SH},N}\left(\mathit{\theta }\right)\equiv \sum _{n=m}^N{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right),\mathrm{SH}}{P}_{n,m}\left(\mathit{\theta }\right),\end{array}$$

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

$$\begin{array}{}\text{(34)}& \underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{P}_{n,m}\left(\mathit{\theta }\right)P_{n^{\prime },m}\left(\mathit{\theta }\right)\sin \mathit{\theta }\mathrm{d}\mathit{\theta }=\left\{\begin{array}{ll}\mathrm{1}\left(\text{or}\mathrm{2}\right)& \text{for}n=n^{\prime },\\ \mathrm{0}& \text{for}n\ne n^{\prime }.\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(35)}& P_{n,m}\left(\mathit{\theta }\right)\text{=}\sum _{\begin{array}{c}l=\mathrm{0}\\ \text{when}n-\left|m\right|-l\text{is\hspace{0.17em} even}\end{array}}^{n-\left|m\right|}{a}_{n,m,l}{\sin }^{\left|m\right|}\mathit{\theta }\cos l\mathit{\theta }.\end{array}$$

Conversely, the functions ${\sin }^{\left|m\right|}\mathit{\theta }\cos \left(n-\left|m\right|\right)\mathit{\theta }\left(n\ge \left|m\right|\right)$ can be expressed as the linear combination of $P_{l,m}\left(\mathit{\theta }\right)\left(l=\left|m\right|,\mathrm{\dots },n\right)$. Substituting Eq. (35) into Eq. (33) gives

$$\begin{array}{}\text{(36)}& T_m^{\mathrm{c}\left(\mathrm{s}\right),\mathrm{SH},N}\left(\mathit{\theta }\right)=\sum _{n=\mathrm{0}}^{N-m}{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right),{\text{SH}}^{\prime }}{\sin }^m\mathit{\theta }\cos n\mathit{\theta },\end{array}$$

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 ${\sin }^{\mathrm{2}}\mathit{\theta }\cos n\mathit{\theta }\left(n=\mathrm{0},\mathrm{\dots },N-\mathrm{2}\right)$ for m=2 and ${\sin }^{\mathrm{3}}\mathit{\theta }\cos n\mathit{\theta }\left(n=\mathrm{0},\mathrm{\dots },N-\mathrm{3}\right)$ for m=3 in Eq. (36) are the linear combinations of $\sin \mathit{\theta }\sin n\mathit{\theta }\left(n=\mathrm{1},\mathrm{\dots },N-\mathrm{1}\right)$ and ${\sin }^{\mathrm{2}}\mathit{\theta }\sin n\mathit{\theta }\left(n=\mathrm{1},\mathrm{\dots },N-\mathrm{2}\right)$ 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 $T_{n,m}^{\mathrm{c},\text{SH}}$ in the SH model is about $N^{\mathrm{2}}/\mathrm{2}$ when M=N. This triangular truncation used in the SH method gives a uniform resolution over the sphere. From Eqs. (8) and (10), the number of the expansion coefficients $T_{n,m}^{\mathrm{c}}$ in the DFS method is about N² 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 $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right),\mathrm{SH}}$ in Eq. (33) are calculated from $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ by the forward Legendre transform as

$$\begin{array}{}\text{(37)}& T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right),\mathrm{SH}}=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{P}_{n,m}\left(\mathit{\theta }\right)T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\sin \mathit{\theta }\mathrm{d}\mathit{\theta },\end{array}$$

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 $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ is replaced with ${\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ (e.g., Sneeuw and Bun, 1996), although the values of $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right),\mathrm{SH}}$ calculated from ${\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ are different from those calculated from $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$. In the new DFS method, the values of $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ calculated from ${\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ in Eq. (29) are the same as those calculated from ${\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ in Eq. (28) (see Eq. D4 in Appendix D).

Equation (37) can be derived using the least-squares method that minimizes the error E^SH (the squared L₂ norm of the residual):

$$\begin{array}{}\text{(38)}& E^{\text{SH}}\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{R}^{\text{SH}}{\left(\mathit{\lambda },\mathit{\theta }\right)}^{\mathrm{2}}\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda },\end{array}$$

where the residual R^SH(λ,θ) is

$$\begin{array}{}\text{(39)}& \begin{array}{rl}{R}^{\text{SH}}\left(\mathit{\lambda },\mathit{\theta }\right)& \equiv \left(\sum _{m=\mathrm{0}}^M{T}_m^{\mathrm{c},\text{SH},N}\left(\mathit{\theta }\right)\cos m\mathit{\lambda }\right.\\ & \left.+\sum _{m=\mathrm{1}}^M{T}_m^{\mathrm{s},\text{SH},N}\left(\mathit{\theta }\right)\sin m\mathit{\lambda }\right)-T\left(\mathit{\lambda },\mathit{\theta }\right).\end{array}\end{array}$$

From Eqs. (38), (39), and (33) and the equations $\partial E^{\text{SH}}/\partial T_{n,m}^{\mathrm{c},\text{SH}}=\mathrm{0}$ and $\partial E^{\text{SH}}/\partial T_{n,m}^{\mathrm{s},\text{SH}}=\mathrm{0}$ used in the least-squares method to minimize E^SH, we derive

$$\begin{array}{}\text{(40a)}& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{P}_{n,m}\left(\mathit{\theta }\right)\cos m\mathit{\lambda }{R}^{\text{SH}}\left(\mathit{\lambda },\mathit{\theta }\right)\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }=\mathrm{0},\text{(40b)}& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{P}_{n,m}\left(\mathit{\theta }\right)\sin m\mathit{\lambda }{R}^{\text{SH}}\left(\mathit{\lambda },\mathit{\theta }\right)\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }=\mathrm{0}.\end{array}$$

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 $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ 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 $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ with dense matrices, which leads to O(N³) operations. When using E, the simultaneous equations with penta-diagonal or tri-diagonal matrices require only O(N²) operations. Therefore, we choose to use E instead of E^′.

The new DFS meridional basis functions S_n,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 $S_{n,m}^{\text{O}}\left(\mathit{\theta }\right)$, which satisfy

$$\begin{array}{}\text{(41)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}^{\text{O}}\left(\mathit{\theta }\right)S_{n^{\prime },m}^{\text{O}}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\left\{\begin{array}{ll}\mathrm{1}& \text{for}n=n^{\prime },\\ \mathrm{0}& \text{for}n\ne n^{\prime }.\end{array}\right.\end{array}$$

This is similar to Eq. (34), but the latitudinal weight is constant. $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ in Eq. (8) are expanded with $S_{n,m}^{\text{O}}\left(\mathit{\theta }\right)$ as

$$\begin{array}{}\text{(42)}& T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\cong T_m^{\mathrm{c}\left(\mathrm{s}\right),\text{O},N}\left(\mathit{\theta }\right)\equiv \sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}{T}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)\text{,O}}{S}_{n,m}^{\text{O}}\left(\mathit{\theta }\right).\end{array}$$

By using the least squares method or the Galerkin method with Eq. (42), we obtain the same equations as Eqs. (24)–(29), except that $T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ and S_n,m(θ) are replaced with $T_m^{\mathrm{c}\left(\mathrm{s}\right),\text{O},N}\left(\mathit{\theta }\right)$ and $S_{n,m}^{\text{O}}\left(\mathit{\theta }\right)$, respectively. From Eq. (29) with $T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ and S_n,m(θ) replaced by $T_m^{\mathrm{c}\left(\mathrm{s}\right),\text{O},N}\left(\mathit{\theta }\right)$ and $S_{n,m}^{\text{O}}\left(\mathit{\theta }\right)$ and Eqs. (41) and (42), we derive

$$\begin{array}{}\text{(43)}& T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)\text{,O}}={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}^{\text{O}}\left(\mathit{\theta }\right){\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }.\end{array}$$

Thus, $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)\text{,O}}$ and $T_m^{\mathrm{c}\left(\mathrm{s}\right),\text{O},N}\left(\mathit{\theta }\right)$ in Eqs. (43) and (42) are calculated uniquely. This unique solution $T_m^{\mathrm{c}\left(\mathrm{s}\right),\text{O},N}\left(\mathit{\theta }\right)$ is the same as $T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ in Eq. (29) obtained by the least-squares method with the non-orthogonal basis functions S_n,m(θ) because $S_{n,m}^{\text{O}}\left(\mathit{\theta }\right)\left(n=N_{\text{min,}m},\mathrm{\dots },N_{\text{max,\hspace{0.17em} }m}\right)$ is the linear combination of $S_{n,m}\left(\mathit{\theta }\right)\left(n=N_{\text{min,}m},\mathrm{\dots },N_{\text{max,\hspace{0.17em} }m}\right)$ 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

$$\begin{array}{}\text{(44a)}& {\displaystyle }u={\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial \mathit{\chi }}{\partial \mathit{\lambda }}}-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\varphi }}}={\displaystyle \frac{\mathrm{1}}{a\sin \mathit{\theta }}}{\displaystyle \frac{\partial \mathit{\chi }}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\theta }}},\text{(44b)}& {\displaystyle }v={\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial \mathit{\chi }}{\partial \mathit{\varphi }}}={\displaystyle \frac{\mathrm{1}}{a\sin \mathit{\theta }}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\lambda }}}-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial \mathit{\chi }}{\partial \mathit{\theta }}},\end{array}$$

where $u=a\cos \mathit{\varphi }\mathrm{d}\mathit{\lambda }/\mathrm{d}t$ is the zonal wind and $v=a\mathrm{d}\mathit{\varphi }/\mathrm{d}t$ is the meridional wind. The scalar variables χ and ψ are expanded like Eqs. (8) and (10) as

$$\begin{array}{}\text{(45)}& {\displaystyle }\begin{array}{rl}\left[\begin{array}{c}\mathit{\chi }\left(\mathit{\lambda },\mathit{\theta }\right)\\ \mathit{\psi }\left(\mathit{\lambda },\mathit{\theta }\right)\end{array}\right]& \cong \sum _{m=\mathrm{0}}^M\left[\begin{array}{c}{\mathit{\chi }}_m^{\mathrm{c}}\left(\mathit{\theta }\right)\\ {\mathit{\psi }}_m^{\mathrm{c}}\left(\mathit{\theta }\right)\end{array}\right]\cos m\mathit{\lambda }\\ & +\sum _{m=\mathrm{1}}^M\left[\begin{array}{c}{\mathit{\chi }}_m^{\mathrm{s}}\left(\mathit{\theta }\right)\\ {\mathit{\psi }}_m^{\mathrm{s}}\left(\mathit{\theta }\right)\end{array}\right]\sin m\mathit{\lambda },\end{array}\text{(46)}& {\displaystyle }\begin{array}{rl}\left[\begin{array}{c}{\mathit{\chi }}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\\ {\mathit{\psi }}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\end{array}\right]& \cong \left[\begin{array}{c}{\mathit{\chi }}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\\ {\mathit{\psi }}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\end{array}\right]\\ & \equiv \sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\left[\begin{array}{c}{\mathit{\chi }}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\mathit{\psi }}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\end{array}\right]S_{n,m}\left(\mathit{\theta }\right),\end{array}\end{array}$$

The truncated variables ${\mathit{\psi }}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ and ${\mathit{\chi }}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ are defined as

$$\begin{array}{}\text{(47)}& \begin{array}{rl}\left[\begin{array}{c}{\mathit{\chi }}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)\\ {\mathit{\psi }}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)\end{array}\right]& \equiv \sum _{m=\mathrm{0}}^M\left[\begin{array}{c}{\mathit{\chi }}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\\ {\mathit{\psi }}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\end{array}\right]\cos m\mathit{\lambda }\\ & +\sum _{m=\mathrm{1}}^M\left[\begin{array}{c}{\mathit{\chi }}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\\ {\mathit{\psi }}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\end{array}\right]\sin m\mathit{\lambda },\end{array}\end{array}$$

From Eqs. (44)–(47), the equations for the wind vector components $u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ and $v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ are derived as

$$\begin{array}{}\text{(48a)}& {\displaystyle }\begin{array}{rl}{u}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)& \equiv {\displaystyle \frac{\mathrm{1}}{a\sin \mathit{\theta }}}{\displaystyle \frac{\partial {\mathit{\chi }}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial {\mathit{\psi }}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial \mathit{\theta }}}\\ & =\sum _{m=\mathrm{0}}^M{u}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\cos m\mathit{\lambda }+\sum _{m=\mathrm{1}}^M{u}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\sin m\mathit{\lambda },\end{array}\text{(48b)}& {\displaystyle }\begin{array}{rl}{u}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)& \equiv {\displaystyle \frac{m{\mathit{\chi }}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial {\mathit{\psi }}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\\ & =\sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\left({\mathit{\chi }}_{n,m}^{\mathrm{s}}{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\right.\\ & \left.+{\mathit{\psi }}_{n,m}^{\mathrm{c}}{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\right),\end{array}\text{(48c)}& {\displaystyle }\begin{array}{rl}{u}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)& \equiv -{\displaystyle \frac{m{\mathit{\chi }}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial {\mathit{\psi }}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\\ & =\sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\left(-{\mathit{\chi }}_{n,m}^{\mathrm{c}}{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\right.\\ & \left.+{\mathit{\psi }}_{n,m}^{\mathrm{s}}{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(49a)}& {\displaystyle }\begin{array}{rl}{v}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)& \equiv {\displaystyle \frac{\mathrm{1}}{a\sin \mathit{\theta }}}{\displaystyle \frac{\partial {\mathit{\psi }}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial \mathit{\lambda }}}\\ & -{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial {\mathit{\chi }}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial \mathit{\theta }}}=\sum _{m=\mathrm{0}}^M{v}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\cos m\mathit{\lambda }\\ & +\sum _{m=\mathrm{1}}^M{v}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\sin m\mathit{\lambda },\end{array}\text{(49b)}& {\displaystyle }\begin{array}{rl}{v}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)& \equiv {\displaystyle \frac{m{\mathit{\psi }}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial {\mathit{\chi }}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\\ & =\sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\left({\mathit{\psi }}_{n,m}^{\mathrm{s}}{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\right.\\ & \left.-{\mathit{\chi }}_{n,m}^{\mathrm{c}}{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\right),\end{array}\text{(49c)}& {\displaystyle }\begin{array}{rl}{v}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)& \equiv -{\displaystyle \frac{m{\mathit{\psi }}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial {\mathit{\chi }}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\\ & =\sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\left(-{\mathit{\psi }}_{n,m}^{\mathrm{c}}{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\right.\\ & \left.-{\mathit{\chi }}_{n,m}^{\mathrm{s}}{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\right).\end{array}\end{array}$$

The vector $\left(u^{N,M},v^{N,M}\right)$ in Eqs. (48) and (49) can also be represented as

$$\begin{array}{}\text{(50)}& \begin{array}{rl}& \left(u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right),v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)\right)=\sum _{m=\mathrm{0}}^M\sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\\ & \left({\mathit{\chi }}_{n,m}^{\mathrm{c}}{\mathit{V}}_{n,m}^{\mathrm{1}}+{\mathit{\chi }}_{n,m}^{\mathrm{s}}{\mathit{V}}_{n,m}^{\mathrm{2}}+{\mathit{\psi }}_{n,m}^{\mathrm{c}}{\mathit{V}}_{n,m}^{\mathrm{3}}+{\mathit{\psi }}_{n,m}^{\mathrm{s}}{\mathit{V}}_{n,m}^{\mathrm{4}}\right),\end{array}\end{array}$$

where we define the new DFS vector basis functions ${\mathit{V}}_{n,m}^{\mathrm{1}}$, ${\mathit{V}}_{n,m}^{\mathrm{2}}$, ${\mathit{V}}_{n,m}^{\mathrm{3}}$, and ${\mathit{V}}_{n,m}^{\mathrm{4}}$ as

$$\begin{array}{}\text{(51a)}& {\displaystyle }{\mathit{V}}_{n,m}^{\mathrm{1}}\left(\mathit{\lambda },\mathit{\theta }\right)\equiv \left(-{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\sin m\mathit{\lambda },-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\cos m\mathit{\lambda }\right),\text{(51b)}& {\displaystyle }{\mathit{V}}_{n,m}^{\text{2}}\left(\mathit{\lambda },\mathit{\theta }\right)\equiv \left({\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\cos m\mathit{\lambda },-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\sin m\mathit{\lambda }\right),\text{(51c)}& {\displaystyle }{\mathit{V}}_{n,m}^{\text{3}}\left(\mathit{\lambda },\mathit{\theta }\right)\equiv \left({\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\cos m\mathit{\lambda },-{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\sin m\mathit{\lambda }\right),\text{(51d)}& {\displaystyle }{\mathit{V}}_{n,m}^{\text{4}}\left(\mathit{\lambda },\mathit{\theta }\right)\equiv \left({\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\sin m\mathit{\lambda },{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\cos m\mathit{\lambda }\right).\end{array}$$

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

$$\begin{array}{}\text{(52)}& \begin{array}{rl}& \left[\begin{array}{c}{u}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\\ v_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\end{array}\right]\\ & =\left\{\begin{array}{l}\sum _{n=\mathrm{1}}^N\left[\begin{array}{c}{u}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ v_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\end{array}\right]\sin n\mathit{\theta }\text{for\hspace{0.17em} }m=\mathrm{0},\\ \sum _{n=\mathrm{0}}^N\left[\begin{array}{c}{u}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ v_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\end{array}\right]\cos n\mathit{\theta }\text{for }m=\mathrm{1},\\ \sum _{n=\mathrm{1}}^N\left[\begin{array}{c}{u}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ v_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\end{array}\right]\sin n\mathit{\theta }\text{for\hspace{0.17em} even }m\ge \mathrm{2},\\ \sum _{n=\mathrm{0}}^N\left[\begin{array}{c}{u}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ v_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\end{array}\right]\cos n\mathit{\theta }\left(=\sum _{n=\mathrm{1}}^{N-\mathrm{1}}\left[\begin{array}{c}{u}_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}\\ v_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}\end{array}\right]\sin \mathit{\theta }\sin n\mathit{\theta }\right)\\ \mathrm{\text{ for odd }}m\ge \mathrm{3},\end{array}\right.\end{array}\end{array}$$

where

$$\begin{array}{}\text{(53a)}& {\displaystyle }{u}_{n,m}^{\mathrm{c}}={\displaystyle \frac{\mathrm{1}}{a}}\left[-n{\mathit{\psi }}_{n,m}^{\mathrm{c}}\right]\left(n=\mathrm{1},\mathrm{\dots },N\right)\text{for\hspace{0.17em} }m=\mathrm{0},\text{(53b)}& {\displaystyle }\begin{array}{rl}{u}_{n,m}^{\mathrm{c}}& ={\displaystyle \frac{\mathrm{1}}{a}}\left[m{\mathit{\chi }}_{n,m}^{\mathrm{s}}+{\displaystyle \frac{n\left({\mathit{\psi }}_{n-\mathrm{1},m}^{\mathrm{c}}-{\mathit{\psi }}_{n+\mathrm{1},m}^{\mathrm{c}}\right)}{\mathrm{2}}}\right]\\ & \left(n=\mathrm{0},\mathrm{\dots },N\right)\text{for }m=\mathrm{1},\\ & \text{except\hspace{0.17em} for\hspace{0.17em} }{u}_{\mathrm{1},m}^{\mathrm{c}}={\displaystyle \frac{\mathrm{1}}{a}}\left[m{\mathit{\chi }}_{\mathrm{1},m}^{\mathrm{s}}+{\displaystyle \frac{\mathrm{2}{\mathit{\psi }}_{\mathrm{0},m}^{\mathrm{c}}-{\mathit{\psi }}_{\mathrm{2},m}^{\mathrm{c}}}{\mathrm{2}}}\right]\\ & \left(n=\mathrm{1}\right),\end{array}\text{(53c)}& {\displaystyle }\begin{array}{rl}{u}_{n,m}^{\mathrm{c}}& ={\displaystyle \frac{\mathrm{1}}{a}}\left[m{\mathit{\chi }}_{n,m}^{\mathrm{s}}+{\displaystyle \frac{n\left({\mathit{\psi }}_{n-\mathrm{1},m}^{\mathrm{c}}-{\mathit{\psi }}_{n+\mathrm{1},m}^{\mathrm{c}}\right)}{\mathrm{2}}}\right]\\ & \left(n=\mathrm{1},\mathrm{\dots },N\right)\text{for\hspace{0.17em} even\hspace{0.17em} }m\ge \mathrm{2},\end{array}\text{(53d)}& {\displaystyle }\begin{array}{rl}{u}_{n,m}^{\mathrm{c}}& ={\displaystyle \frac{\mathrm{1}}{a}}\left[{\displaystyle \frac{m\left(-{\mathit{\chi }}_{n-\mathrm{1},m}^{\mathrm{s}}+{\mathit{\chi }}_{n+\mathrm{1},m}^{\mathrm{s}}\right)}{\mathrm{2}}}\right.\\ & \left.+{\displaystyle \frac{n\left(-{\mathit{\psi }}_{n-\mathrm{2},m}^{\mathrm{c}}+\mathrm{2}{\mathit{\psi }}_{n,m}^{\mathrm{c}}-{\mathit{\psi }}_{n+\mathrm{2},m}^{\mathrm{c}}\right)}{\mathrm{4}}}\right]\\ & \left(n=\mathrm{0},\mathrm{\dots },N\right)\text{for\hspace{0.17em} odd\hspace{0.17em} }m\ge \mathrm{3}\\ & \text{except\hspace{0.17em} for\hspace{0.17em} }{u}_{\mathrm{1},m}^{\mathrm{c}}={\displaystyle \frac{\mathrm{1}}{a}}\left[{\displaystyle \frac{m{\mathit{\chi }}_{\mathrm{2},m}^{\mathrm{s}}}{\mathrm{2}}}\right.\\ & \left.+{\displaystyle \frac{\left(\mathrm{3}{\mathit{\psi }}_{\mathrm{1},m}^{\mathrm{c}}-{\mathit{\psi }}_{\mathrm{3},m}^{\mathrm{c}}\right)}{\mathrm{4}}}\right]\left(n=\mathrm{1}\right).\end{array}\end{array}$$

The equations for $u_{n,m}^{\mathrm{s}}$ are the same as Eqs. (53b)–(53d), except that $u_{n,m}^{\mathrm{c}}$, ${\mathit{\chi }}_{n,m}^{\mathrm{s}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{c}}$ are replaced with $u_{n,m}^{\mathrm{s}}$, $-{\mathit{\chi }}_{n,m}^{\mathrm{c}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$, respectively. The equations for $v_{n,m}^{\mathrm{c}}$ are the same as Eqs. (53a)–(53d), except that $u_{n,m}^{\mathrm{c}}$, ${\mathit{\chi }}_{n,m}^{\mathrm{s}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{c}}$ are replaced with $v_{n,m}^{\mathrm{c}}$, ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$, and $-{\mathit{\chi }}_{n,m}^{\mathrm{c}}$, respectively. The equations for $v_{n,m}^{\mathrm{s}}$ are the same as Eqs. (53b)–(53d), except that $u_{n,m}^{\mathrm{c}}$, ${\mathit{\chi }}_{n,m}^{\mathrm{s}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{c}}$ are replaced with $v_{n,m}^{\mathrm{s}}$, $-{\mathit{\psi }}_{n,m}^{\mathrm{c}}$, and $-{\mathit{\chi }}_{n,m}^{\mathrm{s}}$, respectively.

From Eqs. (52) and (53), it can be seen that $u_m^{\mathrm{c},N}\left(\mathit{\theta }\right)$, $u_m^{\mathrm{s},N}\left(\mathit{\theta }\right)$, $v_m^{\mathrm{c},N}\left(\mathit{\theta }\right)$, and $v_m^{\mathrm{s},N}\left(\mathit{\theta }\right)$ at the poles are finite for m=1 and zero for m≠1. Moreover, the following relations are satisfied for m=1:

$$\begin{array}{}\text{(54a)}& {\displaystyle }\begin{array}{rl}{u}_{m=\mathrm{1}}^{\mathrm{c},N}\left(\mathit{\theta }\right)& =-v_{m=\mathrm{1}}^{\mathrm{s},N}\left(\mathit{\theta }\right)\left(={\displaystyle \frac{\mathrm{1}}{a}}\sum _{n=\mathrm{0}}^{N-\mathrm{1}}\left({\mathit{\chi }}_{n,m=\mathrm{1}}^{\mathrm{s}}+{\mathit{\psi }}_{n,m=\mathrm{1}}^{\mathrm{c}}\right)\right)\\ & \text{at\hspace{0.17em} }\mathit{\theta }=\mathrm{0}\left(\text{North\hspace{0.17em} Pole}\right),\end{array}\text{(54b)}& {\displaystyle }\begin{array}{rl}{u}_{m=\mathrm{1}}^{s,N}\left(\mathit{\theta }\right)& =v_{m=\mathrm{1}}^{\text{c,}N}\left(\mathit{\theta }\right)\left(={\displaystyle \frac{\mathrm{1}}{a}}\sum _{n=\mathrm{0}}^{N-\mathrm{1}}\left(-{\mathit{\chi }}_{n,m=\mathrm{1}}^{\mathrm{c}}+{\mathit{\psi }}_{n,m=\mathrm{1}}^{\mathrm{s}}\right)\right)\\ & \text{at\hspace{0.17em} }\mathit{\theta }=\mathrm{0}\left(\text{North\hspace{0.17em} Pole}\right),\end{array}\text{(54c)}& {\displaystyle }\begin{array}{rl}{u}_{m=\mathrm{1}}^{\mathrm{c},N}\left(\mathit{\theta }\right)& =v_{m=\mathrm{1}}^{\mathrm{s},N}\left(\mathit{\theta }\right)\left(={\displaystyle \frac{\mathrm{1}}{a}}\sum _{n=\mathrm{0}}^{N-\mathrm{1}}{\left(-\mathrm{1}\right)}^n\left({\mathit{\chi }}_{n,m=\mathrm{1}}^{\mathrm{s}}-{\mathit{\psi }}_{n,m=\mathrm{1}}^{\mathrm{c}}\right)\right)\\ & \text{at\hspace{0.17em} }\mathit{\theta }=\mathit{\pi }\left(\text{South\hspace{0.17em} Pole}\right),\end{array}\text{(54d)}& {\displaystyle }\begin{array}{rl}{u}_{m=\mathrm{1}}^{\mathrm{s},N}\left(\mathit{\theta }\right)& =-v_{m=\mathrm{1}}^{\mathrm{c},N}\left(\mathit{\theta }\right)\left(={\displaystyle \frac{\mathrm{1}}{a}}\sum _{n=\mathrm{0}}^{N-\mathrm{1}}{\left(-\mathrm{1}\right)}^n\left(-{\mathit{\chi }}_{n,m=\mathrm{1}}^{\mathrm{c}}-{\mathit{\psi }}_{n,m=\mathrm{1}}^{\mathrm{s}}\right)\right)\\ & \text{at\hspace{0.17em} }\mathit{\theta }=\mathit{\pi }\left(\text{South\hspace{0.17em} Pole}\right).\end{array}\end{array}$$

Thus, it is guaranteed that the wind vector $\left(u^{N,M},v^{N,M}\right)$ 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 ${\mathit{\chi }}_{n,m}^{\mathrm{c}}$, ${\mathit{\chi }}_{n,m}^{\mathrm{s}}$, ${\mathit{\psi }}_{n,m}^{\mathrm{c}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$ in Eqs. (48)–(50) using the least-squares method to minimize the error of $u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ and $v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ from u(λ,θ) and v(λ,θ) due to the meridional wavenumber truncation. First, the wind vector components u and v are expanded zonally as

$$\begin{array}{}\text{(55)}& \begin{array}{rl}\left[\begin{array}{c}u\left(\mathit{\lambda },\mathit{\theta }\right)\\ v\left(\mathit{\lambda },\mathit{\theta }\right)\end{array}\right]& \cong \sum _{m=\mathrm{0}}^M\left[\begin{array}{c}{u}_m^{\mathrm{c}}\left(\mathit{\theta }\right)\\ v_m^{\mathrm{c}}\left(\mathit{\theta }\right)\end{array}\right]\cos m\mathit{\lambda }\\ & +\sum _{m=\mathrm{1}}^M\left[\begin{array}{c}{u}_m^{\mathrm{s}}\left(\mathit{\theta }\right)\\ v_m^{\mathrm{s}}\left(\mathit{\theta }\right)\end{array}\right]\sin m\mathit{\lambda },\end{array}\end{array}$$

where $u_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ and $v_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ are calculated from u(λ,θ) and v(λ,θ) by the forward Fourier transform as

$$\begin{array}{}\text{(56a)}& {\displaystyle }\begin{array}{rl}\left[\begin{array}{c}{u}_m^{\mathrm{c}}\left(\mathit{\theta }\right)\\ v_m^{\mathrm{c}}\left(\mathit{\theta }\right)\end{array}\right]& ={\displaystyle \frac{a}{\mathrm{2}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\cos m\mathit{\lambda }\left[\begin{array}{c}u\left(\mathit{\lambda },\mathit{\theta }\right)\\ v\left(\mathit{\lambda },\mathit{\theta }\right)\end{array}\right]\mathrm{d}\mathit{\lambda },\\ & a\equiv \left\{\begin{array}{l}\mathrm{1}\text{ for}m=\mathrm{0}\\ \mathrm{2}\text{ for}m\ge \mathrm{1},\end{array}\right.\end{array}\text{(56b)}& {\displaystyle }\left[\begin{array}{c}{u}_m^{\mathrm{s}}\left(\mathit{\theta }\right)\\ v_m^{\mathrm{s}}\left(\mathit{\theta }\right)\end{array}\right]={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\sin m\mathit{\lambda }\left[\begin{array}{c}u\left(\mathit{\lambda },\mathit{\theta }\right)\\ v\left(\mathit{\lambda },\mathit{\theta }\right)\end{array}\right]\mathrm{d}\mathit{\lambda }.\end{array}$$

The variables $u_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ and $v_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ are meridionally expanded as

$$\begin{array}{}\text{(57)}& \begin{array}{rl}& \left[\begin{array}{c}{u}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\\ v_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\end{array}\right]\cong \left[\begin{array}{c}{\stackrel{\mathrm{̃}}{u}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\\ {\stackrel{\mathrm{̃}}{v}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\end{array}\right]\\ & \equiv \left\{\begin{array}{l}\sum _{n=\mathrm{1}}^N\left[\begin{array}{c}{\stackrel{\mathrm{̃}}{u}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\end{array}\right]\sin n\mathit{\theta }\text{for\hspace{0.17em} even\hspace{0.17em} }m,\\ \sum _{n=\mathrm{0}}^N\left[\begin{array}{c}{\stackrel{\mathrm{̃}}{u}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\end{array}\right]\cos n\mathit{\theta }\text{for\hspace{0.17em} odd}m,\end{array}\right.\end{array}\end{array}$$

where ${\stackrel{\mathrm{̃}}{u}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and ${\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ are calculated from $u_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ and $v_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ by the forward discrete cosine or sine transform (see Appendix B).

Next, ${\mathit{\chi }}_{n,m}^{\mathrm{c}}$, ${\mathit{\chi }}_{n,m}^{\mathrm{s}}$, ${\mathit{\psi }}_{n,m}^{\mathrm{c}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$ are calculated to minimize the following error F (the squared L₂ norm of the residual vector):

$$\begin{array}{}\text{(58)}& F\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\left(R^u{\left(\mathit{\lambda },\mathit{\theta }\right)}^{\mathrm{2}}+R^v{\left(\mathit{\lambda },\mathit{\theta }\right)}^{\mathrm{2}}\right)\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda },\end{array}$$

where the residual vector $\left(R^u\left(\mathit{\lambda },\mathit{\theta }\right),R^v\left(\mathit{\lambda },\mathit{\theta }\right)\right)$ is defined as

$$\begin{array}{}\text{(59a)}& {\displaystyle }{R}^u\left(\mathit{\lambda },\mathit{\theta }\right)\equiv u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)-u\left(\mathit{\lambda },\mathit{\theta }\right),\text{(59b)}& {\displaystyle }{R}^v\left(\mathit{\lambda },\mathit{\theta }\right)\equiv v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)-v\left(\mathit{\lambda },\mathit{\theta }\right).\end{array}$$

From Eqs. (58) and (59) and the equations $\partial F/\partial {\mathit{\chi }}_{m,n}^{\mathrm{c}}=\mathrm{0}$, $\partial F/\partial {\mathit{\chi }}_{n,m}^{\mathrm{s}}=\mathrm{0}$, $\partial F/\partial {\mathit{\psi }}_{n,m}^{\mathrm{c}}=\mathrm{0}$, and $\partial F/\partial {\mathit{\psi }}_{n,m}^{\mathrm{s}}=\mathrm{0}$ used in the least-squares method, we obtain

$$\begin{array}{}\text{(60a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \left[{\displaystyle \frac{\partial u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\chi }}_{n,m}^{\mathrm{c}}}}{R}^u\left(\mathit{\lambda },\mathit{\theta }\right)\right.\\ & \left.+{\displaystyle \frac{\partial v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\chi }}_{n,m}^{\mathrm{c}}}}{R}^v\left(\mathit{\lambda },\mathit{\theta }\right)\right]\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }=\mathrm{0},\end{array}\text{(60b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \left[{\displaystyle \frac{\partial u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\chi }}_{n,m}^{\mathrm{s}}}}{R}^u\left(\mathit{\lambda },\mathit{\theta }\right)\right.\\ & \left.+{\displaystyle \frac{\partial v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\chi }}_{n,m}^{\mathrm{s}}}}{R}^v\left(\mathit{\lambda },\mathit{\theta }\right)\right]\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }=\mathrm{0},\end{array}\text{(60c)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \left[{\displaystyle \frac{\partial u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\psi }}_{n,m}^{\mathrm{c}}}}{R}^u\left(\mathit{\lambda },\mathit{\theta }\right)\right.\\ & \left.+{\displaystyle \frac{\partial v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\psi }}_{n,m}^{\mathrm{c}}}}{R}^v\left(\mathit{\lambda },\mathit{\theta }\right)\right]\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }=\mathrm{0},\end{array}\text{(60d)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \left[{\displaystyle \frac{\partial u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\psi }}_{n,m}^{\mathrm{s}}}}{R}^u\left(\mathit{\lambda },\mathit{\theta }\right)\right.\\ & \left.+{\displaystyle \frac{\partial v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\psi }}_{n,m}^{\mathrm{s}}}}{R}^v\left(\mathit{\lambda },\mathit{\theta }\right)\right]\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }=\mathrm{0}.\end{array}\end{array}$$

From Eq. (50), we derive

$$\begin{array}{}\text{(61a)}& {\displaystyle }\left({\displaystyle \frac{\partial u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\chi }}_{n,m}^{\mathrm{c}}}},{\displaystyle \frac{\partial v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\chi }}_{n,m}^{\mathrm{c}}}}\right)={\mathit{V}}_{n,m}^{\mathrm{1}}\left(\mathit{\lambda },\mathit{\theta }\right),\text{(61b)}& {\displaystyle }\left({\displaystyle \frac{\partial u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\chi }}_{n,m}^{\mathrm{s}}}},{\displaystyle \frac{\partial v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\chi }}_{n,m}^{\mathrm{s}}}}\right)={\mathit{V}}_{n,m}^{\mathrm{2}}\left(\mathit{\lambda },\mathit{\theta }\right),\text{(61c)}& {\displaystyle }\left({\displaystyle \frac{\partial u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\psi }}_{n,m}^{\mathrm{c}}}},{\displaystyle \frac{\partial v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\psi }}_{n,m}^{\mathrm{c}}}}\right)={\mathit{V}}_{n,m}^{\mathrm{3}}\left(\mathit{\lambda },\mathit{\theta }\right),\text{(61d)}& {\displaystyle }\left({\displaystyle \frac{\partial u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\psi }}_{n,m}^{\mathrm{s}}}},{\displaystyle \frac{\partial v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)}{\partial {\mathit{\psi }}_{n,m}^{\mathrm{s}}}}\right)={\mathit{V}}_{n,m}^{\mathrm{4}}\left(\mathit{\lambda },\mathit{\theta }\right).\end{array}$$

Equations (60) and (61) show that the residual vector $\left(R^u\left(\mathit{\lambda },\mathit{\theta }\right),R^v\left(\mathit{\lambda },\mathit{\theta }\right)\right)$ 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

$$\begin{array}{}\text{(62a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \left[-{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\left(u_m^{\mathrm{s},N}\left(\mathit{\theta }\right)-{\stackrel{\mathrm{̃}}{u}}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\right)\right.\\ & \left.-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\left(v_m^{\mathrm{c},N}\left(\mathit{\theta }\right)-{\stackrel{\mathrm{̃}}{v}}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\right)\right]\mathrm{d}\mathit{\theta }=\mathrm{0},\end{array}\text{(62b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \left[{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\left(u_m^{\mathrm{c},N}\left(\mathit{\theta }\right)-{\stackrel{\mathrm{̃}}{u}}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\right)\right.\\ & \left.-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\left(v_m^{\mathrm{s},N}\left(\mathit{\theta }\right)-{\stackrel{\mathrm{̃}}{v}}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\right)\right]\mathrm{d}\mathit{\theta }=\mathrm{0},\end{array}\text{(62c)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \left[{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\left(u_m^{\mathrm{c},N}\left(\mathit{\theta }\right)-{\stackrel{\mathrm{̃}}{u}}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\right)\right.\\ & \left.-{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\left(v_m^{\mathrm{s},N}\left(\mathit{\theta }\right)-{\stackrel{\mathrm{̃}}{v}}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\right)\right]\mathrm{d}\mathit{\theta }=\mathrm{0},\end{array}\text{(62d)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \left[{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\left(u_m^{\mathrm{s},N}\left(\mathit{\theta }\right)-{\stackrel{\mathrm{̃}}{u}}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\right)\right.\\ & \left.+{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{a\sin \mathit{\theta }}}\left(v_m^{\mathrm{c},N}\left(\mathit{\theta }\right)-{\stackrel{\mathrm{̃}}{v}}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\right)\right]\mathrm{d}\mathit{\theta }=\mathrm{0}.\end{array}\end{array}$$

By substituting Eqs. (52) and (57) into Eq. (62a), the following equations for ${\mathit{\chi }}_{n,m}^{\mathrm{c}}$ and ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$ are derived as shown in Appendix H.

For m=0,

$$\begin{array}{}\text{(63a)}& {\displaystyle \frac{\mathrm{1}}{a}}\left[n{\mathit{\chi }}_{n,m}^{\mathrm{c}}\right]={\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{c}}\left(\mathrm{1}\le n\le N\right).\end{array}$$

The coefficient ${\mathit{\chi }}_{m=\mathrm{0},n=\mathrm{0}}^{\mathrm{c}}$ 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,

$$\begin{array}{}\text{(63b)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{a}}& \left[-{\left(n-\mathrm{1}\right)}^{\mathrm{2}}{\mathit{\chi }}_{n-\mathrm{2},m}^{\mathrm{c}}-\mathrm{2}m{\mathit{\psi }}_{n-\mathrm{1},m}^{\mathrm{s}}\right.\\ & +\left(\mathrm{4}{m}^{\mathrm{2}}+\mathrm{2}{n}^{\mathrm{2}}+\mathrm{2}\right){\mathit{\chi }}_{n,m}^{\mathrm{c}}\\ & \left.-\mathrm{2}m{\mathit{\psi }}_{n+\mathrm{1},m}^{\mathrm{s}}-{\left(n+\mathrm{1}\right)}^{\mathrm{2}}{\mathit{\chi }}_{n+\mathrm{2},m}^{\mathrm{c}}\right]\\ & =\mathrm{2}\left(n-\mathrm{1}\right){\stackrel{\mathrm{̃}}{v}}_{n-\mathrm{1},m}^{\mathrm{c}}-\mathrm{4}m{\stackrel{\mathrm{̃}}{u}}_{n,m}^{\mathrm{s}}-\mathrm{2}\left(n+\mathrm{1}\right){\stackrel{\mathrm{̃}}{v}}_{n+\mathrm{1},m}^{\mathrm{c}}\\ & \left(\mathrm{0}\le n\le N-\mathrm{1}\right),\end{array}\end{array}$$

with the exception of the following underlined values:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{a}}\left[\left(\underset{\mathrm{¯}}{\mathrm{8}}{m}^{\mathrm{2}}+\underset{\mathrm{¯}}{\mathrm{4}}\right){\mathit{\chi }}_{\mathrm{0},m}^{\mathrm{c}}-\underset{\mathrm{¯}}{\mathrm{4}}m{\mathit{\psi }}_{\mathrm{1},m}^{\mathrm{s}}-\underset{\mathrm{¯}}{\mathrm{2}}{\mathit{\chi }}_{\mathrm{2},m}^{\mathrm{c}}\right]\\ {\displaystyle }& {\displaystyle }=-\underset{\mathrm{¯}}{\mathrm{8}}m{\stackrel{\mathrm{̃}}{u}}_{\mathrm{0},m}^{\mathrm{s}}-\underset{\mathrm{¯}}{\mathrm{4}}{\stackrel{\mathrm{̃}}{v}}_{\mathrm{1},m}^{\mathrm{c}}\left(n=\mathrm{0}\right),\\ {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{a}}\left[-\underset{\mathrm{¯}}{\mathrm{4}}m{\mathit{\psi }}_{\mathrm{0},m}^{\mathrm{s}}+\left(\mathrm{4}{m}^{\mathrm{2}}+\mathrm{4}\right){\mathit{\chi }}_{\mathrm{1},m}^{\mathrm{c}}+\mathrm{\dots }\right]=\mathrm{\dots }\left(n=\mathrm{1}\right),\\ {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{a}}\left[-\underset{\mathrm{¯}}{\mathrm{2}}{\mathit{\chi }}_{\mathrm{0},m}^{\mathrm{c}}-\mathrm{2}m{\mathit{\psi }}_{\mathrm{1},m}^{\mathrm{s}}+\mathrm{\dots }\right]=\mathrm{\dots }\left(n=\mathrm{2}\right).\end{array}$$

For even m≥2,

$$\begin{array}{}\text{(63c)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{a}}& \left[-{\left(n-\mathrm{1}\right)}^{\mathrm{2}}{\mathit{\chi }}_{n-\mathrm{2},m}^{\mathrm{c}}-\mathrm{2}m{\mathit{\psi }}_{n-\mathrm{1},m}^{\mathrm{s}}\right.\\ & +\left(\mathrm{4}{m}^{\mathrm{2}}+\mathrm{2}{n}^{\mathrm{2}}+\mathrm{2}\right){\mathit{\chi }}_{n,m}^{\mathrm{c}}-\mathrm{2}m{\mathit{\psi }}_{n+\mathrm{1},m}^{\mathrm{s}}\\ & \left.-{\left(n+\mathrm{1}\right)}^{\mathrm{2}}{\mathit{\chi }}_{n+\mathrm{2},m}^{\mathrm{c}}\right]\\ & =\mathrm{2}\left(n-\mathrm{1}\right){\stackrel{\mathrm{̃}}{v}}_{n-\mathrm{1},m}^{\mathrm{c}}-\mathrm{4}m{\stackrel{\mathrm{̃}}{u}}_{n,m}^{\mathrm{s}}-\mathrm{2}\left(n+\mathrm{1}\right){\stackrel{\mathrm{̃}}{v}}_{n+\mathrm{1},m}^{\mathrm{c}}\\ & \left(\mathrm{1}\le n\le N-\mathrm{1}\right),\end{array}\end{array}$$

with no exception.

For odd m≥3,

$$\begin{array}{}\text{(63d)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{a}}& \left[{\left(n-\mathrm{2}\right)}^{\mathrm{2}}{\mathit{\chi }}_{n-\mathrm{4},m}^{\mathrm{c}}+\mathrm{2}m{\mathit{\psi }}_{n-\mathrm{3},m}^{\mathrm{s}}\right.\\ & +\left(-\mathrm{4}{m}^{\mathrm{2}}-\mathrm{4}{n}^{\mathrm{2}}+\mathrm{8}n-\mathrm{8}\right){\mathit{\chi }}_{n-\mathrm{2},m}^{\mathrm{c}}\\ & -\mathrm{2}m{\mathit{\psi }}_{n-\mathrm{1},m}^{\mathrm{s}}+\left(\mathrm{8}{m}^{\mathrm{2}}+\mathrm{6}{n}^{\mathrm{2}}+\mathrm{8}\right){\mathit{\chi }}_{n,m}^{\mathrm{c}}\\ & -\mathrm{2}m{\mathit{\psi }}_{n+\mathrm{1},m}^{\mathrm{s}}+\left(-\mathrm{4}{m}^{\mathrm{2}}-\mathrm{4}{n}^{\mathrm{2}}-\mathrm{8}n-\mathrm{8}\right){\mathit{\chi }}_{n+\mathrm{2},m}^{\mathrm{c}}\\ & \left.+\mathrm{2}m{\mathit{\psi }}_{n+\mathrm{3},m}^{\mathrm{s}}+{\left(n+\mathrm{2}\right)}^{\mathrm{2}}{\mathit{\chi }}_{n+\mathrm{4},m}^{\mathrm{c}}\right]\\ & =\mathrm{4}\left(n-\mathrm{2}\right){\stackrel{\mathrm{̃}}{v}}_{n-\mathrm{2},m}^{\mathrm{c}}-\mathrm{8}m{\stackrel{\mathrm{̃}}{u}}_{n-\mathrm{1},m}^{\mathrm{s}}-\mathrm{8}n{\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{c}}\\ & +\mathrm{8}m{\stackrel{\mathrm{̃}}{u}}_{n+\mathrm{1},m}^{\mathrm{s}}+\mathrm{4}\left(n+\mathrm{2}\right){\stackrel{\mathrm{̃}}{v}}_{n+\mathrm{2},m}^{\mathrm{c}}\left(\mathrm{1}\le n\le N-\mathrm{2}\right),\end{array}\end{array}$$

with the exception of the following underlined values:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{a}}\left[\left(\underset{\mathrm{¯}}{\mathrm{12}}{m}^{\mathrm{2}}+\underset{\mathrm{¯}}{\mathrm{18}}\right){\mathit{\chi }}_{\mathrm{1},m}^{\mathrm{c}}-\underset{\mathrm{¯}}{\mathrm{4}}m{\mathit{\psi }}_{\mathrm{2},m}^{\mathrm{s}}+\left(-\mathrm{4}{m}^{\mathrm{2}}-\underset{\mathrm{¯}}{\mathrm{21}}\right){\mathit{\chi }}_{\mathrm{3},m}^{\mathrm{c}}+\mathrm{\dots }\right]\\ {\displaystyle }& {\displaystyle }=-\underset{\mathrm{¯}}{\mathrm{16}}m{\stackrel{\mathrm{̃}}{u}}_{\mathrm{0},m}^{\mathrm{s}}-\underset{\mathrm{¯}}{\mathrm{12}}{\stackrel{\mathrm{̃}}{v}}_{\mathrm{1},m}^{\mathrm{c}}+\mathrm{\dots }\left(n=\mathrm{1}\right),\\ {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{a}}\left[-\underset{\mathrm{¯}}{\mathrm{4}}m{\mathit{\psi }}_{\mathrm{1},m}^{\mathrm{s}}+\left(\mathrm{8}{m}^{\mathrm{2}}+\mathrm{32}\right){\mathit{\chi }}_{\mathrm{2},m}^{\mathrm{c}}+\mathrm{\dots }\right]=\mathrm{\dots }\left(n=\mathrm{2}\right),\\ {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{a}}\left[\left(-\mathrm{4}{m}^{\mathrm{2}}-\underset{\mathrm{¯}}{\mathrm{21}}\right){\mathit{\chi }}_{\mathrm{1},m}^{\mathrm{c}}-\mathrm{2}m{\mathit{\psi }}_{\mathrm{2},m}^{\mathrm{s}}+\mathrm{\dots }\right]=\mathrm{\dots }\left(n=\mathrm{3}\right).\end{array}$$

Similarly, from Eq. (62b) we derive the same equations as Eqs. (63b)–(63d), except that χ^c, ψ^s, ${\stackrel{\mathrm{̃}}{v}}^{\mathrm{c}}$, and ${\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}$ are replaced with χ^s, −ψ^c, ${\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}$, and $-{\stackrel{\mathrm{̃}}{u}}^{\mathrm{c}}$, respectively. From Eq. (62c), we derive the same equations as Eqs. (63a)–(63d), except that χ^c, ψ^s, ${\stackrel{\mathrm{̃}}{v}}^{\mathrm{c}}$, and ${\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}$ are replaced with −ψ^c, χ^s, ${\stackrel{\mathrm{̃}}{u}}^{\mathrm{c}}$, and $-{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}$, respectively. From Eq. (62d), we derive the same equations as Eqs. (63b)–(63d), except that χ^c, ψ^s, ${\stackrel{\mathrm{̃}}{v}}^{\mathrm{c}}$, and ${\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}$ are replaced with ψ^s, χ^c, $-{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}$, and $-{\stackrel{\mathrm{̃}}{v}}^{\mathrm{c}}$, respectively.

Equation (63a) is easily solved. From Eq. (63d) and from the same equations as Eq. (63d), except that χ^c, ψ^s, ${\stackrel{\mathrm{̃}}{v}}^{\mathrm{c}}$, and ${\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}$ are replaced with ψ^s, χ^c, $-{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}$, and $-{\stackrel{\mathrm{̃}}{v}}^{\mathrm{c}}$, respectively, we derive the following linear simultaneous equations for m≥3:

$$\begin{array}{}\text{(64)}& \begin{array}{rl}& {\mathbf{E}}_m\left[\begin{array}{l}{\mathit{\chi }}_{\mathrm{1},m}^{\mathrm{c}}\\ {\mathit{\psi }}_{\mathrm{2},m}^{\mathrm{s}}\\ {\mathit{\chi }}_{\mathrm{3},m}^{\mathrm{c}}\\ {\mathit{\psi }}_{\mathrm{4},m}^{\mathrm{s}}\\ :\end{array}\right]={\mathbf{F}}_m\left[\begin{array}{l}{\stackrel{\mathrm{̃}}{u}}_{\mathrm{0},m}^{\mathrm{s}}\\ {\stackrel{\mathrm{̃}}{v}}_{\mathrm{1},m}^{\mathrm{c}}\\ {\stackrel{\mathrm{̃}}{u}}_{\mathrm{2},m}^{\mathrm{s}}\\ {\stackrel{\mathrm{̃}}{v}}_{\mathrm{3},m}^{\mathrm{c}}\\ :\end{array}\right],\\ & {\mathbf{E}}_m\left[\begin{array}{l}{\mathit{\psi }}_{\mathrm{1},m}^{\mathrm{s}}\\ {\mathit{\chi }}_{\mathrm{2},m}^{\mathrm{c}}\\ {\mathit{\psi }}_{\mathrm{3},m}^{\mathrm{s}}\\ {\mathit{\chi }}_{\mathrm{4},m}^{\mathrm{c}}\\ :\end{array}\right]={\mathbf{F}}_m\left[\begin{array}{l}-{\stackrel{\mathrm{̃}}{v}}_{\mathrm{0},m}^{\mathrm{c}}\\ -{\stackrel{\mathrm{̃}}{u}}_{\mathrm{1},m}^{\mathrm{s}}\\ -{\stackrel{\mathrm{̃}}{v}}_{\mathrm{2},m}^{\mathrm{c}}\\ -{\stackrel{\mathrm{̃}}{u}}_{\mathrm{3},m}^{\mathrm{s}}\\ :\end{array}\right],\end{array}\end{array}$$

where the matrices E_m 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 E_m. The simultaneous equations with nona-diagonal or penta-diagonal matrices E_m can be solved in a similar way to Eq. (31), and the expansion coefficients ${\mathit{\chi }}_{n,m}^{\mathrm{c}}$ and ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$ in Eq. (64) can be calculated efficiently. From the same equations as Eqs. (63b)–(63d), except that χ^c, ψ^s, ${\stackrel{\mathrm{̃}}{v}}^{\mathrm{c}}$, and ${\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}$ are replaced with χ^s, −ψ^c, ${\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}$, and $-{\stackrel{\mathrm{̃}}{u}}^{\mathrm{c}}$, respectively, and the same equations as Eqs. (63b)–(63d), except that χ^c, ψ^s, ${\stackrel{\mathrm{̃}}{v}}^{\mathrm{c}}$, and ${\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}$ are replaced with −ψ^c, χ^s, ${\stackrel{\mathrm{̃}}{u}}^{\mathrm{c}}$, and $-{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}$, respectively, the simultaneous equations similar to Eq. (64) are also derived. Thus, the expansion coefficients ${\mathit{\chi }}_{n,m}^{\mathrm{c}}$, ${\mathit{\chi }}_{n,m}^{\mathrm{s}}$, ${\mathit{\psi }}_{n,m}^{\mathrm{c}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$ are calculated from ${\stackrel{\mathrm{̃}}{u}}_{n,m}^{\mathrm{c}}$, ${\stackrel{\mathrm{̃}}{u}}_{n,m}^{\mathrm{s}}$, ${\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{c}}$, and ${\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{s}}$ using Eqs. (63a)–(63d) and the similar equations. The expansion coefficients $u_{n,m}^{\mathrm{c}}$, $u_{n,m}^{\mathrm{s}}$, $v_{n,m}^{\mathrm{c}}$, and $v_{n,m}^{\mathrm{s}}$ are calculated from ${\mathit{\chi }}_{n,m}^{\mathrm{c}}$, ${\mathit{\chi }}_{n,m}^{\mathrm{s}}$, ${\mathit{\psi }}_{n,m}^{\mathrm{c}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$ using Eq. (53) for $u_{n,m}^{\mathrm{c}}$ and the similar equations for $u_{n,m}^{\mathrm{s}}$, $v_{n,m}^{\mathrm{c}}$, and $v_{n,m}^{\mathrm{s}}$.

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=∇²χ and the vorticity ζ=∇²ψ 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

$$\begin{array}{}\text{(65)}& \begin{array}{rl}g\left(\mathit{\lambda },\mathit{\theta }\right)& ={\mathrm{\nabla }}^{\mathrm{2}}f\left(\mathit{\lambda },\mathit{\theta }\right)={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[{\displaystyle \frac{\mathrm{1}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}{\displaystyle \frac{{\partial }^{\mathrm{2}}f}{\partial {\mathit{\lambda }}^{\mathrm{2}}}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial f}{\partial \mathit{\theta }}}\right)\right],\end{array}\end{array}$$

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

$$\begin{array}{}\text{(66)}& \begin{array}{rl}\left[\begin{array}{c}f\left(\mathit{\lambda },\mathit{\theta }\right)\\ g\left(\mathit{\lambda },\mathit{\theta }\right)\end{array}\right]& \cong \sum _{m=\mathrm{0}}^M\left[\begin{array}{c}{f}_m^{\mathrm{c}}\left(\mathit{\theta }\right)\\ g_m^{\mathrm{c}}\left(\mathit{\theta }\right)\end{array}\right]\cos m\mathit{\lambda }\\ & +\sum _{m=\mathrm{1}}^M\left[\begin{array}{c}{f}_m^{\mathrm{s}}\left(\mathit{\theta }\right)\\ g_m^{\mathrm{s}}\left(\mathit{\theta }\right)\end{array}\right]\sin m\mathit{\lambda }.\end{array}\end{array}$$

The variables $f_m^{\mathrm{c}}\left(\mathit{\theta }\right)$, $f_m^{\mathrm{s}}\left(\mathit{\theta }\right)$, $g_m^{\mathrm{c}}\left(\mathit{\theta }\right)$, and $g_m^{\mathrm{s}}\left(\mathit{\theta }\right)$ are expanded meridionally like Eq. (10) as

$$\begin{array}{}\text{(67)}& \begin{array}{rl}\left[\begin{array}{c}{f}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\\ g_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\end{array}\right]& \cong \left[\begin{array}{c}{f}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\\ g_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\end{array}\right]\\ & \equiv \sum _{n=N_{\text{min,}m}}^{N_{\text{max,\hspace{0.17em} }m}}\left[\begin{array}{c}{f}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ g_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\end{array}\right]S_{n,m}\left(\mathit{\theta }\right).\end{array}\end{array}$$

We define the truncated variables f^N,M(θ) and g^N,M(θ) as

$$\begin{array}{}\text{(68)}& \begin{array}{rl}\left[\begin{array}{c}{f}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)\\ g^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)\end{array}\right]& \equiv \sum _{m=\mathrm{0}}^M\left[\begin{array}{c}{f}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\\ g_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\end{array}\right]\cos m\mathit{\lambda }\\ & +\sum _{m=\mathrm{1}}^M\left[\begin{array}{c}{f}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\\ g_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\end{array}\right]\sin m\mathit{\lambda }.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(69)}& \begin{array}{rl}{\mathrm{\nabla }}^{\mathrm{2}}{f}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)& =\sum _{m=\mathrm{0}}^M{\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[{\displaystyle \frac{-m^{\mathrm{2}}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}{f}_m^{\text{c,}N}\left(\mathit{\theta }\right)\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial f_m^{\text{c,}N}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\right)\right]\cos m\mathit{\lambda }\\ & +\sum _{m=\mathrm{1}}^M{\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[{\displaystyle \frac{-m^{\mathrm{2}}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}{f}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial f_m^{\text{s,}N}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\right)\right]\sin m\mathit{\lambda }.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(70)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}\left(\mathit{\theta }\right)\left[\begin{array}{c}\cos m\mathit{\lambda }\\ \sin m\mathit{\lambda }\end{array}\right]R^g\left(\mathit{\lambda },\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }=\mathrm{0},\end{array}$$

where the residual is as follows:

$$\begin{array}{}\text{(71)}& R^g\left(\mathit{\lambda },\mathit{\theta }\right)\equiv g^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)-{\mathrm{\nabla }}^{\mathrm{2}}{f}^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right),\end{array}$$

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

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

$$\begin{array}{}\text{(72)}& H\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{R}^g{\left(\mathit{\lambda },\mathit{\theta }\right)}^{\mathrm{2}}\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }.\end{array}$$

When calculating g by applying the Laplacian operator to a given f, $g_{n,m}^{\mathrm{c}}$ and $g_{n,m}^{\mathrm{s}}$ can also be calculated from $\partial H/\partial g_{n,m}^{\mathrm{c}}$ and $\partial H/\partial g_{n,m}^{\mathrm{s}}$ using the least-squares method. The equations $\partial H/\partial g_{n,m}^{\mathrm{c}}$ and $\partial H/\partial g_{n,m}^{\mathrm{s}}$ give the equivalent equations to Eq. (70). When calculating f from a given g in the Poisson equation, $f_{n,m}^{\mathrm{c}}$ and $f_{n,m}^{\mathrm{s}}$ can also be calculated from $\partial H/\partial f_{n,m}^{\mathrm{c}}$ and $\partial H/\partial f_{n,m}^{\mathrm{s}}$ using the least-squares method. However, the equations derived from $\partial H/\partial f_{n,m}^{\mathrm{c}}$ and $\partial H/\partial f_{n,m}^{\mathrm{s}}$ 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

$$\begin{array}{}\text{(73)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}\left(\mathit{\theta }\right)& \left\{g_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)-{\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[{\displaystyle \frac{-m^{\mathrm{2}}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}{f}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\right.\right.\\ & \left.\left.+{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial f_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\right)\right]\right\}\mathrm{d}\mathit{\theta }=\mathrm{0}.\end{array}\end{array}$$

For m=0, we calculate $g_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ by using

$$\begin{array}{}\text{(74)}& \begin{array}{rl}{g}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)& ={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[{\displaystyle \frac{-m^{\mathrm{2}}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}{f}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial f_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}\right)\right],\end{array}\end{array}$$

instead of Eq. (73) following Yee (1981) and Cheong (2000a) for ease of calculation. For $\mathrm{0}\le m\le \mathrm{3}$, the exact solutions of $g_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ can be obtained from Eq. (74) because the new DFS meridional basis functions for $\mathrm{0}\le m\le \mathrm{3}$ are the linear combination of the associated Legendre functions for $\mathrm{0}\le m\le \mathrm{3}$ and vice versa as described in Sect. 3.4.

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

$$\begin{array}{}\text{(75a)}& \begin{array}{rl}-g_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}& +\mathrm{2}{g}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-g_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & ={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[\left(n-\mathrm{1}\right)\left(n-\mathrm{2}\right)f_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right.\\ & \left.-\mathrm{2}{n}^{\mathrm{2}}{f}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}+\left(n+\mathrm{1}\right)\left(n+\mathrm{2}\right)f_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right]\\ & \left(\mathrm{0}\le n\le N\right),\end{array}\end{array}$$

except for the following underlined values:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{¯}}{\mathrm{1}}{g}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-g_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}=\mathrm{\dots }\left(n=\mathrm{1}\right),\\ {\displaystyle }& {\displaystyle }-\underset{\mathrm{¯}}{\mathrm{2}}{g}_{\mathrm{0},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{2}{g}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}-g_{\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}=\mathrm{\dots }\left(n=\mathrm{2}\right).\end{array}$$

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

$$\begin{array}{}\text{(75b)}& \begin{array}{rl}-g_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}& +\mathrm{2}{g}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-g_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & ={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[\left(n-\mathrm{1}\right)n{f}_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right.\\ & \left.-\left(\mathrm{2}{n}^{\mathrm{2}}+\mathrm{4}{m}^{\mathrm{2}}\right)f_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}+\left(n+\mathrm{1}\right)n{f}_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right]\\ & \left(\mathrm{0}\le n\le N-\mathrm{1}\right),\end{array}\end{array}$$

except for the following underlined values:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{¯}}{\mathrm{1}}{g}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-g_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}=\mathrm{\dots }\left(n=\mathrm{1}\right),\\ {\displaystyle }& {\displaystyle }-\underset{\mathrm{¯}}{\mathrm{2}}{g}_{\mathrm{0},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{2}{g}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}-g_{\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[\underset{\mathrm{¯}}{\mathrm{4}}{f}_{\mathrm{0},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{\dots }\right]\left(n=\mathrm{2}\right).\end{array}$$

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

$$\begin{array}{}\text{(75c)}& \begin{array}{rl}-g_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}& +\mathrm{2}{g}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-g_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & ={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[\left(n-\mathrm{1}\right)n{f}_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right.\\ & \left.-\left(\mathrm{2}{n}^{\mathrm{2}}+\mathrm{4}{m}^{\mathrm{2}}\right)f_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}+\left(n+\mathrm{1}\right)n{f}_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right]\\ & \left(\mathrm{1}\le n\le N-\mathrm{1}\right),\end{array}\end{array}$$

except for the following underlined values:

$$\underset{\mathrm{¯}}{\mathrm{3}}{g}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-g_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}=\mathrm{\dots }\left(n=\mathrm{1}\right).$$

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

$$\begin{array}{}\text{(75d)}& \begin{array}{rl}& g_{n-\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{g}_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{6}{g}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{g}_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}+g_{n+\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & ={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[-\left(n-\mathrm{2}\right)\left(n-\mathrm{1}\right)f_{n-\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right.\\ & +\left(\mathrm{4}{n}^{\mathrm{2}}-\mathrm{6}n+\mathrm{4}+\mathrm{4}{m}^{\mathrm{2}}\right)f_{n-\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & -\left(\mathrm{6}{n}^{\mathrm{2}}+\mathrm{4}+\mathrm{8}{m}^{\mathrm{2}}\right)f_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & +\left(\mathrm{4}{n}^{\mathrm{2}}+\mathrm{6}n+\mathrm{4}+\mathrm{4}{m}^{\mathrm{2}}\right)f_{n+\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & \left.-\left(n+\mathrm{2}\right)\left(n+\mathrm{1}\right)f_{n+\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}\right]\left(\mathrm{1}\le n\le N-\mathrm{2}\right),\end{array}\end{array}$$

except for the following underlined values:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{¯}}{\mathrm{10}}{g}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\underset{\mathrm{¯}}{\mathrm{5}}{g}_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}+g_{\mathrm{5},m}^{\mathrm{c}\left(\mathrm{s}\right)}={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[-\left(\underset{\mathrm{¯}}{\mathrm{12}}+\underset{\mathrm{¯}}{\mathrm{12}}{m}^{\mathrm{2}}\right)f_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{\dots }\right]\\ {\displaystyle }& {\displaystyle }\left(n=\mathrm{1}\right),\\ {\displaystyle }& {\displaystyle }\underset{\mathrm{¯}}{\mathrm{5}}{g}_{\mathrm{2},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{g}_{\mathrm{4},m}^{\mathrm{c}\left(\mathrm{s}\right)}+g_{\mathrm{6},m}^{\mathrm{c}\left(\mathrm{s}\right)}=\mathrm{\dots }\left(n=\mathrm{2}\right),\\ {\displaystyle }& {\displaystyle }-\underset{\mathrm{¯}}{\mathrm{5}}{g}_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{6}{g}_{\mathrm{3},m}^{\mathrm{c}\left(\mathrm{s}\right)}-\mathrm{4}{g}_{\mathrm{5},m}^{\mathrm{c}\left(\mathrm{s}\right)}+g_{\mathrm{7},m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[\left(\underset{\mathrm{¯}}{\mathrm{24}}+\mathrm{4}{m}^{\mathrm{2}}\right)f_{\mathrm{1},m}^{\mathrm{c}\left(\mathrm{s}\right)}+\mathrm{\dots }\right]\left(n=\mathrm{3}\right).\end{array}$$

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

$$\begin{array}{}\text{(76)}& \begin{array}{rl}{\mathbf{A}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}{\mathit{g}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{c\left(s\right)}& ={\mathbf{B}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}{\mathit{f}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{\mathrm{c}\left(\mathrm{s}\right)},\\ {\mathbf{A}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}{\mathit{g}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{\mathrm{c}\left(\mathrm{s}\right)}& ={\mathbf{B}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}{\mathit{f}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{\mathrm{c}\left(\mathrm{s}\right)}\end{array}\end{array}$$

where ${\mathit{g}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and ${\mathit{g}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{\mathrm{c}\left(\mathrm{s}\right)}$ are the vectors whose components are $g_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ (n is even) and $g_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ (n is odd), respectively, and ${\mathit{f}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and ${\mathit{f}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{\mathrm{c}\left(\mathrm{s}\right)}$ are the vectors whose components are $f_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ (n is even) and $f_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ (n is odd), respectively; A_{n_even,m}, B_{n_even,m}, A_{n_odd,m} and B_{n_odd,m} are tri-diagonal or penta-diagonal matrices. ${\mathit{g}}_{\text{n\_even},m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and ${\mathit{g}}_{\text{n\_odd},m}^{\mathrm{c}\left(\mathrm{s}\right)}$ are calculated by

$$\begin{array}{}\text{(77)}& \begin{array}{rl}{\mathit{g}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{\mathrm{c}\left(\mathrm{s}\right)}& ={\mathbf{A}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{-\mathrm{1}}{\mathbf{B}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}{\mathit{f}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{\mathrm{c}\left(\mathrm{s}\right)},\\ {\mathit{g}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{\mathrm{c}\left(\mathrm{s}\right)}& ={\mathbf{A}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{-\mathrm{1}}{\mathbf{B}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}{\mathit{f}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{\mathrm{c}\left(\mathrm{s}\right)},\end{array}\end{array}$$

which can be solved efficiently as in Eq. (31). We have verified that all the eigenvalues of the matrices ${\mathbf{A}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{-\mathrm{1}}{\mathbf{B}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}$ and ${\mathbf{A}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{-\mathrm{1}}{\mathbf{B}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}$ are negative real numbers for several truncation wavenumbers M and N, but we have not yet proven that this is true for all truncation wavenumbers.

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

$$\begin{array}{}\text{(78)}& \begin{array}{rl}{\mathit{f}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{\mathrm{c}\left(\mathrm{s}\right)}& ={\mathbf{B}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{-\mathrm{1}}{\mathbf{A}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}{\mathit{g}}_{\mathrm{n}\mathrm{\_}\mathrm{even},m}^{\mathrm{c}\left(\mathrm{s}\right)},\\ {\mathit{f}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{\mathrm{c}\left(\mathrm{s}\right)}& ={\mathbf{B}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{-\mathrm{1}}{\mathbf{A}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}{\mathit{g}}_{\mathrm{n}\mathrm{\_}\mathrm{odd},m}^{\mathrm{c}\left(\mathrm{s}\right)},\end{array}\end{array}$$

except when m=0 and n is even. For m=0, $f_{n=\mathrm{0},m=\mathrm{0}}^{\mathrm{c}}$ disappears in Eq. (75a). The coefficients $f_{n,m=\mathrm{0}}^{\mathrm{c}}$ (even n≥2) are calculated from $g_{n,m=\mathrm{0}}^{\mathrm{c}}$ (even n≥2) by using Eq. (75a). The value $f_{n=\mathrm{0},m=\mathrm{0}}^{\mathrm{c}}$ is calculated from $f_{n,m=\mathrm{0}}^{\mathrm{c}}$ (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

$$\begin{array}{}\text{(79)}& \begin{array}{rl}f-\mathit{\epsilon }{\mathrm{\nabla }}^{\mathrm{2}}f& =\left\{\mathrm{1}-\mathit{\epsilon }{\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[{\displaystyle \frac{\mathrm{1}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial {\mathit{\lambda }}^{\mathrm{2}}}}\right.\right.\\ & \left.\left.+{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\right)\right]\right\}f=g,\end{array}\end{array}$$

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

$$\begin{array}{}\text{(80)}& \mathbf{A}\mathit{g}=\mathbf{B}\mathit{f},\end{array}$$

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

$$\begin{array}{}\text{(81)}& \mathbf{A}\mathit{f}-\mathit{\epsilon }\mathbf{B}\mathit{f}=\mathbf{A}\mathit{g}.\end{array}$$

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

$$\begin{array}{}\text{(82)}& \mathit{f}={\left(\mathbf{A}-\mathit{\epsilon }\mathbf{B}\right)}^{-\mathrm{1}}\mathbf{A}\mathit{g}.\end{array}$$

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

$$\begin{array}{}\text{(83)}& f-\mathit{\epsilon }{\mathrm{\nabla }}^{\mathrm{2}}f={\mathrm{\nabla }}^{\mathrm{2}}g,\end{array}$$

is represented as

$$\begin{array}{}\text{(84)}& \mathbf{A}\mathit{f}-\mathit{\epsilon }\mathbf{B}\mathit{f}=\mathbf{B}\mathit{g}.\end{array}$$

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

$$\begin{array}{}\text{(85)}& \mathit{f}={\left(\mathbf{A}-\mathit{\epsilon }\mathbf{B}\right)}^{-\mathrm{1}}\mathbf{B}\mathit{g}.\end{array}$$

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

$$\begin{array}{}\text{(86)}& f+\mathit{\epsilon }{\mathrm{\nabla }}^{\mathrm{4}}f=g,\end{array}$$

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

$$\begin{array}{}\text{(87)}& \left(\mathrm{1}+i\sqrt{\mathit{\epsilon }}{\mathrm{\nabla }}^{\mathrm{2}}\right)\left(\mathrm{1}-i\sqrt{\mathit{\epsilon }}{\mathrm{\nabla }}^{\mathrm{2}}\right)f=g,\end{array}$$

where $i=\sqrt{-\mathrm{1}}$. The calculation of Eq. (86) is accomplished by successive calculations of the following Helmholtz equations:

$$\begin{array}{}\text{(88a)}& {\displaystyle }\left(\mathrm{1}+i\sqrt{\mathit{\epsilon }}{\mathrm{\nabla }}^{\mathrm{2}}\right)f^{\prime }=g,\text{(88b)}& {\displaystyle }\left(\mathrm{1}-i\sqrt{\mathit{\epsilon }}{\mathrm{\nabla }}^{\mathrm{2}}\right)f=f^{\prime },\end{array}$$

which are represented as

$$\begin{array}{}\text{(89a)}& {\displaystyle }\left(\mathbf{A}+i\sqrt{\mathit{\epsilon }}\mathbf{B}\right){\mathit{f}}^{\prime }=\mathbf{A}\mathit{g},\text{(89b)}& {\displaystyle }\left(\mathbf{A}-i\sqrt{\mathit{\epsilon }}\mathbf{B}\right)\mathit{f}=\mathbf{A}{\mathit{f}}^{\prime }.\end{array}$$

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

$$\begin{array}{}\text{(90)}& \mathit{f}={\left(\mathbf{A}-i\sqrt{\mathit{\epsilon }}\mathbf{B}\right)}^{-\mathrm{1}}\mathbf{A}{\left(\mathbf{A}+i\sqrt{\mathit{\epsilon }}\mathbf{B}\right)}^{-\mathrm{1}}\mathbf{A}\mathit{g}.\end{array}$$

Here, $\mathbf{A}-i\sqrt{\mathit{\epsilon }}\mathbf{B}$ and $\mathbf{A}+i\sqrt{\mathit{\epsilon }}\mathbf{B}$ are complex matrices and f and g are real column vectors. For efficient computation, two real column vectors can be combined into one complex column vector (Cheong et al., 2004); for example, $\mathit{f}={\mathit{f}}^{\mathrm{c}}-i{\mathit{f}}^{\mathrm{s}}$ and $\mathit{g}={\mathit{g}}^{\mathrm{c}}-i{\mathit{g}}^{\mathrm{s}}$, where the superscript c indicates the zonal cosine component, and the superscript s indicates the zonal sine component.

3.10 Essential summary (cookbook) of the new DFS method

The essential summary for a scalar variable is as follows.

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

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

   • a.

     calculate the coefficients $T_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}$ from $T_{n,m}^{\text{c(s)}}$ by using Eqs. (14),

   • b.

     calculate $T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ in Eq. (13) from $T_{n,m}^{\mathrm{c}{\left(\mathrm{s}\right)}^{\prime }}$ by inverse cosine and sine Fourier transforms in Appendix B,

   • c.

     calculate the grid point values $T^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ in Eq. (15) from $T_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ by inverse Fourier transform.

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

   • a.

     calculate $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ in Eq. (9) from the grid point values T(λ,θ) by forward Fourier transform,

   • b.

     calculate the coefficients ${\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ in Eq. (23) from $T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ by forward cosine and sine transforms in Appendix B,

   • c.

     calculate the coefficients $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ from ${\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ by using Eqs. (30) and (31), where the coefficients $T_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ are calculated so that Eq. (29) derived from the Galerkin method (or the least-squares method) is satisfied.

The essential summary for a vector variable is as follows.

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

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

   • a.

     calculate the coefficients $u_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and $v_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ from ${\mathit{\chi }}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and ${\mathit{\psi }}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ by using Eq. (53) and the similar equations,

   • b.

     calculate $u_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ and $v_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ in Eq. (52) from $u_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and $v_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ by inverse cosine and sine transforms in Appendix B,

   • c.

     calculate the grid point values $u^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ and $v^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ in Eqs. (48a) and (49a) from $u_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ and $v_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)$ by inverse Fourier transform.

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

   • a.

     calculate $u_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ and $v_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ in Eq. (56) from the grid point values u(λ,θ) and v(λ,θ) by forward Fourier transform,

   • b.

     calculate the coefficients ${\stackrel{\mathrm{̃}}{u}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and ${\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ in Eq. (57) from the grid point values $u_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ and $v_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)$ by forward cosine and sine transforms in Appendix B,

   • c.

     calculate the coefficients ${\mathit{\chi }}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and ${\mathit{\psi }}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ from ${\stackrel{\mathrm{̃}}{u}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and ${\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ by using Eqs. (63) and (64), where ${\mathit{\chi }}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ and ${\mathit{\psi }}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ are calculated so that Eq. (62) derived from the Galerkin method (or the least-squares method) is satisfied.

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

Here we examine the error due to the meridional wavenumber truncation when the same continuity conditions at the poles as Eq. (3) are satisfied. In the DFS method of Orszag (1974) using Eq. (2), only $f_{N-\mathrm{1},m}$ and f_N,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 F_m(θ_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 F_m(θ_j).

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

In the DFS method of Orszag, a very large error occurs, especially for odd |m| (≥3) (Fig. 2c), when $f_{N-\mathrm{1},m}$ and f_N,m are modified to satisfy the pole conditions in Eq. (4). Dividing F_m(θ_j) by sin θ before the forward Fourier cosine transform for odd m also contributes to the large error.

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

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

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

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

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

Table 1 Normalized L₂ errors of Laplacian operator calculation (∇²f). We use the old DFS method with Grid [0]; the new DFS methods with Grid [0], Grid [1], and Grid [-1]; and the SH method. J⁰ is the number of latitudinal grid points in Grid [0]. The truncation wavenumber $N\cong \mathrm{2}{J}^{\mathrm{0}}/\mathrm{3}$.

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

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

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

$$\begin{array}{}\text{(91)}& \left(\mathrm{1}-\mathit{\epsilon }{\mathrm{\nabla }}^{\mathrm{2}}\right)f=h.\end{array}$$

Here, we give the function f as

$$\begin{array}{}\text{(92)}& f=\left\{\begin{array}{ll}\frac{H}{\mathrm{4}}{\left(\mathrm{1}+\cos \frac{\mathit{\pi }r}{R}\right)}^{\mathrm{2}}& \text{if }r<R,\\ \mathrm{0}& \text{if }r\ge R,\end{array}\right.\end{array}$$

$$\begin{array}{}\text{(93)}& r=a{\cos }^{-\mathrm{1}}\left[\sin {\mathit{\varphi }}_c\sin \mathit{\varphi }+\cos {\mathit{\varphi }}_c\cos \mathit{\varphi }\cos \left(\mathit{\lambda }-{\mathit{\lambda }}_c\right)\right],\end{array}$$

where $H=\mathrm{1000},R=a/\mathrm{3}$, ϕ is latitude, λ is longitude, a is the radius of the Earth, and r is the distance between (λ,ϕ) and the center $\left({\mathit{\lambda }}_c,{\mathit{\varphi }}_c\right)=\left(\mathrm{3}\mathit{\pi }/\mathrm{2},\mathit{\pi }/\mathrm{2}-\mathrm{0.05}\right)$. The function f is similar to the cosine bell in Williamson test case 1 (Williamson et al., 1992), but $\left(\mathrm{1}+\cos \mathit{\pi }r/R\right)$ is squared so that the second derivative of f is continuous. To easily calculate the exact values of ∇²f, the center is temporarily set to the North Pole, that is, $\left({\mathit{\lambda }}_c,{\mathit{\varphi }}_c\right)=\left(\mathrm{0},\mathit{\pi }/\mathrm{2}\right)$ and $r=a{\cos }^{-\mathrm{1}}\left[\sin \mathit{\varphi }\right]=a\mathit{\theta }$, where θ is colatitude. At this time, g is calculated as follows:

$$\begin{array}{}\text{(94)}& \begin{array}{rl}g& ={\mathrm{\nabla }}^{\mathrm{2}}f={\displaystyle \frac{\mathrm{1}}{a^{\mathrm{2}}}}\left[{\displaystyle \frac{\mathrm{1}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}{\displaystyle \frac{{\partial }^{\mathrm{2}}f}{\partial {\mathit{\lambda }}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\sin \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \mathit{\theta }{\displaystyle \frac{\partial f}{\partial \mathit{\theta }}}\right)\right]\\ & =-{\displaystyle \frac{\cos \mathit{\theta }}{\sin \mathit{\theta }}}{\displaystyle \frac{H}{\mathrm{2}{a}^{\mathrm{2}}}}{\displaystyle \frac{\mathit{\pi }a}{R}}\left[\left(\mathrm{1}+\cos {\displaystyle \frac{\mathit{\pi }r}{R}}\right)\sin {\displaystyle \frac{\mathit{\pi }r}{R}}\right]\\ & +{\displaystyle \frac{H}{\mathrm{2}{a}^{\mathrm{2}}}}{\left({\displaystyle \frac{\mathit{\pi }a}{R}}\right)}^{\mathrm{2}}\left[{\sin }^{\mathrm{2}}{\displaystyle \frac{\mathit{\pi }r}{R}}-\left(\mathrm{1}+\cos {\displaystyle \frac{\mathit{\pi }r}{R}}\right)\cos {\displaystyle \frac{\mathit{\pi }r}{R}}\right].\end{array}\end{array}$$

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

$$\begin{array}{}\text{(95)}& h=\left(\mathrm{1}-\mathit{\epsilon }{\mathrm{\nabla }}^{\mathrm{2}}\right)f=f-\mathit{\epsilon }g,\end{array}$$

where ε=0.01a² 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 ∇²f 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 ∇²f in Eq. (94). Here, the exact values of ∇²f 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 L₂ error between the calculated values and the exact values, which is normalized by the L₂ 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 ∇²f. The exact value of the global mean of ∇²f is zero. In Table 2, the global means calculated with each method are very close to zero. The global means of ∇²f 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 L₂ error between the calculated values and the exact values. The differences in error between the methods are small, and which method is better depends on the resolution and the arrangement of the grid points.

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

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

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

6.1 Shallow-water equations on a sphere

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

$$\begin{array}{}\text{(96)}& {\displaystyle \frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}}=-\mathrm{2}{\left(\mathbf{\Omega }\times \mathit{v}\right)}_{\text{H}}-g\mathrm{\nabla }h,\text{(97)}& {\displaystyle \frac{\mathrm{d}\left(h-h_{\mathrm{s}}\right)}{\mathrm{d}t}}=-\left(h-h_{\mathrm{s}}\right)\mathrm{\nabla }\cdot \mathit{v},\end{array}$$

where t is time, v is the horizontal wind vector, h is the height, h_s 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

$$\begin{array}{}\text{(98)}& {\displaystyle \frac{\mathrm{d}\left(\mathit{v}+\mathrm{2}\mathbf{\Omega }\times \mathit{r}\right)}{\mathrm{d}t}}=-g\mathrm{\nabla }h,\end{array}$$

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

$$\begin{array}{}\text{(99)}& {\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}=-\left(h-h_{\mathrm{s}}\right)\mathrm{\nabla }\cdot \mathit{v}+\mathit{v}\cdot \mathrm{\nabla }{h}_{\mathrm{s}}.\end{array}$$

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 J⁰=64, Δt=1800 s at about 120 km resolution with J⁰=160, Δt=1200 s at about 60 km resolution with J⁰=320, Δt=600 s at about 20 km resolution with J⁰=960, and Δt=90 s at about 1.3 km resolution with J⁰=15360, where J⁰ is the number of latitudinal grid points in Grid [0]. The number of latitudinal grid points J is J⁰ in Grid [0] (and in the Gaussian grid), J⁰+1 in Grid [1], and J⁰−1 in Grid [−1] (see Sect. 2). The number of longitudinal grid points I is 2J⁰. 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

$$\begin{array}{}\text{(100)}& {\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}={\displaystyle \frac{\partial h}{\partial t}}+\mathit{v}\cdot \mathrm{\nabla }h.\end{array}$$

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 M₀=20 is used in the DFS semi-Lagrangian shallow-water models. However, the larger the value M₀ is, the higher the longitudinal resolution around the pole becomes. Because of this, when the Eulerian scheme is used and M₀ is large, a time step must be very short due to the CFL condition. Therefore M₀ should be as small as possible. We have tested M₀=0, but this degrades the result of Williamson test case 1. We have also tested M₀=1, and this result is good. Therefore, we use M₀=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 J⁰=64. The meridional truncation wavenumber N and the zonal truncation wavenumber M are set as $N=M=\mathrm{42}\cong \mathrm{2}{J}^{\mathrm{0}}/\mathrm{3}$ because the two-thirds rule (Orszag, 1971) is used in order to avoid aliasing in the nonlinear advection term. The time step is 30 min. The angle between the solid body rotation and the polar axis α is  $\mathit{\pi }/\mathrm{2}-\mathrm{0.05}$. The results for DFS [0], DFS [1], DFS [−1], and SH are very similar. Instability occurs in the old DFS model without horizontal diffusion. This is probably because of the appearance of high-wavenumber oscillations due to the wavenumber truncation with $N\cong \mathrm{2}{J}^{\mathrm{0}}/\mathrm{3}$ for even m (≠0) in the old DFS method, as shown in Sect. 4. Table 4 shows the normalized L₂ 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 J⁰=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 M₀=1 than in the SH model. When the fourth-order horizontal diffusion in Eq. (86) with $\mathit{\epsilon }=a^{\mathrm{4}}\mathrm{\Delta }t/\left(\mathrm{7.2}\times \mathrm{3600}\times {\mathrm{107}}^{\mathrm{2}}{N}^{\mathrm{2}}\right)$ 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 4 Normalized L₂ errors of the predicted height after a 12 d integration in Williamson test case 1 when using the Eulerian advection models. The truncation wavenumber $N\cong \mathrm{2}{J}^{\mathrm{0}}/\mathrm{3}$.

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

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 J⁰= 64, the errors in the semi-Lagrangian models are larger than those in the Eulerian models, but at the resolutions J⁰=160, 320, and 960, the errors in the semi-Lagrangian models are smaller than those in the Eulerian models.

Table 5 The same as Table 4 but using the semi-Lagrangian models and the truncation wavenumber $N\cong J^{\mathrm{0}}-\mathrm{1}$.

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

$$\begin{array}{}\text{(101)}& \begin{array}{rl}\text{2}\mathbf{\Omega }\times \mathit{r}& =\left(\mathrm{2}\mathrm{\Omega }a\left[\cos \mathit{\theta }\cos \mathit{\alpha }+\cos \mathit{\lambda }\sin \mathit{\theta }\sin \mathit{\alpha }\right],\right.\\ & \left.-\mathrm{2}\mathrm{\Omega }a\sin \mathit{\lambda }\sin \mathit{\alpha }\right).\end{array}\end{array}$$

Figure 5 shows the time series of forecast errors of the height for a 5 d integration in Williamson test case 2 with $\mathit{\alpha }=\mathit{\pi }/\mathrm{2}-\mathrm{0.05}$ in the models at the resolution J⁰=64 and N=63 (DFS) or N=62 (SH) using no horizontal diffusion. The normalized L₁, L₂, 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 L₂ errors of the predicted height after a 5 d integration. The errors are almost the same between the old DFS, new DFS, and SH models at each resolution.

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

Figure 5 Time series of prediction error of height (m) for 5 d (120 h) integration in Williamson test case 2 $\left(\mathit{\alpha }=\mathit{\pi }/\mathrm{2}-\mathrm{0.05}\right)$. The number of longitudinal grid points is I=128. The number of latitudinal grid points in Grid [0] is J⁰=64. The truncation wavenumber is N=63. Solid, dashed, and dotted lines represent normalized L₁, L₂, and L_∞ errors, respectively. The colors blue, green, red, purple, and orange represent the models using SH, old DFS with Grid [0], new DFS with Grid [0], new DFS with Grid [1], and new DFS with Grid [−1], respectively.

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 h₀=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 J⁰=64. Table 7 shows the normalized L₂ errors of the predicted height after a 15 d integration. The errors are almost the same between the old DFS, new DFS, and SH models at each resolution. The errors do not decrease when the resolution increases, which is different from the results in the other test cases. This may be because the mountain topography is not a differentiable function, and the mountain is added impulsively on to an initially balanced flow (Galewsky et al., 2004).

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

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

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

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 J⁰=64 and J⁰=960. While the zonal wavenumber 1 component is dominant in the new DFS model at the resolution J⁰=64 and N=63, high zonal wavenumber noise appears in the old DFS model at the same resolution. One possible reason is that the latitudinal derivative at the grid points can be unrealistically large in the old DFS method, even when $N^{\prime }=J^{\mathrm{0}}$ 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 ²θsin nθ. In the old DFS model at high resolution with J⁰=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 J⁰=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 J⁰=960, the high-wavenumber components are a little larger than in the other models, but the differences are slight.

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

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

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

Table 8 The same as Table 7 but with $N\cong \mathrm{2}{J}^{\mathrm{0}}/\mathrm{3}$.

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

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 $N\cong \mathrm{2}{J}^{\mathrm{0}}/\mathrm{3}$. At the resolution J⁰=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 J⁰=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 L₂ errors of the predicted height after a 14 d integration. The error in the new DFS model using Grid [1] is the smallest, and the error in the new DFS model using Grid [−1] is the second smallest at each resolution. The errors in the old and new DFS models using Grid [0] and in the SH model are very close.

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

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

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 J⁰=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.

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

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 J⁰=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 J⁰=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(N³) operations, whereas the Fourier cosine and sine transforms used in the DFS model require only O(N²log N) operations. If the fast Legendre transform, which requires only N²(log N)³ 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.

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

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

7 Conclusions and perspectives

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

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

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

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

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

The high zonal wavenumber noise of the meridional wind appears near the poles in the old DFS model, but not in the new DFS models in Williamson test case 5. One possible reason is that the latitudinal derivative at the grid points can be unrealistically large in the old DFS method even when the truncation wavenumber N^′ for even m (≠0) is equal to the number of latitudinal grid points J, while the new DFS expansion method with the least-squares method does not have this problem. In the old DFS model, $N^{\prime }<J$ for even m (≠0) causes numerical instability. In the new DFS model, an arbitrary meridional wavenumber truncation N (<J) can be used without the stability problem because the error due to meridional wavenumber truncation is small when using the new DFS expansion method with the least-squares method. This is one of the merits of the new DFS method because quadratic truncation $\left(N\cong \mathrm{2}J/\mathrm{3}\right)$ or cubic truncation $\left(N\cong J/\mathrm{2}\right)$ are usually used in the Eulerian model and are also starting to be used in the semi-Lagrangian model instead of the linear truncation $\left(N\cong J-\mathrm{1}\right)$ 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(N²log N) operations, while the Legendre transform in the SH model requires O(N³) operations. We have executed our shallow-water models on Intel CPUs. The execution on GPUs is one important topic, but we have not tested our models in this way because execution on GPUs is not an easy task. MPI parallelization is another important topic. However, in our shallow-water models, we use only OpenMP parallelization and not MPI parallelization due to the simplicity of the source code.

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

Appendix A: Trigonometric identities

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

$$\begin{array}{}\text{(A1a)}& {\displaystyle }\sin n\mathit{\theta }\cos n^{\prime }\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\sin \left(n+n^{\prime }\right)\mathit{\theta }+\sin \left(n-n^{\prime }\right)\mathit{\theta }\right],\text{(A1b)}& {\displaystyle }\cos n\mathit{\theta }\sin n^{\prime }\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\sin \left(n+n^{\prime }\right)\mathit{\theta }-\sin \left(n-n^{\prime }\right)\mathit{\theta }\right],\text{(A1c)}& {\displaystyle }\cos n\mathit{\theta }\cos n^{\prime }\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\cos \left(n+n^{\prime }\right)\mathit{\theta }+\cos \left(n-n^{\prime }\right)\mathit{\theta }\right],\text{(A1d)}& {\displaystyle }\sin n\mathit{\theta }\sin n^{\prime }\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[-\cos \left(n+n^{\prime }\right)\mathit{\theta }+\cos \left(n-n^{\prime }\right)\mathit{\theta }\right].\end{array}$$

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

$$\begin{array}{}\text{(A2a)}& {\displaystyle }\sin \mathit{\theta }\cos n\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\sin \left(n+\mathrm{1}\right)\mathit{\theta }-\sin \left(n-\mathrm{1}\right)\mathit{\theta }\right],\text{(A2b)}& {\displaystyle }\sin \mathit{\theta }\sin n\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[-\cos \left(n+\mathrm{1}\right)\mathit{\theta }+\cos \left(n-\mathrm{1}\right)\mathit{\theta }\right],\text{(A2c)}& {\displaystyle }\begin{array}{rl}{\sin }^{\mathrm{2}}\mathit{\theta }\sin n\mathit{\theta }& ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[-\sin \left(n-\mathrm{2}\right)\mathit{\theta }+\mathrm{2}\sin n\mathit{\theta }\right.\\ & \left.-\sin \left(n+\mathrm{2}\right)\mathit{\theta }\right],\end{array}\text{(A2d)}& {\displaystyle }\begin{array}{rl}{\sin }^{\mathrm{2}}\mathit{\theta }\cos n\mathit{\theta }& ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[-\cos \left(n-\mathrm{2}\right)\mathit{\theta }+\mathrm{2}\cos n\mathit{\theta }\right.\\ & \left.-\cos \left(n+\mathrm{2}\right)\mathit{\theta }\right],\end{array}\text{(A2e)}& {\displaystyle }\begin{array}{rl}{\sin }^{\mathrm{4}}\mathit{\theta }\sin n\mathit{\theta }& ={\displaystyle \frac{\mathrm{1}}{\mathrm{16}}}\left[\sin \left(n-\mathrm{4}\right)\mathit{\theta }-\mathrm{4}\sin \left(n-\mathrm{2}\right)\mathit{\theta }\right.\\ & +\mathrm{6}\sin n\mathit{\theta }-\mathrm{4}\sin \left(n+\mathrm{2}\right)\mathit{\theta }\\ & \left.+\sin \left(n+\mathrm{4}\right)\mathit{\theta }\right].\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A3a)}& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\cos m\mathit{\lambda }\cos m^{\prime }\mathit{\lambda }\mathrm{d}\mathit{\lambda }=\left\{\begin{array}{ll}\mathrm{2}\mathit{\pi }& \text{ for}m=m^{\prime }=\mathrm{0},\\ \mathit{\pi }& \text{ for}m=m^{\prime }\ne \mathrm{0},\\ \mathrm{0}& \text{ for}m\ne m^{\prime },\end{array}\right.\text{(A3b)}& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\cos m\mathit{\lambda }\sin m^{\prime }\mathit{\lambda }\mathrm{d}\mathit{\lambda }=\mathrm{0},\text{(A3c)}& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\sin m\mathit{\lambda }\sin m^{\prime }\mathit{\lambda }\mathrm{d}\mathit{\lambda }=\left\{\begin{array}{ll}\mathit{\pi }& \text{ for }m=m^{\prime }\ne \mathrm{0},\\ \mathrm{0}& \text{ for }m\ne m^{\prime }.\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(A4a)}& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\cos n\mathit{\theta }\cos n^{\prime }\mathit{\theta }\mathrm{d}\mathit{\theta }=\left\{\begin{array}{ll}\mathit{\pi }& \text{ for}n=n^{\prime }=\mathrm{0},\\ \frac{\mathrm{1}}{\mathrm{2}}\mathit{\pi }& \text{ for}n=n^{\prime }\ne \mathrm{0},\\ \mathrm{0}& \text{ for}n\ne n^{\prime },\end{array}\right.\text{(A4b)}& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\sin n\mathit{\theta }\sin n^{\prime }\mathit{\theta }\mathrm{d}\mathit{\theta }=\left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\pi }& \text{ for }n=n^{\prime }\ne \mathrm{0},\\ \mathrm{0}& \text{ for }n\ne n^{\prime }.\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(A5a)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left({\sin }^l\mathit{\theta }\cos n\mathit{\theta }\right)={\displaystyle \frac{n+l}{\mathrm{2}}}{\sin }^{l-\mathrm{1}}\mathit{\theta }\cos \left(n+\mathrm{1}\right)\mathit{\theta }\\ & -{\displaystyle \frac{n-l}{\mathrm{2}}}{\sin }^{l-\mathrm{1}}\mathit{\theta }\cos \left(n-\mathrm{1}\right)\mathit{\theta },\end{array}\text{(A5b)}& {\displaystyle }\begin{array}{rl}& \sin \mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left[\sin \mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left({\sin }^l\mathit{\theta }\cos n\mathit{\theta }\right)\right]\\ & ={\displaystyle \frac{\left(n+l\right)\left(n+l+\mathrm{1}\right)}{\mathrm{4}}}{\sin }^l\mathit{\theta }\cos \left(n+\mathrm{2}\right)\mathit{\theta }\\ & -{\displaystyle \frac{\mathrm{2}{n}^{\mathrm{2}}-\mathrm{2}{l}^{\mathrm{2}}+\mathrm{2}l}{\mathrm{4}}}{\sin }^l\mathit{\theta }\cos n\mathit{\theta }\\ & +{\displaystyle \frac{\left(n-l\right)\left(n-l-\mathrm{1}\right)}{\mathrm{4}}}{\sin }^l\mathit{\theta }\cos \left(n-\mathrm{2}\right)\mathit{\theta },\end{array}\text{(A5c)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left({\sin }^l\mathit{\theta }\sin n\mathit{\theta }\right)={\displaystyle \frac{n+l}{\mathrm{2}}}{\sin }^{l-\mathrm{1}}\mathit{\theta }\sin \left(n+\mathrm{1}\right)\mathit{\theta }\\ & -{\displaystyle \frac{n-l}{\mathrm{2}}}{\sin }^{l-\mathrm{1}}\mathit{\theta }\sin \left(n-\mathrm{1}\right)\mathit{\theta },\end{array}\text{(A5d)}& {\displaystyle }\begin{array}{rl}& \sin \mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left[\sin \mathit{\theta }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left({\sin }^l\mathit{\theta }\sin n\mathit{\theta }\right)\right]\\ & ={\displaystyle \frac{\left(n+l\right)\left(n+l+\mathrm{1}\right)}{\mathrm{4}}}{\sin }^l\mathit{\theta }\sin \left(n+\mathrm{2}\right)\mathit{\theta }\\ & -{\displaystyle \frac{\mathrm{2}{n}^{\mathrm{2}}-\mathrm{2}{l}^{\mathrm{2}}+\mathrm{2}l}{\mathrm{4}}}{\sin }^l\mathit{\theta }\sin n\mathit{\theta }\\ & +{\displaystyle \frac{\left(n-l\right)\left(n-l-\mathrm{1}\right)}{\mathrm{4}}}{\sin }^l\mathit{\theta }\sin \left(n-\mathrm{2}\right)\mathit{\theta }.\end{array}\end{array}$$

Appendix B: Discrete Fourier cosine and sine transforms in latitude

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

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

$$\begin{array}{}\text{(B1a)}& {\displaystyle }g\left({\mathit{\theta }}_j\right)=\sum _{n=\mathrm{0}}^{J^{\mathrm{0}}-\mathrm{1}}{g}_n\cos n{\mathit{\theta }}_j,\text{(B1b)}& {\displaystyle }\begin{array}{rl}{g}_n& ={\displaystyle \frac{b}{J^{\mathrm{0}}}}\sum _{j=\mathrm{0}}^{J^{\mathrm{0}}-\mathrm{1}}g\left({\mathit{\theta }}_j\right)\cos n{\mathit{\theta }}_j,\\ & b\equiv \left\{\begin{array}{ll}\mathrm{1}& \text{ for }n=\mathrm{0}\\ \mathrm{2}& \text{ for }\mathrm{1}\le n\le J^{\mathrm{0}}-\mathrm{1}.\end{array}\right.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B2a)}& {\displaystyle }h\left({\mathit{\theta }}_j\right)=\sum _{n=\mathrm{1}}^{J^{\mathrm{0}}}{h}_n\sin n{\mathit{\theta }}_j,\text{(B2b)}& {\displaystyle }\begin{array}{rl}{h}_n& ={\displaystyle \frac{b}{J^{\mathrm{0}}}}\sum _{j=\mathrm{0}}^{J^{\mathrm{0}}-\mathrm{1}}h\left({\mathit{\theta }}_j\right)\sin n{\mathit{\theta }}_j,\\ & b\equiv \left\{\begin{array}{ll}\mathrm{1}& \text{ for}n=J^{\mathrm{0}}\\ \mathrm{2}& \text{ for}\mathrm{1}\le n\le J^{\mathrm{0}}-\mathrm{1}.\end{array}\right.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B3a)}& {\displaystyle }g\left({\mathit{\theta }}_j\right)=\sum _{n=\mathrm{0}}^{J^{\mathrm{0}}}{g}_n\cos n{\mathit{\theta }}_j,\text{(B3b)}& {\displaystyle }\begin{array}{rl}{g}_n& ={\displaystyle \frac{b}{J^{\mathrm{0}}}}\sum _{j=\mathrm{0}}^{J^{\mathrm{0}}}cg\left({\mathit{\theta }}_j\right)\cos n{\mathit{\theta }}_j,\\ & b\equiv \left\{\begin{array}{ll}\mathrm{1}& \text{ for}n=\mathrm{0},J^{\mathrm{0}}\\ \mathrm{2}& \text{ for}\mathrm{1}\le n\le J^{\mathrm{0}}-\mathrm{1},\end{array}\right.\\ & c\equiv \left\{\begin{array}{ll}\mathrm{1}/\mathrm{2}& \text{ for}j=\mathrm{0},J^{\mathrm{0}}\\ \mathrm{1}& \text{ for}\mathrm{1}\le j\le J^{\mathrm{0}}-\mathrm{1}.\end{array}\right.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B4a)}& {\displaystyle }h\left({\mathit{\theta }}_j\right)=\sum _{n=\mathrm{1}}^{J^{\mathrm{0}}-\mathrm{1}}{h}_n\sin n{\mathit{\theta }}_j,h\left({\mathit{\theta }}_{\mathrm{0}}\right)=h\left({\mathit{\theta }}_{J^{\mathrm{0}}}\right)=\mathrm{0},\text{(B4b)}& {\displaystyle }{h}_n={\displaystyle \frac{\mathrm{2}}{J^{\mathrm{0}}}}\sum _{j=\mathrm{1}}^{J^{\mathrm{0}}-\mathrm{1}}h\left({\mathit{\theta }}_j\right)\sin n{\mathit{\theta }}_j\left(\mathrm{1}\le n\le J^{\mathrm{0}}-\mathrm{1}\right).\end{array}$$

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:

$$\begin{array}{}\text{(B5)}& g\left({\mathit{\theta }}_j\right)=\sum _{n=\mathrm{0}}^{J^{\mathrm{0}}-\mathrm{2}}{g}_n\cos n{\mathit{\theta }}_j,\end{array}$$

where n is from 0 to $J^{\mathrm{0}}-\mathrm{2}\left(=J-\mathrm{1}\right)$ because the number of the meridional grid points is $J^{\mathrm{0}}-\mathrm{1}\left(=J\right)$ in Grid [−1]. However, the forward discrete cosine transform cannot be performed like Eq. (B3b). We can calculate the expansion coefficients g_n from g(θ_j) in the following way. Equation (B5) is multiplied by sin θ_j, and we define $\widehat{g}\left({\mathit{\theta }}_j\right)$ as follows:

$$\begin{array}{}\text{(B6)}& \widehat{g}\left({\mathit{\theta }}_j\right)\equiv g\left({\mathit{\theta }}_j\right)\sin {\mathit{\theta }}_j=\sum _{n=\mathrm{0}}^{J^{\mathrm{0}}-\mathrm{2}}{g}_n\sin {\mathit{\theta }}_j\cos n{\mathit{\theta }}_j.\end{array}$$

We can expand $\widehat{g}\left({\mathit{\theta }}_j\right)$ as follows:

$$\begin{array}{}\text{(B7)}& \widehat{g}\left({\mathit{\theta }}_j\right)=\sum _{n=\mathrm{1}}^{J^{\mathrm{0}}-\mathrm{1}}{\widehat{g}}_n\sin n{\mathit{\theta }}_j.\end{array}$$

The expansion coefficients ${\widehat{g}}_n$ can be obtained from $\widehat{g}\left({\mathit{\theta }}_j\right)$ in the same way as in Eq. (B4b) by forward discrete sine transform:

$$\begin{array}{}\text{(B8)}& {\widehat{g}}_n={\displaystyle \frac{\mathrm{2}}{J^{\mathrm{0}}}}\sum _{j=\mathrm{1}}^{J^{\mathrm{0}}-\mathrm{1}}\widehat{g}\left({\mathit{\theta }}_j\right)\sin n{\mathit{\theta }}_j.\end{array}$$

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

$$\begin{array}{}\text{(B9)}& \sum _{n=\mathrm{0}}^{J^{\mathrm{0}}-\mathrm{2}}{g}_n\sin \mathit{\theta }\cos n\mathit{\theta }=\sum _{n=\mathrm{1}}^{J^{\mathrm{0}}-\mathrm{1}}{\widehat{g}}_n\sin n\mathit{\theta }.\end{array}$$

By using Eq. (A2a), we obtain

$$\begin{array}{}\text{(B10)}& \begin{array}{rl}\sum _{n=\mathrm{0}}^{J^{\mathrm{0}}-\mathrm{2}}& g_n\sin \mathit{\theta }\cos n\mathit{\theta }=\left(g_{\mathrm{0}}-{\displaystyle \frac{g_{\mathrm{2}}}{\mathrm{2}}}\right)\sin \mathit{\theta }\\ & +\sum _{n=\mathrm{2}}^{J^{\mathrm{0}}-\mathrm{3}}\left({\displaystyle \frac{g_{n-\mathrm{1}}}{\mathrm{2}}}-{\displaystyle \frac{g_{n+\mathrm{1}}}{\mathrm{2}}}\right)\sin n\mathit{\theta }\\ & +{\displaystyle \frac{g_{J^{\mathrm{0}}-\mathrm{3}}}{\mathrm{2}}}\sin \left(J^{\mathrm{0}}-\mathrm{2}\right)\mathit{\theta }\\ & +{\displaystyle \frac{g_{J^{\mathrm{0}}-\mathrm{2}}}{\mathrm{2}}}\sin \left(J^{\mathrm{0}}-\mathrm{1}\right)\mathit{\theta }.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B11)}& {\widehat{g}}_n=\left\{\begin{array}{ll}{g}_{\mathrm{0}}-\frac{g_{\mathrm{2}}}{\mathrm{2}}& \text{ for}n=\mathrm{1},\\ \frac{g_{n-\mathrm{1}}}{\mathrm{2}}-\frac{g_{n+\mathrm{1}}}{\mathrm{2}}& \text{ for}\mathrm{2}\le n\le J^{\mathrm{0}}-\mathrm{3},\\ \frac{g_{J^{\mathrm{0}}-\mathrm{3}}}{\mathrm{2}}& \text{ for}n=J^{\mathrm{0}}-\mathrm{2},\\ \frac{g_{J^{\mathrm{0}}-\mathrm{2}}}{\mathrm{2}}& \text{ for}n=J^{\mathrm{0}}-\mathrm{1}.\end{array}\right.\end{array}$$

We can calculate $\widehat{g}\left({\mathit{\theta }}_j\right)$ from g(θ_j) using Eq. (B6), calculate ${\widehat{g}}_n$ from $\widehat{g}\left({\mathit{\theta }}_j\right)$ using Eq. (B8), and calculate g_n from ${\widehat{g}}_n$ using Eq. (B11).

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

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

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

${\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}$ in Eq. (23) is calculated by the forward Fourier cosine or sine transform as follows:

$$\begin{array}{}\text{(D1)}& {\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}=\left\{\begin{array}{l}\frac{b}{\mathit{\pi }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\cos n\mathit{\theta }{T}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta },\\ b\equiv \left\{\begin{array}{ll}\mathrm{1}& \text{for}n=\mathrm{0}\\ \mathrm{2}& \text{for}n\ne \mathrm{0},\end{array}\right.\text{for\hspace{0.17em} even\hspace{0.17em} }m,\\ \frac{\mathrm{2}}{\mathit{\pi }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\sin n\mathit{\theta }{T}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\text{for\hspace{0.17em} odd}m.\\ \end{array}\right.\end{array}$$

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:

$$\begin{array}{}\text{(D2a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{b}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \cos n\mathit{\theta }{\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }={\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & \left(n=\mathrm{0},\mathrm{\dots },N\right)\text{for\hspace{0.17em} even}m\end{array}\text{(D2b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \sin n\mathit{\theta }{\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }={\stackrel{\mathrm{̃}}{T}}_{n,m}^{\mathrm{c}\left(\mathrm{s}\right)}\\ & \left(n=\mathrm{1},\mathrm{\dots },N\right)\text{for\hspace{0.17em} odd}m.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(D3a)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \cos n\mathit{\theta }{T}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\cos n\mathit{\theta }{\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & \left(n=\mathrm{0},\mathrm{\dots },N\right)\text{for\hspace{0.17em} even}m,\end{array}\text{(D3b)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \sin n\mathit{\theta }{T}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\sin n\mathit{\theta }{\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & \left(n=\mathrm{1},\mathrm{\dots },N\right)\text{for\hspace{0.17em} odd}m.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(D4)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& S_{n,m}\left(\mathit{\theta }\right)T_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}\left(\mathit{\theta }\right){\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & \left(n=N_{\text{min,\hspace{0.17em} }m},\mathrm{\dots },N_{\text{max,}m}\right).\end{array}\end{array}$$

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):

$$\begin{array}{}\text{(D5a)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \sin n\mathit{\theta }{u}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\sin n\mathit{\theta }{\stackrel{\mathrm{̃}}{u}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & \left(n=\mathrm{1},\mathrm{\dots },N\right)\text{for\hspace{0.17em} even}m,\end{array}\text{(D5b)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& \cos n\mathit{\theta }{u}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\cos n\mathit{\theta }{\stackrel{\mathrm{̃}}{u}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & \left(n=\mathrm{0},\mathrm{\dots },N\right)\text{for\hspace{0.17em} odd}m\end{array}\end{array}$$

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

$$\begin{array}{}\text{(D6a)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& {\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{\sin \mathit{\theta }}}{u}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{\displaystyle \frac{m{S}_{n,m}\left(\mathit{\theta }\right)}{\sin \mathit{\theta }}}{\stackrel{\mathrm{̃}}{u}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & \left(n=N_{\text{min,\hspace{0.17em} }m},\mathrm{\dots },N_{\text{max,}m}\right),\end{array}\text{(D6b)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}& {\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}{u}_m^{\mathrm{c}\left(\mathrm{s}\right)}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{\displaystyle \frac{\partial S_{n,m}\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}{\stackrel{\mathrm{̃}}{u}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & \left(n=N_{\text{min,\hspace{0.17em} }m},\mathrm{\dots },N_{\text{max,}m}\right).\end{array}\end{array}$$

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

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

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

$$\begin{array}{}\text{(E1)}& {\displaystyle }\begin{array}{rl}& \left(\text{l.h.s\hspace{0.17em} of\hspace{0.17em} Eq.\hspace{0.17em} 29\hspace{0.17em} for\hspace{0.17em} odd\hspace{0.17em} }m\ge \mathrm{3}\right)\\ & =\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}\left(\mathit{\theta }\right)T_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & =\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{\sin }^{\mathrm{2}}\mathit{\theta }\sin n\mathit{\theta }\sum _{n^{\prime }=\mathrm{1}}^{N-\mathrm{2}}{T}_{n^{\prime },m}^{\mathrm{c}}{\sin }^{\mathrm{2}}\mathit{\theta }\sin n^{\prime }\mathit{\theta }\mathrm{d}\mathit{\theta }\\ & =\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\sin n\mathit{\theta }\sum _{n^{\prime }=\mathrm{1}}^{N-\mathrm{2}}{\displaystyle \frac{T_{n^{\prime },m}^{\mathrm{c}}}{\mathrm{16}}}\left[\sin \left(n^{\prime }-\mathrm{4}\right)\mathit{\theta }-\mathrm{4}\sin \left(n^{\prime }-\mathrm{2}\right)\mathit{\theta }\right.\\ & \left.+\mathrm{6}\sin n^{\prime }\mathit{\theta }-\mathrm{4}\sin \left(n^{\prime }+\mathrm{2}\right)\mathit{\theta }+\sin \left(n^{\prime }+\mathrm{4}\right)\mathit{\theta }\right]\mathrm{d}\mathit{\theta }\\ & =\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\sin n\mathit{\theta }\left[{\displaystyle \frac{\mathrm{10}{T}_{\mathrm{1},m}^{\mathrm{c}}-\mathrm{5}{T}_{\mathrm{3},m}^{\mathrm{c}}+T_{\mathrm{5},m}^{\mathrm{c}}}{\mathrm{16}}}\sin \mathit{\theta }\right.\\ & +{\displaystyle \frac{\mathrm{5}{T}_{\mathrm{2},m}^{\mathrm{c}}-\mathrm{4}{T}_{\mathrm{4},m}^{\mathrm{c}}+T_{\mathrm{6},m}^{\mathrm{c}}}{\mathrm{16}}}\sin \mathrm{2}\mathit{\theta }\\ & +{\displaystyle \frac{-\mathrm{5}{T}_{\mathrm{1},m}^{\mathrm{c}}+\mathrm{6}{T}_{\mathrm{3},m}^{\mathrm{c}}-\mathrm{4}{T}_{\mathrm{5},m}^{\mathrm{c}}+T_{\mathrm{7},m}^{\mathrm{c}}}{\mathrm{16}}}\sin \mathrm{3}\mathit{\theta }\\ & +\sum _{n^{\prime }=\mathrm{4}}^{N+\mathrm{2}}{\displaystyle \frac{T_{n^{\prime }-\mathrm{4},m}^{\mathrm{c}}-\mathrm{4}{T}_{n^{\prime }-\mathrm{2},m}^{\mathrm{c}}+\mathrm{6}{T}_{n^{\prime },m}^{\mathrm{c}}-\mathrm{4}{T}_{n^{\prime }+\mathrm{2},m}^{\mathrm{c}}+T_{n^{\prime }+\mathrm{4},m}^{\mathrm{c}}}{\mathrm{16}}}\\ & \left.\sin n^{\prime }\mathit{\theta }\right]\mathrm{d}\mathit{\theta }\end{array}\text{(E2)}& {\displaystyle }\begin{array}{rl}& \left(\text{r.h.s\hspace{0.17em} of\hspace{0.17em} Eq.\hspace{0.17em} (29)\hspace{0.17em} for\hspace{0.17em} odd\hspace{0.17em} }m\ge \mathrm{3}\right)=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{S}_{n,m}\left(\mathit{\theta }\right){\stackrel{\mathrm{̃}}{T}}_m^{\mathrm{c}\left(\mathrm{s}\right),N}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & =\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{\sin }^{\mathrm{2}}\mathit{\theta }\sin n\mathit{\theta }\sum _{n^{\prime }=\mathrm{1}}^N{\stackrel{\mathrm{̃}}{T}}_{n^{\prime },m}^{\mathrm{c}}\sin n^{\prime }\mathit{\theta }\mathrm{d}\mathit{\theta },\\ & =\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\sin n\mathit{\theta }\sum _{n^{\prime }=\mathrm{1}}^N{\displaystyle \frac{{\stackrel{\mathrm{̃}}{T}}_{n^{\prime },m}^{\mathrm{c}}}{\mathrm{4}}}\left[-\sin \left(n^{\prime }-\mathrm{2}\right)\mathit{\theta }+\mathrm{2}\sin n^{\prime }\mathit{\theta }\right.\\ & \left.-\sin \left(n^{\prime }+\mathrm{2}\right)\mathit{\theta }\right]\mathrm{d}\mathit{\theta }\\ & =\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\sin n\mathit{\theta }\left[{\displaystyle \frac{\mathrm{3}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{1},m}^{\mathrm{c}}-{\stackrel{\mathrm{̃}}{T}}_{\mathrm{3},m}^{\mathrm{c}}}{\mathrm{4}}}\sin \mathit{\theta }+{\displaystyle \frac{\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{2},m}^{\mathrm{c}}-{\stackrel{\mathrm{̃}}{T}}_{\mathrm{4},m}^{\mathrm{c}}}{\mathrm{4}}}\sin \mathrm{2}\mathit{\theta }\right.\\ & \left.+\sum _{n^{\prime }=\mathrm{3}}^{N+\mathrm{2}}{\displaystyle \frac{-{\stackrel{\mathrm{̃}}{T}}_{n^{\prime }-\mathrm{2},m}^{\mathrm{c}}+\mathrm{2}{\stackrel{\mathrm{̃}}{T}}_{n^{\prime },m}^{\mathrm{c}}-{\stackrel{\mathrm{̃}}{T}}_{n^{\prime }+\mathrm{2},m}^{\mathrm{c}}}{\mathrm{4}}}\sin n^{\prime }\mathit{\theta }\right]\mathrm{d}\mathit{\theta },\end{array}\end{array}$$

where $\mathrm{1}\le n\le N-\mathrm{2}$. From Eqs. (29), (E1), (E2), and (A4b), Eq. (30d) is derived.

Appendix F: Zonal Fourier filter

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

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

$$\begin{array}{}\text{(F1)}& M_{\text{f}}\left({\mathit{\theta }}_j\right)=min\left(M,M_{\mathrm{0}}+M\sin \left({\mathit{\theta }}_j\right)\right).\end{array}$$

The values of $T_m^{\mathrm{c}}\left({\mathit{\theta }}_j\right)$ and $T_m^{\mathrm{s}}\left({\mathit{\theta }}_j\right)$ in Eq. (8) are set to zero for m>M_f(θ_j) during the spectral transform. We use the value M₀=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 M₀=1 as described in Sect. 6.3.

Appendix G: Calculation of global mean and latitudinal area weight

The global mean value of $T^{N,M}\left(\mathit{\lambda },\mathit{\theta }\right)$ in Eq. (15) can be calculated in spectral space by the following equation (Cheong 2000a):

$$\begin{array}{}\text{(G1)}& \begin{array}{rl}G& ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\left(\sum _{m=\mathrm{0}}^M{T}_m^{\mathrm{c},N}\left(\mathit{\theta }\right)\cos m\mathit{\lambda }+\sum _{m=\mathrm{1}}^M{T}_m^{\mathrm{s},N}\left(\mathit{\theta }\right)\sin m\mathit{\lambda }\right)\\ & \sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\lambda }\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\sum _{n=\mathrm{0}}^N{T}_{n,m=\mathrm{0}}^{\mathrm{c}}\cos n\mathit{\theta }\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\\ & =\sum _{\begin{array}{c}n=\mathrm{0}\\ \text{when}n\text{ is\hspace{0.17em} even}\end{array}}^N{\displaystyle \frac{T_{n,m=\mathrm{0}}^{\mathrm{c}}}{\mathrm{1}-n^{\mathrm{2}}}},\end{array}\end{array}$$

where Eq. (A2a) is used.

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

1. The latitudinal distribution of $T_{m=\mathrm{0}}^{\mathrm{c}\left(j\right)}\left({\mathit{\theta }}_k\right)$ for each j is given as

   $$\begin{array}{}\text{(G2)}& T_{m=\mathrm{0}}^{\mathrm{c}\left(j\right)}\left({\mathit{\theta }}_k\right)=\left\{\begin{array}{ll}\mathrm{1}& \text{ for}k=j\\ \mathrm{0}& \text{ for}k\ne j,\end{array}\right.\end{array}$$

   where $\mathrm{0}\le k\le J-\mathrm{1}$ in Grid [0] and Grid [1] and $\mathrm{1}\le k\le J$ in Grid [−1] (see Sect. 2).

2. From $T_{m=\mathrm{0}}^{\mathrm{c}\left(j\right)}\left({\mathit{\theta }}_k\right)$, the meridional expansion coefficients $T_{n,m=\mathrm{0}}^{\mathrm{c}\left(j\right)}\left(\mathrm{0}\le n\le N\right)$ are calculated by forward discrete cosine transform as described in Appendix B.

3. The value of G calculated from $T_{n,m=\mathrm{0}}^{\mathrm{c}\left(j\right)}$ using Eq. (G1) is considered the latitudinal area weight w(θ_j) at colatitude θ_j.

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

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

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

Here we describe the derivation of Eq. (63d) for odd m (≥3) from Eq. (62a). Equations (63b) and (63c) can also be derived similarly. By using Eqs. (52), (57), (11), (A2b), (A2c), and the same equations as Eq. (53), except that $u_{n,m}^{\mathrm{c}}$, ${\mathit{\chi }}_{n,m}^{\mathrm{s}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{c}}$ are replaced with $u_{n,m}^{\mathrm{s}}$, $-{\mathit{\chi }}_{n,m}^{\mathrm{c}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$, respectively, and the same equations as Eq. (53), except that $u_{n,m}^{\mathrm{c}}$, ${\mathit{\chi }}_{n,m}^{\mathrm{s}}$, and ${\mathit{\psi }}_{n,m}^{\mathrm{c}}$ are replaced with $v_{n,m}^{\mathrm{c}}$, ${\mathit{\psi }}_{n,m}^{\mathrm{s}}$, and $-{\mathit{\chi }}_{n,m}^{\mathrm{c}}$, respectively, Eq. (62a) is converted into the following equation for odd m≥3:

$$\begin{array}{}\text{(H1)}& \begin{array}{rl}& \underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\left\{-m{\displaystyle \frac{-\cos \left(n+\mathrm{1}\right)\mathit{\theta }+\cos \left(n-\mathrm{1}\right)\mathit{\theta }}{\mathrm{2}}}\left[{\displaystyle \frac{-m{\mathit{\chi }}_{\mathrm{1},m}^{\mathrm{c}}}{\mathrm{2}a}}-{\stackrel{\mathrm{̃}}{u}}_{\mathrm{0},m}^{\mathrm{s}}\right.\right.\\ & +\left({\displaystyle \frac{-m{\mathit{\chi }}_{\mathrm{2},m}^{\mathrm{c}}}{\mathrm{2}a}}+{\displaystyle \frac{\mathrm{3}{\mathit{\psi }}_{\mathrm{1},m}^{\mathrm{s}}-{\mathit{\psi }}_{\mathrm{3},m}^{\mathrm{s}}}{\mathrm{4}a}}-{\stackrel{\mathrm{̃}}{u}}_{\mathrm{1},m}^{\mathrm{s}}\right)\cos \mathit{\theta }\\ & +\sum _{n^{\prime }=\mathrm{2}}^N\left({\displaystyle \frac{m\left({\mathit{\chi }}_{n^{\prime }-\mathrm{1},m}^{\mathrm{c}}-{\mathit{\chi }}_{n^{\prime }+\mathrm{1},m}^{\mathrm{c}}\right)}{\mathrm{2}a}}\right.\\ & \left.\left.+{\displaystyle \frac{n\left(-{\mathit{\psi }}_{n^{\prime }-\mathrm{2},m}^{\mathrm{s}}+\mathrm{2}{\mathit{\psi }}_{n^{\prime },m}^{\mathrm{s}}-{\mathit{\psi }}_{n^{\prime }+\mathrm{2},m}^{\mathrm{s}}\right)}{\mathrm{4}a}}-{\stackrel{\mathrm{̃}}{u}}_{n^{\prime },m}^{\mathrm{s}}\right)\cos n^{\prime }\mathit{\theta }\right]\\ & -\left(-{\displaystyle \frac{n-\mathrm{2}}{\mathrm{4}}}\cos \left(n-\mathrm{2}\right)\mathit{\theta }+{\displaystyle \frac{\mathrm{2}n}{\mathrm{4}}}\cos n\mathit{\theta }-{\displaystyle \frac{n+\mathrm{2}}{\mathrm{4}}}\cos \left(n+\mathrm{2}\right)\mathit{\theta }\right)\\ & \left[{\displaystyle \frac{m{\mathit{\psi }}_{\mathrm{1},m}^{\mathrm{s}}}{\mathrm{2}a}}-{\stackrel{\mathrm{̃}}{v}}_{\mathrm{0},m}^{\mathrm{c}}+\left({\displaystyle \frac{m{\mathit{\psi }}_{\mathrm{2},m}^{\mathrm{s}}}{\mathrm{2}a}}+{\displaystyle \frac{-\mathrm{3}{\mathit{\chi }}_{\mathrm{1},m}^{\mathrm{c}}+{\mathit{\chi }}_{\mathrm{3},m}^{\mathrm{c}}}{\mathrm{4}a}}-{\stackrel{\mathrm{̃}}{v}}_{\mathrm{1},m}^{\mathrm{c}}\right)\cos \mathit{\theta }\right.\\ & +\sum _{n^{\prime }=\mathrm{2}}^N\left({\displaystyle \frac{m\left(-{\mathit{\psi }}_{n^{\prime }-\mathrm{1},m}^{\mathrm{s}}+{\mathit{\psi }}_{n^{\prime }+\mathrm{1},m}^{\mathrm{s}}\right)}{\mathrm{2}a}}\right.\\ & \left.\left.\left.+{\displaystyle \frac{n^{\prime }\left({\mathit{\chi }}_{n^{\prime }-\mathrm{2},m}^{\mathrm{c}}-\mathrm{2}{\mathit{\chi }}_{n^{\prime },m}^{\mathrm{c}}+{\mathit{\chi }}_{n^{\prime }+\mathrm{2},m}^{\mathrm{c}}\right)}{\mathrm{4}a}}-{\stackrel{\mathrm{̃}}{v}}_{n^{\prime },m}^{\mathrm{c}}\right)\cos n^{\prime }\mathit{\theta }\right]\right\}\\ & \mathrm{d}\mathit{\theta }=\mathrm{0}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(H2)}& \begin{array}{rl}& \underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\left\{{\displaystyle \frac{m}{\mathrm{2}}}\cos \left(n+\mathrm{1}\right)\mathit{\theta }\left({\displaystyle \frac{m\left({\mathit{\chi }}_{n,m}^{\mathrm{c}}-{\mathit{\chi }}_{n+\mathrm{2},m}^{\mathrm{c}}\right)}{\mathrm{2}a}}\right.\right.\\ & \left.+{\displaystyle \frac{\left(n+\mathrm{1}\right)\left(-{\mathit{\psi }}_{n-\mathrm{1},m}^{\mathrm{s}}+\mathrm{2}{\mathit{\psi }}_{n+\mathrm{1},m}^{\mathrm{s}}-{\mathit{\psi }}_{n+\mathrm{3},m}^{\mathrm{s}}\right)}{\mathrm{4}a}}-{\stackrel{\mathrm{̃}}{u}}_{n+\mathrm{1},m}^{\mathrm{s}}\right)\cos \left(n+\mathrm{1}\right)\mathit{\theta }\\ & -{\displaystyle \frac{m}{\mathrm{2}}}\cos \left(n-\mathrm{1}\right)\mathit{\theta }\left({\displaystyle \frac{m\left({\mathit{\chi }}_{n-\mathrm{2},m}^{\mathrm{c}}-{\mathit{\chi }}_{n,m}^{\mathrm{c}}\right)}{\mathrm{2}a}}\right.\\ & \left.+{\displaystyle \frac{\left(n-\mathrm{1}\right)\left(-{\mathit{\psi }}_{n-\mathrm{3},m}^{\mathrm{s}}+\mathrm{2}{\mathit{\psi }}_{n-\mathrm{1},m}^{\mathrm{s}}-{\mathit{\psi }}_{n+\mathrm{1},m}^{\mathrm{s}}\right)}{\mathrm{4}a}}-{\stackrel{\mathrm{̃}}{u}}_{n-\mathrm{1},m}^{\mathrm{s}}\right)\cos \left(n-\mathrm{1}\right)\mathit{\theta }\\ & +{\displaystyle \frac{n-\mathrm{2}}{\mathrm{4}}}\cos \left(n-\mathrm{2}\right)\mathit{\theta }\left({\displaystyle \frac{m\left(-{\mathit{\psi }}_{n-\mathrm{3},m}^{\mathrm{s}}+{\mathit{\psi }}_{n-\mathrm{1},m}^{\mathrm{s}}\right)}{\mathrm{2}a}}\right.\\ & \left.+{\displaystyle \frac{\left(n-\mathrm{2}\right)\left({\mathit{\chi }}_{n-\mathrm{4},m}^{\mathrm{c}}-\mathrm{2}{\mathit{\chi }}_{n-\mathrm{2},m}^{\mathrm{c}}+{\mathit{\chi }}_{n,m}^{\mathrm{c}}\right)}{\mathrm{4}a}}-{\stackrel{\mathrm{̃}}{v}}_{n-\mathrm{2},m}^{\mathrm{c}}\right)\cos \left(n-\mathrm{2}\right)\mathit{\theta }\\ & -{\displaystyle \frac{\mathrm{2}n}{\mathrm{4}}}\cos n\mathit{\theta }\left({\displaystyle \frac{m\left(-{\mathit{\psi }}_{n-\mathrm{1},m}^{\mathrm{s}}+{\mathit{\psi }}_{n+\mathrm{1},m}^{\mathrm{s}}\right)}{\mathrm{2}a}}\right.\\ & \left.+{\displaystyle \frac{n\left({\mathit{\chi }}_{n-\mathrm{2},m}^{\mathrm{c}}-\mathrm{2}{\mathit{\chi }}_{n,m}^{\mathrm{c}}+{\mathit{\chi }}_{n+\mathrm{2},m}^{\mathrm{c}}\right)}{\mathrm{4}a}}-{\stackrel{\mathrm{̃}}{v}}_{n,m}^{\mathrm{c}}\right)\cos n\mathit{\theta }\\ & +{\displaystyle \frac{n+\mathrm{2}}{\mathrm{4}}}\cos \left(n+\mathrm{2}\right)\mathit{\theta }\left({\displaystyle \frac{m\left(-{\mathit{\psi }}_{n+\mathrm{1},m}^{\mathrm{s}}+{\mathit{\psi }}_{n+\mathrm{3},m}^{\mathrm{s}}\right)}{\mathrm{2}a}}\right.\\ & \left.+{\displaystyle \frac{\left(n+\mathrm{2}\right)\left({\mathit{\chi }}_{n,m}^{\mathrm{c}}-\mathrm{2}{\mathit{\chi }}_{n+\mathrm{2},m}^{\mathrm{c}}+{\mathit{\chi }}_{n+\mathrm{4},m}^{\mathrm{c}}\right)}{\mathrm{4}a}}-{\stackrel{\mathrm{̃}}{v}}_{n+\mathrm{2},m}^{\mathrm{c}}\right)\\ & \left.\cos \left(n+\mathrm{2}\right)\mathit{\theta }\right\}\mathrm{d}\mathit{\theta }=\mathrm{0}.\end{array}\end{array}$$

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

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

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

$$\begin{array}{}\text{(I1)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\left(\mathit{v}+\mathrm{2}\mathbf{\Omega }\times \mathit{r}\right)-{\left(\mathit{v}+\mathrm{2}\mathbf{\Omega }\times \mathit{r}\right)}_{\text{D}}^{\mathrm{0}}}{\mathrm{\Delta }t}}=-{\displaystyle \frac{g\left(\mathrm{\nabla }{h}_{\text{D}}^{\left(+\right)}+\mathrm{\nabla }{h}^{\mathrm{0}}\right)}{\mathrm{2}}}\\ & +{\mathit{\beta }}_{\mathit{v}}{\displaystyle \frac{g\left(\mathrm{\nabla }{h}_{\text{D}}^{\left(+\right)}+\mathrm{\nabla }{h}^{\mathrm{0}}\right)}{\mathrm{2}}}-{\mathit{\beta }}_{\mathit{v}}{\displaystyle \frac{g\left(\mathrm{\nabla }{h}_{\text{D}}^{\mathrm{0}}+\mathrm{\nabla }{h}^{+}\right)}{\mathrm{2}}},\end{array}\text{(I2)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{h^{+}-h_{\text{D}}^{\mathrm{0}}}{\mathrm{\Delta }t}}=-{\displaystyle \frac{{\left[\left(h-h_{\mathrm{s}}\right)D\right]}_{\text{D}}^{\left(+\right)}+{\left[\left(h-h_{\mathrm{s}}\right)D\right]}^{\mathrm{0}}}{\mathrm{2}}}\\ & +{\displaystyle \frac{{\left[\mathit{v}\cdot \mathrm{\nabla }{h}_{\mathrm{s}}\right]}_{\text{D}}^{\left(+\right)}+{\left[\mathit{v}\cdot \mathrm{\nabla }{h}_{\mathrm{s}}\right]}^{\mathrm{0}}}{\mathrm{2}}}\\ & +{\mathit{\beta }}_h{\displaystyle \frac{{\left[\overline{h}D\right]}_{\text{D}}^{\left(+\right)}+{\left[\overline{h}D\right]}^{\mathrm{0}}}{\mathrm{2}}}-{\mathit{\beta }}_h{\displaystyle \frac{{\left[\overline{h}D\right]}_{\text{D}}^{\mathrm{0}}+{\left[\overline{h}D\right]}^{+}}{\mathrm{2}}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(I3)}& D\equiv \mathrm{\nabla }\cdot \mathit{v}={\displaystyle \frac{\mathrm{1}}{a}}\left[{\displaystyle \frac{\mathrm{1}}{\cos \mathit{\varphi }}}{\displaystyle \frac{\partial u}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\mathrm{1}}{\cos \mathit{\varphi }}}{\displaystyle \frac{\partial v\cos \mathit{\varphi }}{\partial \mathit{\varphi }}}\right],\end{array}$$

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, $h^{\left(+\right)}=\mathrm{2}{h}^{\mathrm{0}}-h^{-}$; the subscript D means the departure point; the absence of the subscript D means the arrival point; $\overline{h}$ is a constant value of height for semi-implicit linear terms; β_v and β_h are second-order decentering parameters (Yukimoto et al., 2011). Using β_v and β_h larger than 1.0 (e.g., 1.2) increases the effect of the semi-implicit scheme improving computational stability, but ${\mathit{\beta }}_{\mathit{v}}={\mathit{\beta }}_h=\mathrm{1.0}$ is used here because $\overline{h}$ larger than h is enough for stable calculations in the shallow-water model. The departure point x_D 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 x_D 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, x_D is calculated using

$$\begin{array}{}\text{(I4)}& {\mathit{x}}_{\text{D}}\text{=}\mathit{x}-{\displaystyle \frac{{\mathit{v}}_{\text{D}}^{\left(+\right)}+{\mathit{v}}^{\mathrm{0}}}{\mathrm{2}}}\mathrm{\Delta }t.\end{array}$$

However, when Δt is longer than 30 min, using ${\mathit{v}}_{\text{D}}^{\left(+\right)}$ extrapolated in time to calculate x_D causes numerical instability in our experiments. To avoid instability when Δt is 1 h, here we use

$$\begin{array}{}\text{(I5a)}& {\displaystyle }{\mathit{x}}_{\text{D}}\text{=}\mathit{x}-{\displaystyle \frac{{\mathit{v}}_{\text{D}}^{\mathrm{0}}+{\mathit{v}}^{\prime +}}{\mathrm{2}}}\mathrm{\Delta }t,\text{(I5b)}& {\displaystyle }{\mathit{v}}^{\prime +}\equiv {\mathit{v}}_{\text{D}}^{\mathrm{0}}+{\left(\text{2}\mathbf{\Omega }\times \mathit{r}\right)}_{\text{D}}-\mathrm{2}\mathbf{\Omega }\times \mathit{r}-{\displaystyle \frac{g\left(\mathrm{\nabla }{h}_{\text{D}}^{\left(+\right)}+\mathrm{\nabla }{h}^{\mathrm{0}}\right)}{\mathrm{2}}}\mathrm{\Delta }t,\end{array}$$

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

$$\begin{array}{}\text{(I6)}& \begin{array}{rl}{\mathit{x}}_{\text{D}}& =\mathit{x}-\mathrm{\Delta }t\left[{\left({\mathit{v}}^{\mathrm{0}}+\mathbf{\Omega }\times \mathit{r}-{\displaystyle \frac{g\mathrm{\Delta }t\mathrm{\nabla }{h}^{\left(+\right)}}{\mathrm{4}}}\right)}_{\text{D}}\right.\\ & \left.-\mathbf{\Omega }\times \mathit{r}-{\displaystyle \frac{g\mathrm{\Delta }t\mathrm{\nabla }{h}^{\mathrm{0}}}{\mathrm{4}}}\right].\end{array}\end{array}$$

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

Equations (I1) and (I2) are converted into

$$\begin{array}{}\text{(I7a)}& {\displaystyle }{\mathit{v}}^{+}+{\displaystyle \frac{{\mathit{\beta }}_{\mathit{v}}\mathrm{\Delta }t}{\mathrm{2}}}g\mathrm{\nabla }{h}^{+}={\mathit{R}}_{\mathit{v}},\text{(I7b)}& {\displaystyle }\begin{array}{rl}{\mathit{R}}_{\mathit{v}}& \equiv {\left[{\mathit{v}}^{\mathrm{0}}+\mathrm{2}\mathbf{\Omega }\times \mathit{r}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}g\left(\mathrm{\nabla }{h}^{\left(+\right)}-{\mathit{\beta }}_{\mathit{v}}\mathrm{\nabla }{h}^{\left(+\right)}+{\mathit{\beta }}_{\mathit{v}}\mathrm{\nabla }{h}^{\mathrm{0}}\right)\right]}_{\text{D}}\\ & -\text{2}\mathbf{\Omega }\times \mathit{r}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}g\left(\mathrm{\nabla }{h}^{\mathrm{0}}-{\mathit{\beta }}_{\mathit{v}}\mathrm{\nabla }{h}^{\mathrm{0}}\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(I8a)}& {\displaystyle }{h}^{+}+{\displaystyle \frac{{\mathit{\beta }}_h\mathrm{\Delta }t}{\mathrm{2}}}\overline{h}{D}^{+}=R_h,\text{(I8b)}& {\displaystyle }\begin{array}{rl}& R_h\equiv \left\{h^{\mathrm{0}}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}\left[-\left(h-h_{\mathrm{s}}\right)D^{\left(+\right)}+\mathit{v}\cdot \mathrm{\nabla }{h}_{\mathrm{s}}^{\left(+\right)}\right.\right.\\ & {\left.\left.+{\mathit{\beta }}_h\overline{h}{D}^{\left(+\right)}-{\mathit{\beta }}_h\overline{h}{D}^{\mathrm{0}}\right]\right\}}_{\text{D}}\\ & +{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}\left[-\left(h-h_{\mathrm{s}}\right)D^{\mathrm{0}}+\mathit{v}\cdot \mathrm{\nabla }{h}_{\mathrm{s}}^{\mathrm{0}}+{\mathit{\beta }}_h\overline{h}{D}^{\mathrm{0}}\right].\end{array}\end{array}$$

In Eqs. (I7b) and (I8b), the values at x_D 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

$$\begin{array}{}\text{(I9)}& {\displaystyle }& {\displaystyle }{D}^{+}+{\displaystyle \frac{{\mathit{\beta }}_{\mathit{v}}\mathrm{\Delta }t}{\mathrm{2}}}g{\mathrm{\nabla }}^{\mathrm{2}}{h}^{+}=R_D,\text{(I10)}& {\displaystyle }& {\displaystyle }{\mathit{\zeta }}^{+}=R_{\mathit{\zeta }},\end{array}$$

where

$$\begin{array}{}\text{(I11)}& \mathit{\zeta }\equiv \mathit{k}\cdot \mathrm{\nabla }\times \mathit{v}\text{=}{\displaystyle \frac{\mathrm{1}}{a}}\left[{\displaystyle \frac{\mathrm{1}}{\cos \mathit{\varphi }}}{\displaystyle \frac{\partial v}{\partial \mathit{\lambda }}}-{\displaystyle \frac{\mathrm{1}}{\cos \mathit{\varphi }}}{\displaystyle \frac{\partial u\cos \mathit{\varphi }}{\partial \mathit{\varphi }}}\right],\end{array}$$

is vorticity, $\mathit{k}\equiv \mathit{r}/\left|\mathit{r}\right|$ is the vertical unit vector, $R_D\equiv \mathrm{\nabla }\cdot {\mathit{R}}_{\mathit{v}}$, and $R_{\mathit{\zeta }}\equiv \mathit{k}\cdot \mathrm{\nabla }\times {\mathit{R}}_{\mathit{v}}$.

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

1. The scalar variable R_h is transformed from grid space to spectral space using Eqs. (9), (23), (30), and (31). The components of the vector variable ${\mathit{R}}_{\mathit{v}}=\left(R_u,R_v\right)$ 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 R_v, respectively.

2. R_D and R_ζ are calculated by

   $$\begin{array}{}\text{(I12)}& {\displaystyle }{R}_D={\mathrm{\nabla }}^{\mathrm{2}}{R}_{\mathit{\chi }},\text{(I13)}& {\displaystyle }{R}_{\mathit{\zeta }}={\mathrm{\nabla }}^{\mathrm{2}}{R}_{\mathit{\psi }},\end{array}$$

   using Eqs. (75) and (77). ζ⁺ is obtained from R_ζ using Eq. (I10).

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

   $$\begin{array}{}\text{(I14)}& \begin{array}{rl}{D}^{+}& -{\left({\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}\right)}^{\mathrm{2}}{\mathit{\beta }}_{\mathit{v}}{\mathit{\beta }}_{h}g\overline{h}{\mathrm{\nabla }}^{\mathrm{2}}{D}^{+}\\ & ={\mathrm{\nabla }}^{\mathrm{2}}\left(R_{\mathit{\chi }}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\mathit{\beta }}_{\mathit{v}}g{R}_h\right).\end{array}\end{array}$$

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

4. h⁺ is calculated from D⁺ and R_h using Eq. (I8).

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

   $$\begin{array}{}\text{(I15)}& {\displaystyle }{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\chi }}^{+}=D^{+},\text{(I16)}& {\displaystyle }{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\psi }}^{+}={\mathit{\zeta }}^{+},\end{array}$$

   using Eqs. (75) and (78).

6. ${\mathit{v}}^{+}=\left(u^{+},v^{+}\right)$ is calculated from χ⁺ and ψ⁺ using Eq. (53) for $u_{n,m}^{\mathrm{c}}$ and the similar equations for $u_{n,m}^{\mathrm{s}}$, $v_{n,m}^{\mathrm{c}}$, and $v_{n,m}^{\mathrm{s}}$.

7. u⁺, v⁺, h⁺, D⁺, and ∇h⁺ in spectral space are transformed to grid space. h⁺ and D⁺ are transformed meridionally using Eqs. (14) and (13). u⁺ and v⁺ are transformed meridionally using Eq. (52). $\mathrm{\nabla }{h}^{+}=\left(h_{\mathit{\lambda }}^{+},h_{\mathit{\theta }}^{+}\right)$ is transformed meridionally using Eqs. (18)–(21). $h_{\mathit{\lambda }}^{+}$ can also be calculated from $h_m^{+\mathrm{c},N}\left({\mathit{\theta }}_j\right)$ and $h_m^{+\mathrm{s},N}\left({\mathit{\theta }}_j\right)$ at the latitudinal grid points using Eq. (16) and additionally using Eq. (22) at the poles when using Grid [1], which is more efficient than using Eqs. (18) and (19) because the meridional inverse discrete cosine and sine transforms of $h_{\mathit{\lambda }}^{+}$ become unnecessary.