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

The computational cost of a spectral model using spherical harmonics (SH) increases significantly at high resolution because the transform method with SH requires O N operations, where N is the truncation wavenumber. One way to solve this problem is to use double Fourier series (DFS) instead of SH, which requires 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 10 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. In the semi-implicit semi-Lagrangian shallow water model using the new DFS method, the Williamson test cases 2 and 5 and the Galewsky test case give stable results without the appearance of high-wavenumber noise near the poles, even without using horizontal diffusion and a zonal Fourier filter. The new DFS model is faster than the SH model, especially at high resolutions, and gives almost the same results. 15


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., 2011Yukimoto et al., , 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 accuracy in horizontal derivatives is 20 good, 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 timesteps 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 operations and O memory 25 usage, where is the truncation wavenumber. One way to reduce the operation count and the memory usage at high resolutions with large is to use the fast Legendre transform (Suda, 2005;Tygert, 2008;Wedi et al., 2013), which requires only O log operations, although the accuracy is compromised to reduce the operation count. 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 used to improve the Legendre transform is on-the-fly computation of the associated Legendre functions (Schaeffer, 2013;Ishioka, 2018), which still requires O operations but requires only O memory usage. This small memory usage also contributes to speeding up calculations by taking advantage of the cache memory.
Another way to reduce the operation count and the memory usage in the global spectral model is to use double Fourier series (DFS) as basis functions. In the DFS model, the fast Fourier transform (FFT; Cooley and Tukey, 1965;Swarztrauber, 1982) 5 is used not only in the longitudinal (zonal) direction but also in the latitudinal (meridional) direction. The FFT requires only O log operations and O or O memory usage, and it is much faster than the fast Legendre transform.
In DFS models (and also in SH models), the scalar variable , is zonally expanded as , ≅ , 1 where is longitude, is colatitude, and is the zonal truncation wavenumber. Several methods have been proposed for 10 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 with 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. 15 Orszag (1974) and Boyd (1978)  where is the meridional truncation wavenumber. The coefficients , for odd are calculated from the forward Fourier cosine transform of sin ⁄ . Orszag (1974) imposed the following conditions at the poles: 20 0 0 and 0 for | | 2, 3 which can be expressed in terms of the expansion coefficients , as , 0 and , 0 for | | 2. 4 Satisfying the above conditions ensures that the scalar variable , and its gradient ∇ are continuous at the poles. In Orszag (1974), only , and , were modified to satisfy Eq. 4 , but this is not the best way to satisfy the same conditions 25 as Eq. 3 or Eq. 4 , as will be shown in Sect. 3. Yee (1981) and Layton and Spotz (2003) expanded as https://doi.org /10.5194/gmd-2021-168 Preprint. Discussion started: 8 July 2021 c Author(s) 2021. CC BY 4.0 License.

