{'fraction_non_alphanumeric': 0.038931886735099494, 'fraction_numerical': 0.01606378749970772, 'mean_word_length': 4.763881401617251, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A nestable, multigrid-friendly grid on a sphere for global spectral models based on Clenshaw-Curtis quadrature

A new grid system on a sphere is proposed which allows for straightforward implementation of both spherical‐harmonics‐based spectral methods and gridpoint‐based multigrid methods. The latitudinal gridpoints in the new grid are equidistant and spectral transforms in the latitudinal direction are performed using Clenshaw–Curtis quadrature. The spectral transforms with this new grid and quadrature are shown to be exact within machine precision provided that the grid truncation is such that there are at least 2N + 1 latitudinal gridpoints for the total truncation wavenumber of N. The new grid and quadrature is implemented and tested on a shallow‐water equations model and the hydrostatic dry dynamical core of the global NWP model JMA‐GSM. The integration results obtained with the new quadrature are shown to be almost identical to those obtained with the conventional Gaussian quadrature on a Gaussian grid. Only minor code changes are required to adapt any Gaussian‐based spectral models to employ the proposed quadrature.


Introduction
Global spectral atmospheric models that are in use today almost universally adopt Gaussian quadrature in the meridional direction to perform forward (gridpoint-to-wavenumber) spherical harmonics transform. Gaussian quadrature is an optimal quadrature rule in the sense of maximising the degree of polynomials that can be integrated exactly for a given number of quadrature points (or nodes). Given J nodes, Gaussian quadrature is exact for integrand polynomials of up to as high as 2J − 1 degrees. This optimality is achieved, however, with several inconveniences (e.g. Clenshaw and Curtis , 1960): First, the nodes and weights are not given in an explicit analytic form and necessitates (some iterative) solution of an algebraic equation of high degrees. Second, the nodes do not nest, i.e., the nodes for the J-point rule do not contain any J ′ -point nodes as their subset for any J ′ < J. This second limitation imposes inflexibility to the future evolution of a global spectral dynamical core. In particular, the unnestable grid alignment makes it difficult to combine the current spectral dynamical core with a multigrid approach for solving an elliptic boundary value problem that arises from (semi-)implicit time discretisation.
As the horizontal resolution continues to increase, it becomes necessary to properly represent nature's non-hydrostatic aspects in the model. In a non-hydrostatic system, a semi-implicit time stepping results in a Helmholtz problem with spatially-variable coefficients that needs to be solved iteratively, even with horizontal spectral discretisation (e.g., Bénard , 2003). This is in contrast to the hydrostatic case where the resultant Helmholtz problem has constant coefficients and thus can be solved without iteration in the spectral space (Hoskins and Simmons , 1975).
The multigrid approach is an attractive strategy for iterative solution of the non-hydrostatic implicit Helmholtz problem. It accelerates convergence of the iterative algorithm by first solving the problem at a lower resolution and then gradually increasing the resolution, ingesting the solution from the previous (lower) resolution as the initial guess to the next (higher) resolution. Such an approach has been applied and proved effective in the context of grid-based atmospheric models (e.g. Heikes et al., 2013;Sandbach et al. , 2015), but a multigrid approach is difficult to implement on current spectral models because of the aforementioned non-nested nature of Gaussian quadrature. Another spectrally accurate quadrature rule exists, however, that, unlike the conventional Gaussian quadrature, has nestable nodes that would allow for straightforward implementation of multigrid approach. This quadrature rule, introduced by Clenshaw and Curtis (1960), has been well known in the field of Numerical Analysis and is shown by some authors to have several advantages over the classical Gaussian quadrature. Although Clenshaw-Curtis quadrature is not optimal in the sense of giving exact integration for polynomials of up to only J − 1 (as opposed to the Gaussian 2J −1) degrees given J nodes, it has been shown to be practically as accurate as Gaussian quadrature in many applications (e.g. Trefethen , 2008, and the references therein). Despite gaining popularity in numerical analysis, to the authors' best knowledge, Clenshaw-Curtis quadrature appears not to have been used in spectral transform models, at least in atmospheric modelling.
This short paper aims to show that the Clenshaw-Curtis quadrature, with its associated nodes (which turn out to be just equispaced latitude grids, see (13)), can be used as an alternative to the classical Gaussian quadrature with the Gaussian latitude grids in atmospheric global spectral models. Although our ultimate goal is to investigate the effectiveness of a multigrid approach in a global non-hydrostatic spectral model defined on the multigrid-friendly Clenshaw-Curtis grid, we limit the scope of the present paper to only showing numerical soundness of the quadrature applied to the associated Legendre functions and the equivalence of the Gaussian and Clenshaw-Curtis quadrature implemented on a shallow-water equations (SWE) model and a three-dimensional dry hydrostatic primitive equations (HPE) model.
The rest of the paper is structured as follows: Section 2 reminds the reader with how the quadrature is used in global spectral atmospheric models. Section 3 introduces Clenshaw-Curtis quadrature in the context of spectral transform. Section 4 examines the numerical orthonormality of the associated Legendre functions evaluated with Clenshaw-Curtis quadrature with a comment on aliasing of quadratic and cubic terms. Sections 5 and 6 document the implementation of Clenshaw-Curtis quadrature to the spherical SWE model and the HPE global dynamical core, and compare the results produced with Gaussian and Clenshaw-Curtis quadrature in the context of standardised test cases. Section 7 concludes the paper with an outlook for our future plans.

Discrete spherical harmonics transforms and Gaussian quadrature
Global spectral models represent atmospheric state variables in both gridpoint space and wavenumber (spectral) space and transform them back and forth during the course of time integration (Orszag , 1970). The discrete inverse spherical harmonics transform, or synthesis, transforms the state variable from the spectral representation X m n to its gridpointspace representation X(λ i , φ j ) by the following formula: where N is the truncation total wavenumber, n and m are, respectively, the total and zonal wavenumber, λ is the longitude, φ is the latitude, Y m n (λ, φ) is the spherical harmonic function with the total and zonal wavenumbers of n and m, respectively. Throughout this manuscript, the triangular truncation as in (1) is assumed. In the Japan Meteorological Agency (JMA)'s Global Spectral Model (JMA-GSM), as in other spectral models, the transform of (1) is implemented via two steps by first performing the associated Legendre transform in the meridional direction from X m n to and then performing inverse discrete Fourier transform (DFT) in the zonal direction from X m (φ j ) to X(λ i , φ j ) := N m=−N X m (φ j )e √ −1λ i . Here,P m n (sin φ j ) denotes the associated Legendre polynomial of degree n and order m, normalised to satisfy the following orthonormality: The longitudinal grids {λ i }, i = 0, · · · , I − 1 are chosen so that the DFT can be computed efficiently via a Fast Fourier Transform (FFT) algorithm, resulting in equispaced grids λ i = 2πi/I.
The direct spherical harmonics transform, or analysis, recovers spectral coefficients X m n from the gridpoint-space representation X(λ i , φ j ), by exploiting the orthonormality of the spherical harmonics: The integrations in (5-6) have to be evaluated numerically with some quadrature rules. For the zonal direction, direct DFT gives exact integration to an integrand with wavenumbers up to I given I zonal grids. Thus, for each meridional gridpoint φ j , (6) can be discretised as with the constraint to avoid aliasing introduced by quadrature error (note that the integrand of (6) contains zonal wavenumbers of at most 2N because X(λ, φ) is a truncated sum of Y m n with |m| being at most N). For the meridional direction, any quadrature rule can be used as long as it gives exact integration to polynomials of sin φ to a desired degree. In practice, however, Gaussian quadrature (with the ordinary Legendre polynomials as the basis set) is almost always employed due to its optimality in the sense described in the first paragraph of Section 1. With Gaussian quadrature, (5) is discretised as with the Gaussian weights given by and the latitude grids (the Gaussian latitudes) {φ j }, j = 1, · · · , J given as the J roots of P J (sin φ) = 0 where P J (sin φ) is the ordinary Legendre polynomial of degree J (NIST , 2016). In order to avoid aliasing due to quadrature errors, there should be at least latitude grids for the given truncation total wavenumber N . This is because the integrand of (5) is a polynomial of sin φ of at most 2N degree 1 . This condition, along with (8), allows exact discrete transforms with linear truncation (i.e., J ≈ N).

Clenshaw-Curtis quadrature
In numerically evaluating the integration in (5), it is possible, though apparently never have been tried with global spectral models, to use a quadrature rule other than the classical Gaussian formula (9-10). In this paper we focus on using Clenshaw-Curtis quadrature (Clenshaw and Curtis , 1960) in place of the standard Gaussian quadrature, motivated by its nestable property that we described in Introduction.

Formulation
For convenience, we denote the integrand of (5) by where we have introduced the colatitude θ = π 2 −φ. Clenshaw-Curtis quadrature first expands g(cos θ) by the Chebyshev polynomials T l (cos θ) = cos lθ of degrees l = 0, 1, · · · , J − 1 and then analytically integrates them termby-term. Expansion of g(cos θ) into Chebyshev polynomials is performed via the type-I discrete cosine transform (DCT) without any errors but rounding. From this constraint, the collocation points θ j (or φ j ) are chosen as 1 Note that, in spectral methods, X m (φ) is assumed to be a linear combination of P m n ′ (sin φ) (|m| ≤ n ′ ≤ N ) and thatP m n (sin φ) takes the form of 1 − sin 2 φ |m|/2 × (polynomial of sin φ of degree at most n − m).
Here the poles (j = 0 and j = J + 1), which are required to be in the collocation points to perform type-I DCT, are omitted in (13) because, as we eventually see in (16), the quadrature weights given to the poles are zero. Boyd (1987) showed that, instead of using DCT code, the quadrature can be cast in the following standard quadrature form: with the weights w CC j given by Boyd (1987) did not provide explicit derivation of (16); we give a concise derivation in Appendix. Equation (15) and its Gaussian counterpart (9) take the same form, which means that only minor code adaptation is necessary to implement Clenshaw-Curtis quadrature on any model based on Gaussian quadrature.
An important difference from Gaussian quadrature is that the number of nodes, given the truncation total wavenumber N, required for alias-free exact meridional integration, is which is about twice larger than in the Gaussian case (11). This is because Clenshaw-Curtis quadrature expands the integrand into Chebyshev polynomials of up to (J − 1)-th degree and thus J − 1 needs to be no less than 2N, the maximum degree of the integrand polynomial (see the footnote below (11)). Conventionally, in the context of spectral modelling, truncation rules with J ≈ 2N, J ≈ 3 2 N and J ≈ N are referred to, respectively, as "cubic," "quadratic" and "linear" truncation, which are so named because the quadratic truncation, for instance, avoids aliasing when a quadratic term computed in the grid space is transformed back to spectral space (and likewise for cubic and linear truncation). With this nomenclature, the condition (17) requires Clenshaw-Curtis quadrature to be used with "cubic" (or higherorder) truncation to assure exact spectral transforms on linear terms. Strictly speaking, it is perhaps an abuse of terminology to apply this nomenclature to Clenshaw-Curtis quadrature since, for example, "linear" truncation does not guarantee alias-free transform on linear terms. In this manuscript, we nevertheless adhere to this convention since, as we show in section 4.2, aliasing errors that arise from inexact quadrature are orders of magnitude smaller than those that arise from sub-sampling in grid space and hence pose little problem in practice.
It may seem that having to use cubic truncation is too much a restriction for Clenshaw-Curtis quadrature to be usable in an atmospheric model. Recent findings in the context of very high-resolution atmospheric simulation suggest, however, that this is perhaps not really a serious restriction: Wedi (2014) for instance showed, through experimentation using the spectral atmospheric model of the European Centre for Medium-Range Weather Forecast (ECMWF), that the aliasing errors coming from the quadratic and cubic (or even higher-order) terms on the right-hand-side of the governing equations get larger as the horizontal resolution increases, so that, at a resolution as high as ∼ 10 km grid spacing, cubic truncation (e.g., Tc1023), which automatically filters out quadratic and cubic aliasing, yields more accurate forecasts than the linear truncation with the same grid spacing (e.g., Tl2047) does. ECMWF in fact adopted cubic truncation at the total wavenumber of 1279 (∼ 8 km grid spacing) into their operational suite in 2016 (Malardel et al. , 2016). It is thus likely that future high-resolution global spectral models adopt cubic or higher-order truncation, in which case the condition (17) is not a particular disadvantage.

Lower resolution spherical harmonics transforms at a reduced cost
The advantage of using Clenshaw-Curtis quadrature is that the quadrature points nest. This not only allows for straightforward multigrid implementation in gridpoint space but also allows spherical harmonics transforms with a lower truncation wavenumber to be performed economically. Suppose, for instance, the model has Tc479 (triangular cubic truncation at total wavenum-ber 479) horizontal resolution with Clenshaw-Curtis grid, and one wishes to compute spectral coefficients of total wavenumber up to 239 (the lower half of the spectral space) from the physical data on the Tc479 grid. For convenience let N full + 1 = 480, N half + 1 = 240 and assume that, for a truncation wavenumber N, there are J(N) = 2N + 1 latitudinal gridpoints. Then, the colatitudes θ j of Tc479 grid are with j = 1, 2, · · · , 2N full + 1. If we take every other latitudinal gridpoints starting from j = 2, the colatitudes of this subset are with j ′ = 1, 2, · · · , 2N half + 1, which are identical to the colatitudes of Tc239 Clenshaw-Curtis grids. The same nested property also applies to the equispaced longitudinal grid as long as the number of longitudinal gridpoints is chosen to be multiple of two. Thus, with Clenshaw-Curtis quadrature, Tc479 grid contains Tc239 grid as its complete subset, enabling us to compute spectral coefficients of the total wavenumber up to 239 using direct (grid-to-wave) transform code of Tc239 model.
4 Numerical accuracy of spectral transforms based on Clenshaw-Curtis quadrature 4.1 Numerical orthonormality of the associated Legendre functions Spectral transform methods (Orszag , 1970) transform spectral and gridpoint space back and forth on each time step by using the direct and inverse spherical harmonics transforms. Because the transforms are repeated many times during the course of time integration, discrete transforms need to be nearly exact (i.e., the errors must be within the range of rounding errors). In particular, the orthonormality of the associated Legendre functions, (3), with the integral evaluated by a quadrature rule, needs to strictly hold for all m and n within the truncation limit. In the previous section we have seen that, in theory, the strict quadrature should be guaranteed by imposing the condition (17). To verify the numerical exactness on an actual implementation, this subsection examines the numerical orthonormality of the associated Legendre functions (3) evaluated with Clenshaw-Curtis quadrature implemented on JMA-GSM. Here, we define normality error and orthogonality error, respectively, as: where φ j and w j are the nodes and weights defined for each quadrature rule. We examine ε m n,N and ε m n,O for different combinations of quadrature and truncation rules. Clenshaw-Curtis quadrature requires cubic or higherorder truncation (17). We thus focus mainly on cubic truncation with J = 2N + 1 meridional gridpoints for a given truncation wavenumber N. Given the observations from the literature (e.g., Trefethen , 2008), that Clenshaw-Curtis quadrature is nearly as accurate as the Gaussian in many practical applications, it is still interesting to look at linear (J ≈ N) and quadratic (J ≈ 3N/2) truncation.
Throughout this section, the discrete Legendre transforms are computed using the code of JMA-GSM (Miyamoto , 2006). The normalised associated Legendre functionsP m n (sin φ j ) are computed by the three-point recurrence (Eqs. (5-9) of Enomoto , 2015) with quadruple precision arithmetic and then cast to double precision. Computation of (20) and (21) is performed using double precision. JMA-GSM employs the reduced spectral transform (Miyamoto , 2006;Juang , 2004) in which the terms involvingP m n (sin φ j ) whose absolute value is below 10 −16 are neglected from the summations in (2) and (9). We have repeated all computations with this option on or off and found that this option does not affect the results shown in this section. Here we only show the results obtained with this option switched on. Figure 1 shows the normality error and orthogonality error for Clenshaw-Curtis quadrature with the resolutions of Tc479 (triangular cubic truncation at N = 479 total wavenumber), Tq639 (triangular quadratic truncation at N = 639 total wavenumber) and Tl959 (triangular linear truncation at N = 959 total wavenumber), all with J = 959 meridional gridpoints. As the theory predicts, with cubic truncation, both normality error (panel a) and orthogonality error (panel d) are within the range of double-precision rounding (≤ 10 −16 ) for all combinations of m and n. In contrast, with lower-order truncation (panels b and c), normality error exceeds the level of rounding error for total wavenumbers n > (J −1)/2 = 479, with a tendency to become larger for larger n and smaller m. This is consistent with the condition (17) because the integrandP m n (sin φ) 2 is a degree-2n polynomial of cos θ since it takes the form sin 2m θ×[polynomial of cos θ of degree 2(n − m)] = (1 − cos 2 θ) m × [polynomial of cos θ of degree 2(n − m)] = [polynomial of cos θ of degree 2n]. Interestingly, however, the normality error remains < 10 −3 even when the quadrature limit (17) is violated.
Clenshaw-Curtis quadrature guarantees orthogonality only to total wavenumbers n ≤ J − N − 1 because the degree of the integrandP m n (sin φ)P m n ′ (sin φ) in (21), which is n + n ′ ≤ n + N, need to be less than J. Consistently, orthogonality is exact with cubic truncation (panel d) but not with quadratic truncation (panel e) for the total wavenumbers n ≥ J − N + 1 = 320. With linear truncation (panel f) it is not exact even for n = 0. Interestingly again, the orthogonality error nevertheless remains < 10 −3 for any pairs of (n, m) that violate the quadrature limit.
From these results we can conclude that Clenshaw-Curtis quadrature assures exact spectral transforms if it is used with cubic (or higher-order) truncation. This conclusion holds to all other resolutions that we tested.

Aliasing of quadratic and cubic terms evaluated on cubic Clenshaw-Curtis grids
Clenshaw-Curtis quadrature with cubic truncation does not guarantee aliasfree spectral transforms for quadratic and cubic terms. As such, before applying Clenshaw-Curtis quadrature with cubic truncation to nonlinear models, we need to quantitatively assess the degree of aliasing for nonlinear terms.
In this subsection we quantify aliasing errors for quadratic and cubic terms associated with spherical harmonics transforms and compare the Clenshaw-Curtis quadrature with cubic truncation and the Gaussian quadrature with linear truncation. We quantify the aliasing errors by the following procedure. Consider the following series of operations: 1. Produce three scalar fields U, V , and W on the sphere by randomly choosing their spectral representations U m n , V m n and W m n . The spectral coefficients are chosen such that their (m, n)-component, 0 ≤ n ≤ N trunc = 479, |m| ≤ n, has both its real and imaginary part taken as an independent pseudo-random draw from a uniform distribution over [1 × (n + 1) −1/3 , 2 × (n + 1) −1/3 ].
2. Transform the scalar fields U m n , V m n , W m n in spectral space into grid space by (1) , and then compute, in grid space, quadratic and cubic terms 3. Finally transform the quadratic and cubic terms in grid space X q and X c back to spectral space by (5-6) to yield X m q,n and X m c,n .
The scaling of the pseudo-random numbers by (n + 1) −1/3 in the step 1 is intended to mimic the power spectra of physical quantities typically encountered in geophysical applications (for instance, the power spectra of kinetic energy in three-dimensional turbulence in the inertial range follow n −5/3 -law, which corresponds to the spectral coefficients of divergence or vorticity being proportional to ∼ n −1/3 ). This choice is arbitrary and subjective, but we confirmed that it does not sensitively affect the conclusion by repeating the calculation changing the power from −1/3 to −3, −2, −1, −1/2 and 0. We first perform the steps 1-3 on Tl1439 Gaussian grid (N = 1439, J = 1440) using Gaussian quadrature and a linear truncation at the total wavenumber of N ref = 1439(≈ 3N trunc ). This resolution is high enough to produce aliasing-free reference solutions of both quadratic and cubic terms, which we denote respectively by X m,ref q,n and X m,ref c,n . We then repeat the steps 1-3 using the test resolutions and quadrature rules (Tl479 Gaussian or Tc479 Clenshaw-Curtis), using the same values of U m n , V m n and W m n as used to produce the reference solutions. The spectral coefficients of the quadratic and cubic terms computed with the test resolution and quadrature rule, denoted X m,test The results for Gaussian linear grid (Tl479; N = 479, J = 480) and Clenshaw-Curtis cubic grid (Tc479; N = 479, J = 959) are plotted in Figure  2. Clenshaw-Curtis quadrature performed on cubic grid (Figure 2, right panels) produces normalised aliasing errors below 10 −5 (in the quadratic case, Figure 2b) or 10 −8 (in the cubic case, Figure 2d) for most pairs of (n, m) except in a small region with large n/m ratio ( 10). This is in contrast to the case of Gaussian quadrature performed on linear grid ( Figure  2, left panel), where the normalised errors are almost anywhere greater than 10 −4 (in the quadratic case, Figure 2a) or 10 −3 (in the cubic case, Figure 2c), for some pairs of (n, m) even exceeding 1 in the quadratic case.
These results indicate that the aliasing errors that arise from the quadrature rule not being exact (as seen in Figure 2b,d) are orders of magnitude smaller than the aliasing errors that result from sub-sampling of waves in gridpoint space (as seen in Figure 2a,c). From this we may deduce that the aliasing from nonlinear terms in a model with Clenshaw-Curtis quadrature should remain under a controllable level as long as the quadrature is used with cubic truncation. This conclusion is further corroborated by the results shown in the next two sections.

Implementation to a shallow-water equation model
In implementing a new scheme, it is convenient to test it with a simpler model before introducing it to the full three-dimensional model. As a testbed, we first implemented Clenshaw-Curtis quadrature to a semi-Lagrangian shallowwater equations (SWE) model. The main focus of this paper is to examine the consistency between the models with Clenshaw-Curtis and Gaussian quadrature. Accordingly, the validations that follow are not designed to test the performance of the model itself but rather to detect any discrepancies arising from the use of different quadrature rules.

The model
The model used in this section is built by adapting the code of JMA-GSM. The governing equations are the advective form of the SWE on a rotating sphere described in section 2.2 of Williamson et al. (1992) but with the Coriolis terms incorporated in the left-hand side of the momentum equations as in Temperton (1997). The governing equations in advective form where r is a positive integer. This is numerically solved implicitly in spectral space by applying the following operation at the end of each time step: (X m n ) after = 1 1 + 2K∆t n(n+1) a 2 r (X m n ) before .
where n, m, a and ∆t are, respectively, the total and zonal wavenumbers, the Earth's radius and the time step. We modify the model to allow an option to use Clenshaw-Curtis quadrature and compare the model runs produced with this option and the runs produced with the default Gaussian quadrature, under the framework of two standard test cases. The models with the two different quadrature rules run on different grids. To allow for accurate comparison, the Clenshaw-Curtis version of the model is adapted to output model states also on the Gaussian grids by converting the model state variables in spectral space to grid space using the inverse spherical harmonics transform (1) defined on the Gaussian grids. To ensure that the two versions start from exactly identical initial conditions, the initial conditions for the Clenshaw-Curtis version are produced by first computing them in grid space on the Gaussian grids, converting them to spectral space by the direct spherical harmonics transform (4) defined on the Gaussian grids, and then converting them back to grid space (but this time on the Clenshaw-Curtis grids) by the inverse spherical harmonics transform (1) defined on the Clenshaw-Curtis grids.

Test cases
The first test case examined is the evolution of the initially zonal flow over an isolated mountain proposed by Williamson et al. (1992) as Case #5. This test case is repeated for the resolutions of Tc31, Tc63, Tc127 and Tc255, both for Clenshaw-Curtis and Gaussian version of the model. No explicit numerical diffusion of the form (24) is used in any of the integrations for this case.
The solution of Williamson et al. (1992) test case #5 appears to be dominated by relatively low wavenumber components, so it may not be suitable for assessing model behaviour at higher wavenumbers. The second test case we examine, the barotropic jet instability test case proposed by Galewsky et al. (2004), allows us to assess the models in a more realistic, multi-scale flow situations. Galewsky et al. (2004) reported that weak explicit dissipation in the governing equations is necessary to obtain converged solutions, and suggested to apply harmonic diffusion (r = 1 and K = 1.0 × 10 5 m 2 in (24)) to the governing equations to facilitate model comparison. We followed this suggestion and performed the test with a relatively high resolution of Tc479 (∆x ∼ 20 km at the Equator).

Results: flow over an isolated mountain
In interpreting the results, it should be noted that, unlike models with Eulerian advection scheme, exact agreement between Clenshaw-Curtis and Gaussian versions of the model is not expected in our semi-Lagrangian model since the two versions use different grids that result in differences in departure point interpolation. Nevertheless, we found, by plotting overlaid contour maps of relative vorticity (or any other) fields from the two versions of the model (not shown), that the two quadrature rules result in visually indistinguishable integrations up to at least 288 hours, for any of the tested resolutions. To quantitatively assess the difference, we computed the normalised L 2 -differences defined as: for zonal wind (ξ = u), meridional wind (ξ = v) and height field (ξ = h), for integration times t from 0 to 288 hours at 6-hourly intervals. Here the superscripts CC and G represent the integration, respectively, by Clenshaw-Curtis and Gaussian versions of the model, and the the spherical integral I(ξ) of the square of a variable ξ is numerically evaluated with Gaussian quadrature using ξ CC and ξ G both defined on the Gaussian grid. The normalised L 2 -differences of the zonal wind fields (Figure 3) are generally very small. They grow as the integration length gets longer but seem to saturate by ∼ 10-day integration. Higher resolutions result in smaller differences, presumably due to reduced interpolation errors in semi-Lagrangian advection. At saturation, the normalised differences are ∼ 10 −3 with Tc31 resolution and ∼ 10 −5 with Tc255 resolution, which are in practice negligibly small. Similar results were also observed for v and h fields (not shown).

Results: Barotropic jet instability
With sharp frontal structures present in the solution, in the barotropic jet instability test, one may expect to see a larger discrepancy between the two versions of the model with different quadrature rules and the associated grids. We observed, however, that the two versions yielded results that are visually indistinguishable (Figure 4a), even at integration length as long as 144 hours (6 days) by which time multiple rolled-up vortical structures develop. The difference between the two can be found mainly in the frontal regions with strong vorticity gradients (Figure 4b), but the differences are smaller than the raw values by more than five orders of magnitude (note the different contour intervals in panels (a) and (b) of Figure 4). The normalised L 2 -difference for any of u, v and h field (not shown) were below 10 −7 , again indicating that the discrepancies arising from use of different quadrature rules and the associated grids are negligible in practice.  6 Implementation to a hydrostatic primitive equations model

The model
Clenshaw-Curtis quadrature is implemented to the dry dynamical core of JMA-GSM. It is an HPE model with Semi-Implicit Semi-Lagrangian (SISL) time discretisation. The horizontal discretisation is spectral with spherical harmonics as the basis functions, and the vertical discretisation employs finite differencing on σ − p hybrid coordinate (Simmons and Burridge , 1981) with 100 layers extending from the surface up to 0.01 hPa. Detailed description of the model's dynamical core can be found in Section 3.1 of Yukimoto et al. (2011) and Section 3.2.2 of JMA (2013).
As with most other SISL spectral models, by default JMA-GSM uses linear truncation. Since Clenshaw-Curtis quadrature requires cubic (or higher) grid truncation, in this study JMA-GSM is modified to allow for cubic truncation. All the results shown in this section are produced with Tc159 horizontal resolution (∆x ∼ 60 km at the Equator) and a time step ∆t = 1200 s.

Test case specification
To assess any discrepancies in the forecasts that result from the use of Clenshaw-Curtis quadrature instead of the standard Gaussian quadrature, we perform the standardised test case proposed by Jablonowski and Williamson (2006) and adopted by numerous studies that implement new schemes to global atmospheric dynamical cores. The Jablonowski-Williamson test case comprises two parts. The first part, the steady-state test, initialises the model with an analytic, baroclinically-unstable steady-state solution to the HPE, and tests the model's ability to maintain this steady-state initial condition. The second part, the baroclinic wave instability test, initialises the model with the same steady-sate solution but with a localised small-amplitude perturbation in the zonal wind field superimposed on the westerly jet axis located in the Northern Hemisphere mid-latitude. The small perturbation eventually develops into a train of unstable baroclinic waves that experiences nonlinear break-down by ∼ day 9. This second test allows us to examine the model's behaviour in a situation akin to a typical mid-latitude synoptic condition with sharp fronts and fine-scale eddies.

Steady-state test
We first examine the impact of quadrature choice on the model's ability to maintain the steady-state solution by comparing, between the Clenshaw-Curtis and Gaussian versions of the model, the two measures of steady-state maintenance proposed by Jablonowski and Williamson (2006). The first measure, denoted as l 2 (u(t) −ū(t)) and defined as the square root of the three-dimensional integral of zonally asymmetric part of the zonal wind u field (see Eq. (14) of Jablonowski and Williamson (2006)), quantifies to what extent the model is able to keep the zonally symmetric structure of the steady initial state. The second measure, denoted as l 2 (ū(t) −ū(t = 0)) and defined as the square root of the three-dimensional integral of the departure of the zonal componentū(t) of the zonal wind at time t from that at the initial time (see Eq. (15) of Jablonowski and Williamson (2006)), complements the first measure and quantifies how much the zonal mean of u field deviates from the initial zonally symmetric steady state over the course of integration.
As shown in Figure 5, the Clenshaw-Curtis and Gaussian versions of the model result in virtually identical l 2 error measures. The relative difference between the two versions of the model (defined as the difference normalised by the error of the Gaussian version; drawn with dotted lines and to be read on the right axis) is at most 2% for the first measure ( Figure 5a) and is well below ∼1 % for the second measure (Figure 5b), indicating that the choice of the quadrature rule and the associated grid does not affect the model's ability to maintain the steady-state initial condition.

Baroclinic instability wave test
As in the barotropic jet instability test for SWE model (section 5.4, Figure 4), with sharp fronts present in the solution, we can expect the baroclinic instability wave test to better reveal discrepancies (if any) in the numerical results arising from the use of different quadrature rules and the associated grids. Just like in section 5.4, however, again we confirmed a high degree of agreement between the solutions from Clenshaw-Curtis and Gaussian versions of the model. As an example, we show the snapshots of 850 hPa temperature field at day-9 integration in Figure 6a. The results from Clenshaw-Curtis version (thick gray contours) and Gaussian version (thin black contours) are indistinguishable on this plot, and the difference between the two (Figure 6a exhibit noisy patterns localised along the fronts and with very small amplitude (at most 0.13 K), indicating that no systematic difference results from the use of different quadrature rules.

Conclusion and future directions
Traditional spectral global atmospheric models adopt Gaussian quadrature, and the non-nested property of its associated irregular latitudinal quadrature nodes (the Gaussian grid) makes it difficult to apply multigrid approaches to global spectral models. To facilitate straightforward use of a multigrid method in a global spectral model in the future, in this study we proposed to adopt Clenshaw-Curtis quadrature in evaluating direct (grid-to-wave) Legendre transform. With Clenshaw-Curtis quadrature, the latitudinal grids (or quadrature nodes) are equispaced, so that the model grid for an arbitrary horizontal resolution includes the model grid for half its horizontal resolution as its complete subset, allowing for implementation of a multigrid method without a need for any off-grid interpolation.
Theoretical consideration shows, and numerical computation demonstrated, that Clenshaw-Curtis quadrature guarantees exact (up to machine precision) spherical harmonics transforms if it is used with cubic (or higher-order) grid truncation. One may argue that having to move from linear or quadratic truncation, which are the norm for the current semi-Lagrangian or Eulerian spectral models, respectively, is a serious limitation, but that is not the case for high-resolution modelling since, as shown by recent studies at ECMWF (Wedi, 2014), higher-order truncation is anyhow necessary to suppress aliasing errors and the resultant spectral blocking at a resolution as high as ∆x < 10 km.
Comparison of the Clenshaw-Curtis and Gaussian versions of the shallowwater model and dry hydrostatic primitive equations model, both with cubic grid and with semi-Lagrangian advection scheme, performed under the framework of idealised standard test cases, demonstrated that Clenshaw-Curtis quadrature gives model predictions that are almost identical to those obtained with the standard Gaussian quadrature.
Our ultimate goal is to ensure computational efficiency of a global spectral non-hydrostatic atmospheric model by enabling a multigrid approach, and this study serves as a first step in this direction. The success (or not) of such an approach is not clear, particularly because the future evolution Figure 6: (a) Snapshot of 850-hPa temperature at day 9 from the baroclinic instability wave test case produced with Clenshaw-Curtis (thick gray lines) and Gaussian (thin black lines) versions of JMA-GSM dynamical core run at Tc159 resolution, and (b) the difference between the two. Contour intervals are 4 K in (a) and 0.02 K in (b). In (b), positive and negative values are drawn, respectively, with solid and dashed contours, and the zero contours are omitted. Note that in (a) the two sets of contour lines are completely superimposed.
of the high-performance computing (HPC) architecture is still in a state of flux. It has long been predicted that spectral models will face a serious challenge on a massively parallel architecture because their current algorithms require global inter-node communications on every time step. To maintain admissible scalability on future architectures, therefore, global models would have to abandon, or at least reduce reliance on, spectral transforms. In this line, ECMWF has recently introduced a structured grid system, called 'cubic octahedral reduced Gaussian grid' that assures exact spectral transforms via Gaussian quadrature (Malardel et al. , 2016) while at the same time allowing for straightforward grid-based discretisation of horizontal derivatives that only requires neighbouring-node communications , whereby paving a path for a smooth transition from spectral to grid-based modelling (or hybridisation of the two methods). Their approach of allowing both spectral and grid-based discretisations on a single grid system is ingenious and appealing, and we assert that Clenshaw-Curtis quadrature, with its equispaced latitude grids, can be combined to yield cubic octahedral Clenshaw-Curtis grid (Figure 7a).
Compared to the Gaussian octahedral grid, the Clenshaw-Curtis octahedral grid has the advantage of allowing for a straightforward use of multigrid in gridpoint-space solution of elliptic equations that result from a (semi-)implicit time stepping. Moreover, it covers the globe entirely with triangular meshes (Figure 7a), unlike the Gaussian counterpart which has rectangular cells along the Equator. Its latitudinally-uniform grid spacing may also help to make the derivative stencils somewhat simpler.
Alternative grid alignments based on reduced Clenshaw-Curtis grid are also possible that cover the globe with triangular meshes. Figure 7b shows one such possibility, which may be called cubic icositetrahedral (24-face polyhedral) Clenshaw-Curtis grid, where the globe is divided into 24 triangular panels of equal area, 6 covering the Northern Hemisphere extratropics (from 90°N down to 30°N), 12 covering the tropical 30°N-30°S belt, and the remaining 6 covering the rest (from 90°S up to 30°S). Compared to the octahedral grid, having to shift the phase in zonal Fourier transforms in the tropical belt introduces additional complexity in the code, but the icositetrahedral grid gives more uniform grid spacing across latitudes and is more compatible (Figure 8) with the grid reduction rule of Miyamoto (2006) which guarantees numerically exact spherical harmonics transform by requiring, for each latitude φ j , that there be at least 4M j + 1 longitudinal points (with cubic truncation) where M j is defined as the smallest zonal wavenumber for which P m n (sin φ j ) < 10 −16 (∼ double precision machine epsilon) holds for any 0 ≤ n ≤ N and |m| ≤ M j . The effectiveness in terms of accuracy, computational efficiency on a massively parallel machine architecture, and ease of transition from the current spectral model, of using octahedral and/or icositetrahedral Clenshaw-Curtis grids, along with gridpoint-space-based multigrid methods, will be explored in our future project. Appendix: Derivation of the Clenshaw-Curtis quadrature weights w CC j Let g(cos θ) be the integrand as defined in (12). In Clenshaw-Curtis quadrature, g(cos θ) is expanded by the Chebyshev polynomials of degrees up to J − 1, giving g(cos θ) sin θ = J−1 l=0 a l T l (cos θ) sin θ = J−1 l=0 a l cos lθ sin θ. (A1) The original formulation by Clenshaw and Curtis (1960) obtains the expansion coefficients a l by performing DCT on g(cos θ) and then analytically integrates each summand in the rightmost-hand side of (A1). The alternative derivation by Boyd (1987) notes that {cos lθ sin θ} J−1 l=0 spans the same space as {sin pθ} J p=1 since cos lθ sin θ = 1 2 {sin (l + 1)θ − sin (l − 1)θ} (note that p = 0 is dropped because the function sin 0θ is identically zero), and expands g(cos θ) sin θ rather than g(cos θ) itself into a linear combination of "trigonometric cardinal functions" {C j (θ)} J j=1 : g(cos θ) sin θ = J j=1 {g(cos θ j ) sin θ j } C j (θ) where {C j (θ)} J j=1 are the transformed basis of the space spanned by {sin pθ} J p=1 that satisfies for any i, j ∈ {1, 2, · · · , J} with the nodes θ j chosen as in (13). As shown in the Appendix of Boyd (1987), by exploiting the discrete orthonormality of {sin pθ} J p=1 on the points in (13), the cardinal functions C j (θ) can be expressed as C j (θ) = 2 J + 1 J p=1 sin pθ j sin pθ.
The explicit expression of the quadrature weights w CC j as shown in (16) follows from plugging (A4) into (A2) and then carrying out the integration from θ = 0 to θ = π.