6
The meridional basis functions sin sin for even 0 are different from Eq. 5 . The coefficients , for even 0 are calculated by forward Fourier sine transform of sin ⁄ . The basis functions in Eq. 6 automatically satisfy the same conditions at the poles as Eq. 3 for even , and guarantee the continuity of the scalar variable at the poles, which is 10 an advantage compared with the basis functions in Eq. 5 . However, Eq. 6 does not automatically satisfy the conditions in Eq. 3 for odd , and does not guarantee the continuity of ∇ at the poles. The shallow water model and the vorticity equation model using a semi-implicit Eulerian scheme ran stably using high-order horizontal diffusion with O 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 15 et al., 2013). However, the computational results of these models appear to be a little different from (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 ≅ 2 3 ⁄ or 2 ⁄ for even 0 (See Sect. 3), and the use of strong high-order horizontal diffusion to smooth out the oscillation, where 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 20 hydrostatic and nonhydrostatic atmospheric models using the DFS basis functions of Cheong in Eq. 6 . These models used meridional truncation with ≅ , and sin and sin (instead of sin ⁄ and sin ⁄ ) were transformed from grid space to spectral space, where is the zonal wind and is the meridional wind. These models used the same horizontal diffusion as the SH models, and did not require the strong high-order horizontal diffusion. The results of these models were very similar to those of the SH models. However, we found the following two problems in these models: 25 https://doi.org /10.5194/gmd-2021-168 Preprint.  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-5 squares method (or the Galerkin method) to minimize the error due to the meridional wavenumber truncation.
2 New DFS basis functions that automatically satisfy the pole conditions in Eq. 3 , which guarantee 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 and the shallow water equations. 10 Section 2 describes the details of the new DFS method using the new DFS expansion method and the new DFS basis functions. Section 3 examines the error due to the wavenumber truncation in the new DFS method, Orszag's DFS method, and Cheong's DFS method. Section 4 describes how to integrate the semi-implicit semi-Lagrangian shallow water model using the new DFS method. Section 5 compares the results of the model using the new DFS method with those using the old DFS method of Yoshimura and Matsumura (2005), and with those using the SH method. Section 6 presents conclusions and perspectives. 15 2 Improved double Fourier series on the sphere

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.

. The scalar variable
, is expanded zonally as  Dividing by sin reduces the numerical stability of the model more significantly than dividing by sin .
Here we propose a new method to calculate expansion coefficients using the least-squares method to minimize the error due to the meridional wavenumber truncation. This method also avoids dividing by sin or sin before the forward 10 cosine or sine transforms. The coefficients , and , in Eq. 8 are calculated as follows. First, and in Eq.   24d with the exception of the following underlined values: 20 From Eq. 24d , two linear simultaneous equations with penta-diagonal matrices, To solve , we solve with forward substitution first, and then solve with backward substitution. There are also other methods to solve Eq. 25 . For example, the method using LU decomposition considering penta-diagonal matrices as 2 2 block tri-diagonal matrices makes SIMD operations more effective. The method using cyclic reduction for block tri-diagonal matrices (e.g., Gander and Golub, 1997

Relation between the least-squares method and Galerkin method for a scalar variable
Here we discuss the relation between the least-squares method described above and the Galerkin method when calculating the expansion coefficients of a scalar variable.

Comparison of new DFS with SH
Here we compare the new DFS method with the SH method to see the difference between them. In the SH method, and in Eq. 7 are expanded with the associated Legendre functions , as 5 where 0. The functions , satisfy the following orthogonality relations for each : , , sin 1 for , 0 for . 31 By the modified Robert expansion (Merilees, 1973a;Orszag, 1974) 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 Sect. 2.11) 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 , , and , , in Eq. 30 are usually calculated from the grid-point values of and , respectively, by using Gaussian quadrature or Clenshaw-Curtis quadrature (e.g., Hotta and Ujiie, 2018). They can also 5 be calculated from , and , in Eq. 19 instead of and at the grid points as follows (e.g., Sneeuw and Bun, 1996):  20 (and Eq. 27 ), which is another difference between the SH and the new DFS methods. In the DFS method, the constant 20 latitudinal weight is used in Eq. 20 , although the latitudinal area weight described below in Appendix B is usually used as the latitudinal weight at the grid points.
When calculating the coefficients , (and , ) in Eq. 8 , we can also consider the least-squares method, not using in Eq. 20 but using with latitudinal weight sin like Eq. 35 . However, minimizing derives the simultaneous equations for calculating , with dense matrices, which leads to O operations. When using , the simultaneous equations with 25 penta-diagonal matrices require only O operations. Therefore, we choose to use instead of ′ .

Application of the new basis functions to a wind vector
The velocity potential and the stream function can be converted into the wind vector components and using the equations Thus, it is guaranteed that the wind vector , , , in Eqs. 45 and 46 is continuous at the poles.

New method to calculate expansion coefficients for a wind vector
We propose a new method that calculates the expansion coefficients , , , , , and , using the least-squares 5 odd . 60c From Eq. 60 , and from the same equations as Eqs. 49b, c, d , except that , , , , and , are replaced with , , , , and , , respectively, and the same equations as Eqs. 49a, b, c, d , except that , , , , and , are replaced with , , , , and , , respectively, we derive the following equations for , and , . For 0, 20 where and are nine-diagonal matrices. We also derive two similar linear simultaneous equations with penta-diagonal matrices for 1 and each even 2. The simultaneous equations with nine-diagonal or penta-diagonal matrices can be solved in a similar way to Eq. 25 , and the expansion coefficients , and , in Eq. 62 can be solved efficiently. From the same equations as Eqs. 61b, c, d , except that , , , and are replaced with , , , and , respectively, 5 and the same equations as Eqs. 61b, c, d , except that , , , and are replaced with , , , and , respectively, two similar linear simultaneous equations with nine-diagonal matrices for each 3 and two linear simultaneous equations with penta-diagonal matrices for 1 and each even 2 are also derived. Thus, the expansion coefficients , , , , , , and , are obtained from , , , , , , and , using Eqs. 61a, b, c, d and the similar equations. The expansion coefficients , , , , , , and , are obtained from , , , , , , and , using Eq. 49 for 10 , and the similar equations for , , , , and , .

Relation between the least-squares method and the Galerkin method for the wind vector
Here we discuss the relation between the least-squares method described above and the Galerkin method when calculating the expansion coefficients related to the wind vector.    , 1989;Temperton, 1991;Swarztrauber, 1993), where the SH expansion coefficients of the divergence ∇ and the 10 vorticity ∇ are calculated from the grid-point values of and using the Galerkin spectral method with the orthogonal vector SH basis functions.

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  Merilees (1973b), Orszag (1974), Cheong (2000aCheong ( , 2000b, and Yoshimura and Matsumura (2005). Grid [1] was used, for example, in the DFS expansion in Yee (1981)

Zonal Fourier filter
In regular longitude-latitude grids, 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. In this study, we set the largest zonal wavenumber at each latitude as 5 min , sin , 76 where we use the value 20 to make the resolution similar to that in the reduced grid of Miyamoto (2006). The values of and in Eq. 7 are set to zero for during the spectral transform. The use of a reduced grid (Hortal and Simmons, 1991;Juang, 2004;Miyamoto, 2006) has a similar effect to the zonal Fourier filter.

Laplacian operator and Poisson equation 10
The calculation of the Laplacian operator and the Poisson equation in the new DFS method is described in this section. In the equation where , and , are the vectors whose components are , ( is even) and , ( is odd), respectively, and 15 , _ and , _ are the vectors whose components are , ( is even) and , ( is odd), respectively; which can be solved efficiently as in Eq. 25 .
By using Eq. 87b instead of Eq. 87a , we obtain the equations to calculate , from , , which are the same as Eqs. 20 89 to 91 , except that the superscript c is replaced with the superscript s.
We have verified that all the eigenvalues of the matrices

Horizontal diffusion
The horizontal diffusion is calculated in the same way as in Cheong et al. (2004)  ; for example, and , where the superscript c indicates the zonal cosine component, and the superscript s indicates the zonal sine component.

The error due to meridional wavenumber truncation in DFS expansion methods 5
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), only , and , are modified to satisfy Eq. 4 equivalent to Eq. 3 . The DFS meridional basis functions of Cheong in Eq. 6 automatically satisfy the pole conditions in Eq. 3 for even , but not for odd . The new DFS meridional basis functions in Eq. 8 automatically satisfy the condition in Eq. 3 for both even and odd . We compare the error due to the wavenumber truncation among these DFS methods. 10 In the method of Orszag using Eq. 2 , a very large error occurs especially for odd | | 3 (Fig. 2) when , and , are modified to satisfy the pole conditions in Eq. 4 . Dividing by sin before the forward Fourier cosine transform for odd | | 3 also contributes to the large error.
In the method of Cheong using Eq. 6 , large high wavenumber oscillations appear for even 0 in Fig. 2 In the new DFS method described in Sect. 2, the usual small oscillations from the Gibbs phenomenon appear in Fig. 2, but the error is small because the expansion coefficients are calculated using the least-squares method (or the Galerkin method) to 30 https://doi.org /10.5194/gmd-2021-168 Preprint.  where is time, is the horizontal wind vector, ℎ is the height, ℎ is the surface height, is the acceleration due to gravity, is the 3-dimensional angular velocity of the earth's rotation, and the subscript H indicates the horizontal component. Equation   104 is transformed for the advective treatment of the Coriolis term (Temperton, 1997) where is the 3-dimensional position vector from the Earth's center. Equation 105 is transformed for the spatially averaged 15 Eulerian treatment of mountains (Ritchie and Tanguay, 1996) into ℎ ℎ ℎ ∇ • • ∇ℎ . 107

Time integration method
A two-time-level semi-implicit semi-Lagrangian scheme (e.g., Temperton et al., 2001)  is horizontal divergence; ∆ is a timestep; the superscripts , 0, and mean past time ∆ , present time , and future time ∆ , respectively, and the superscript means future time ∆ extrapolated in time, for example, ℎ 2ℎ ℎ ; the subscript D means the departure point, and the absence of the subscript D means the arrival point; ℎ is a constant value of height for semi-implicit linear terms; and are second-order decentering parameters (Yukimoto et al., 2011). 5 Using and larger than 1.0 (e.g., 1.2) increases the effect of the semi-implicit scheme improving computational stability, but 1.0 is used here because ℎ larger than ℎ is enough for stable calculations in the shallow water model. The  This method using a provisional future value to calculate is similar to the method in Gospodinov et al., (2001). Since the value with the subscript D depends on , is calculated iteratively from Eq. 113 (e.g., Ritchie, 1995;Temperton et al., 2001). Since is not generally on the grid point, the value at is calculated by spatial interpolation from nearby grid points.
In the right-hand side of Eq. 113 , the value at with the subscript D is calculated by third-order Lagrange interpolation. We calculate ℎ and using the spectral transform method and the Galerkin method with the new DFS method as follows. 7. , , ℎ , , and ∇ℎ in spectral space are transformed to grid space. ℎ and are transformed meridionally using Eqs. 9 and 10 . and are transformed meridionally using Eq. 48 . ∇ℎ ℎ , ℎ is transformed meridionally using Eqs. 14 to 17 . ℎ can also be calculated from ℎ , and ℎ , at the latitudinal grid points using Eq.
12 , and additionally using Eq. 18 at the poles when using Grid [1], which is more efficient than using Eqs. 14) and 15 because the meridional inverse discrete cosine and sine transforms of ℎ become unnecessary. 5

Models
We ran Williamson test cases 2 and 5 (Williamson et al., 1992) and the Galewsky test case (Galewsky et al., 2004) in the semi-implicit semi-Lagrangian shallow water model using the new improved DFS method described in Sect. 2 (hereafter the new DFS model). We also ran the same test cases in the semi-implicit semi-Lagrangian shallow water model using the DFS 10 method of Yoshimura and Matsumura (2005) with the basis functions of Cheong (2000aCheong ( , 2000b  The zonal Fourier transforms in all of the models and the meridional Fourier cosine and sine transforms in the DFS models are calculated using the Netlib BIHAR library, which is 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. 25

Williamson test case 2
The Williamson test case 2 simulates a steady state non-linear 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 become  In the old DFS model at high resolution with 1920  960 grid points, the high wavenumber noise is not seen in Fig. 5. The higher the resolution, the smaller the high wavenumber noise becomes. Figure 6 shows the kinetic energy spectra of the horizontal winds (Lambert, 1984) after a 15-day 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 firstly 20 calculating the Gaussian grid-point values from the DFS coefficients using Eq. 8 for the new DFS method and Eq. 6 for the old DFS method, and secondly calculating the SH expansion coefficients from the Gaussian grid-point values by using a forward Legendre transform. In the old DFS model with 128  64 grid points, 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. 5. In the old DFS model with 1920  960 grid points, the high wavenumber components are a little larger than in the other models, but the differences 25 are slight. Figure 7 shows the predicted height after a 15-day integration in Williamson test case 5, which is the same as Fig. 4 except for the truncation wavenumber . In our semi-implicit semi-Lagrangian models, we usually use satisfying ≅ 1 ( is the number of latitudinal grid points), which is called linear truncation. However, here is determined to satisfy ≅ 2 1 3 ⁄ to eliminate aliasing errors with quadratic nonlinearity (Orszag, 1971), which is called quadratic truncation. occurs because of the high-wavenumber oscillations due to the quadratic wavenumber truncation for even 0 , as explained in Sect. 3. The results for the new DFS models are almost the same as for the SH model. Figure 8 shows the kinetic energy spectrum of the horizontal winds after a 15-day integration in Williamson test case 5, which is the same as Fig. 6 except for the truncation wavenumber . At the resolution 42 with 128  64 grid points, the high wavenumber components are a little larger in the SH model than in the new DFS model. At the resolution 639 with 1920  960 grid points, small 5 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 and divided by sin are transformed from grid space to spectral space (Ritchie, 1988;Temperton, 1991), which seems to be the cause of the small oscillation in the high wavenumber region. Another way to transform and from grid space to spectral space in the SH model is to use the vector harmonic transform (see Sect. 2.8), which avoids dividing and by sin and improves the stability of the model (Swarztrauber, 2004). This approach is similar to the expansion method 10 for and using the least-squares method in the new DFS method described in Sects. 2.7 and 2.8, and probably solves the problem with the high wavenumber components in the SH model. Alternatively, using and instead of and as prognostic variables may mitigate this problem.

Galewsky test case
The Galewsky test case simulates a barotropically unstable mid-latitude jet. Figure 9 shows the predicted vorticity after a 6-15 day integration in the Galewsky test case for the models at 1.3 km resolution with 30720  15360 grid points and the quadratic truncation 10239, 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 Fig. 7. Figure 10 shows the kinetic energy spectrum of horizontal winds after a 6-day 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 small oscillations appear near the truncation wavenumber in the 20 SH model. This is probably for the same reason as in Williamson test case 5 in Fig. 8. Figure 11 shows the elapsed time for the 15-

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 5 Fourier cosine or sine transform.
2. New DFS basis functions that guarantee 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 and the shallow water equations. 10 To test the new DFS method, we conducted experiments for the Williamson test cases 2 and 5, and the Galewsky test case in semi-implicit semi-Lagrangian shallow water models using the new DFS method with the three types of equally spaced latitudinal grids with or without the poles. We compared the results of the new DFS models using the new DFS method with the old DFS model using the method of Yoshimura and Matsumura (2005), and with 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 15 DFS models. This is because the new DFS expansion method with the least-squares method improves the model's stability. In the old DFS model, a truncation wavenumber lower than the number of latitudinal grid points for even 0 causes numerical instability. In the new DFS model, an arbitrary meridional wavenumber truncation 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 the quadratic truncation 20 ≅ 2 1 3 ⁄ or the cubic truncation ≅ 1 2 ⁄ is usually used in the Eulerian model and is also becoming to be used in the semi-Lagrangian model instead of the linear truncation ≅ 1 for stability and efficiency at high resolutions (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. 25 The results of the new DFS shallow water model are almost the same as the SH shallow water model. But in the SH model without horizontal diffusion, 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 and divided by sin are transformed from grid space to spectral space in the SH model. This problem with the SH model can probably be solved by using the vector harmonic transform, which is similar to the expansion method for and using the least-squares method in the new DFS 30 model. 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, which will be described in another paper after improving the nonhydrostatic 5 dynamical core as needed. license. These models utilize the Netlib BIHAR library and the ISPACK library. The Netlib BIHAR library is available at 10 https://www.netlib.org/bihar/ and is also included in the Supplement. The ISPACK library is available at https://www. . B1 The latitudinal area weight at each latitude is calculated as follows: 1. The latitudinal distribution of for each is given as 1 for 0 for 0 1 . B2 2. From , the meridional expansion coefficients , 0 are calculated by forward discrete cosine 5 transform described in Sect. 2.10.
3. The value of calculated from , using Eq. C1 is considered as the latitudinal area weight at latitude .
The latitudinal area weight is used, for example, to calculate the global mean in the grid space.