## Abstract

A dynamical core of a general circulation model with the spectral method using double Fourier series (DFS) as basis functions is presented. The model uses the hydrostatic balance approximation and sigma coordinate system in the vertical direction and includes no topography. The model atmosphere is divided into 25 layers with equal sigma depths. Prognostic equations for the vorticity, divergence, temperature, and logarithmic surface pressure are solved by the DFS spectral-transform method with the Fourier filtering at middle and high latitudes. A semi-implicit time-stepping procedure, which deals with the eigendecomposition and inversion of the 3D Helmholtz equation associated with the gravity wave terms, is incorporated for the gravity wave–related terms. The DFS model is tested in terms of the solution of the 3D Helmholtz equation, balanced initial state, developing baroclinic waves, and short- and long-term Held–Suarez–Williamson simulations for T42, T62, T84, and T106 resolutions. It is found that the DFS model is stable and accurate and produces almost the same results as the spherical harmonics method (SHM). The normalized difference (i.e., *L*_{2} norm error) measured from the results of highest-resolution SHM-T106 showed a desirable convergence of the DFS solution with the resolution. The convergence property, however, varies with the test case and prognostic variables. The total mass (or global integrated surface pressure) is conserved to a good approximation in the long-term simulations. Computing on the high-performance computer NEC SX-5 (parallel-vector architecture) indicated that DFS is more efficient than the SHM and the efficiency increases with the resolution, for example, by factors of 2.09 and 7.68 for T212 and T1022, respectively.

## 1. Introduction

The spectral method that incorporates the associated Legendre functions as an orthogonal basis is widely used in the global weather forecast models (e.g., Gates 1992; LLNL 2001) because of its many advantages such as accuracy, stability, uniform resolvability over the sphere, and easiness of semi-implicit time stepping (Hoskins and Simmons 1975; Haltiner and Williams 1980; Kanamitsu 1989; Spotz et al. 1998). As the resolution of the model increases, however, the computational cost for the spectral transform becomes dominant over other parts of the spectral model because it requires *O*(*N* ^{3}) operations, with *N* being the truncation of the model. Including this point, some detailed description of the minor points inherent in the traditional spectral method can be found in Williamson et al. (1992) and Spotz et al. (1998). In spite of the computational cost, the spectral method is applied to the operational forecast models with a very high resolution, for example, European Centre for Medium-Range Weather Forecasts medium-range forecast models [refer to White (2002) for T230 as of April 2002 and the Web site http://www.ecmwf.int/products/forecasts/guide for the current model with T512]. Although an excessive longitudinal resolution at high latitudes is inevitable for the spectral method, the use of reduced grids at high latitudes contributes to a reduction of the computational time (Hortal and Simmons 1991) without deteriorating the accuracy. In fact, White (2002) found that even for a high-resolution model such as T230 the computational efficiency of the spectral method is comparable to the semi-Lagrangian method (Côté and Staniforth 1990; Layton and Spotz 2003), which allows a large time step.

Recently, a new class of dynamical cores has been developed (Ringler et al. 2000; Majewski et al. 2002; Fournier et al. 2004; Tomita and Satoh 2004) that perform with accuracy and reduced computational cost at high resolution. The most notable feature of these models is to adopt a grid system with almost uniformly distributed grids over the sphere without any pole problem arising from the spherical coordinate system. Ringler et al. (2000) and Majewski et al. (2002) introduced the 3D atmosphere model in which the icosahedral–hexagonal grid system is used for discretizing the differential equations on a spherical surface with the finite-difference method. Fournier et al. (2004) applied the spectral element method for a cubed sphere (Taylor et al. 1997; Thomas and Loft 2002) to the 3D primitive equations with hydrostatic approximation. The spectral element method can provide a local mesh refinement on the spherical surface, which is quite useful for forecasting multiscale phenomenon under a single dynamics with reduced computational cost. In their model, the scale-dependent viscosity in the form of the Laplacian operator is applied to all prognostic variables including surface pressure, and the vertical smoothing is done by fourth-order finite difference to stabilize the model. Tomita and Satoh (2004) developed a nonhydrostatic three-dimensional atmosphere model that incorporates the finite-element method with an icosahedral grid system. They showed, through a variety of test methods, that the model can be used for the prediction of 3D atmospheric flows. The model incorporates a biharmonic diffusion in the horizontal and a sixth-order finite-difference smoothing in the vertical as a numerical stabilizer.

The biharmonic horizontal diffusion is commonly used in numerical models either to control the small-scale noises resulting from nonlinear terms or to represent horizontal mixing (Kanamitsu 1989; Ringler et al. 2000; Majewski et al. 2002; White 2002; Tomita and Satoh 2004). In some cases even higher-order viscosity (hyperviscosity) is incorporated for this purpose (Jukes and McIntyre 1987; Tung and Orlando 2003). A spectral filter that incorporates a function of total wavenumber other than the integer power of the Laplacian operator is also used in numerical models (Boer and Denis 1997; Gelb and Gleeson 2001; Fournier et al. 2004). When a direct implementation of the biharmonic or higher-order Laplacian operator is not feasible, another type of filtering in terms of gridpoint values may be applied (Shapiro 1970; Nihei and Ishii 2003). It is not likely that a too-strong viscosity in the numerical model reproduces a realistic kinetic energy spectrum in the wavenumber domain.

In this study, the double Fourier series (DFS) spectral method (Cheong 2000a; Cheong et al. 2004) is extended to the three-dimensional atmosphere model with the hydrostatic approximation. The DFS retains most of the advantages available for the traditional spectral method and it performs with reduced computational cost characterized by *O*(*N* ^{2}log_{2}*N*). The DFS provides an efficient and simple high-order filtering method, which is essential for the stable time integration of nonlinear partial-differential equations for the atmospheric flow without affecting the physics of the phenomenon (Cheong et al. 2004). The DFS model is tested in terms of various dynamical phenomena, which include 3D elliptic-equation solver, balanced initial state, developing baroclinic waves, and short- and long-term simulations with simplified climatological forcing. Also addressed in this paper is the computational efficiency of the DFS dynamical core over the spherical harmonics method (SHM), which will be tested with the high-performance computer NEC SX-5, a parallel-vector architecture.

This paper is organized as follows. The next section presents the governing equations along with a brief description of numerical methods used in this study. Section 3 deals with six test cases to judge the accuracy, stability, and the computational efficiency of DFS dynamical core with the resolutions of T42, T62, T84, and T106. Discussion and conclusions are given in the final section.

## 2. Governing equations, numerical methods, and experiment set

### a. Governing equations

The primitive equations with vertical sigma coordinate on the spherical domain are used in this study (Hoskins and Simmons 1975; Haltiner and Williams 1980; Kanamitsu 1989), where the time and length are scaled by the rotation rate *Ω* and radius *a* of the earth, and temperature *T* and pressure *p* are scaled by *a*^{2}*Ω*^{2}/*R* (where *R* is the gas constant) and 1000 hPa, respectively,

and

where

The meaning of symbols is as usual: *λ* and *θ* are longitude and latitude, respectively; *u* and *υ* are the velocity components in longitude and latitude, respectively; (*U*, *V*) = (*u*, *υ*) cos*θ*, *σ* = *p*/*p _{s}*, and

*ω*,

*σ̇*, and

*f*(≡2 sin

*θ*) are the vertical pressure velocity, sigma velocity, and Coriolis parameter, respectively;

*ζ*,

*D*, and Φ represent the vorticity, divergence, and geopotential, respectively;

*x*=

*R*/

*C*, with

_{p}*C*being the specific heat capacity at constant pressure;

_{p}*T*(

*T*′) is the time-invariant layer-mean temperature (deviation from the layer mean), and

*T*is the prescribed reference temperature. Here

^{r}*ν*and

*a*are the Rayleigh damping and Newtonian cooling coefficients, respectively. The terms

*S*in the vorticity, divergence, and temperature equations denote the horizontal scale-dependent viscosity and diffusion, which are given in the form of fourth-order (∇

^{8}) spherical Laplacian. The choice of this type is to efficiently filter out the high-wavenumber components without affecting the dynamics of the flow (e.g., Cheong et al. 2004). Throughout the paper, the coefficient of the fourth-order (∇

^{8}) viscosity is given as 1.4961 × 10

^{−13}in nondimensional units for all models in order to compare the model performances in the same condition (e.g., Giraldo and Rosmond 2004; Polvani et al. 2004).

The model atmosphere extends from *p* = 0 to *p* = *p _{s}* and the vertical domain is divided in

*K*(≡25) layers with equal sigma depth. All variables are placed at the center of each layers, that is,

*σ*= Δ

*σ*(

*k*− 0.5) with Δ

*σ*= 0.04 and

*k*= 1, 2, . . . ,

*K*, except the vertical velocity

*σ̇*, which is placed at the layer boundary. The vertical boundary condition is given as

*σ̇*= 0 at

*p*= 0 and

*p*=

*p*.

_{s}The spectral method using the spherical harmonics functions to solve the Eqs. (1a)–(1e) with a semi-implicit time stepping is well illustrated in Hoskins and Simmons (1975) and Haltiner and Williams (1980), so we do not repeat the procedure in this study. The vertical differentiation is replaced by the finite-difference form of Arakawa (1972). The equations can also be solved with the DFS spectral method used for the solution of the shallow-water equations (Cheong 2000b). Unlike the spherical harmonics method, the semi-implicit time-stepping procedure for the DFS spectral shallow water model requires the inversion of tridiagonal matrices because the basis functions are not the eigensolution of the spherical Laplacian operator. The difference in basis functions also brings about the difference in the detailed procedure for the three-dimensional semi-implicit time stepping, as will be given in the next section.

### b. Spectral methods with spherical harmonics and DFS

The SHM and DFS dynamical cores described below are the research models of Pukyong National University used to study large-scale atmospheric dynamics and climate sensitivity. Common to both the SHM and DFS spectral models, a horizontal function is represented as a Fourier expansion in zonal direction. Let *X _{m}*(

*θ*) be the coefficient of zonal Fourier transform of a horizontal function

*X*(

*λ*,

*θ*) with the subscript

*m*being the zonal wavenumber,

where *M* represents the maximum zonal wavenumber. Then, for the spherical harmonics method the zonal Fourier transformed variable is expanded with the Legendre functions *P*^{|m|}_{l}(*μ*) with *μ* = sin*θ*:

For the double Fourier method, the half-ranged sine or cosine functions are used for spectral expansion in meridional direction:

where *ϕ* = *θ* + *π*/2. The spectral coefficients are obtained by the Gauss–Legendre quadrature formula for the SHM (e.g., Haltiner and Williams 1980) and the double Fourier transform for the DFS (Cheong 2000a), respectively. In Cheong (2000b), Legendre polynomials were used to solve the elliptic equation for *n* equal to even of zonal wavenumber zero. In this study, however, all differential equations associated with the dynamical core are solved with DFS as in Cheong et al. (2004).

Detailed procedures associated with the spectral representation of the governing equations can be found in Hoskins and Simmons (1975) and Haltiner and Williams (1980) for the SHM, and Cheong (2000a, b) and Cheong et al. (2004) for the DFS. For both models the transform method is used for evaluating the nonlinear terms (Orszag 1970; Cheong 2000a), but the detailed procedure for the DFS model is slightly different from the SHM model because the series expansion in (5) includes a parity function. To enhance the computational efficiency and conservation property of first-order variable, a spectral transform method using the recurrence of associated Legendre functions is introduced for SHM, as shown in Rivier et al. (2002). A Fourier filtering is applied to the coefficients of zonal Fourier transform, *X _{m}* in (4) and (5) (Hortal and Simmons 1991; Spotz et al. 1998; Cheong 2000a). The Fourier filtering may not be necessary for the SHM because of the uniform resolvability over the sphere. Nevertheless, this is done for both models to compare the results in the same condition. For the fast Fourier transform (FFT) between the grid space and wave space in one or two dimensions, the prime-factor algorithm FFT is used (Temperton 1985).

The simulations are carried out for four horizontal resolutions, which are T42 (*M* = 42, 128 × 64 transform grids), T62 (*M* = 62, 192 × 96), T84 (*M* = 84, 256 × 128), and T106 (*M* = 106, 320 × 160). The prognostic equations are time integrated with the leapfrog method, and the Asselin time filter with a coefficient of 0.064 is applied at every time step (Asselin 1972). The relaxation forcing *T ^{r}* in (1c), which is incorporated in Held–Suarez–Williamson simulations, is a function of the pressure, and since the sigma coordinate is used in the vertical direction the relaxation forcing should be calculated as a function of sigma at every time step.

As a measure of difference between two gridpoint variables, either root-mean-square difference (rmsd) or *L*_{2} error norm (*E*) is used:

where *q* is a variable to be compared, *q _{r}* is the reference solution, and

*N*is the total number of grids used for the error evaluation. The summation is done for one-, two-, or three-dimensional grid space.

_{t}## 3. Tests of DFS spectral dynamical core and comparison with SHM

### a. Case I: Inversion of 3D Helmholtz equation for the semi-implicit time-stepping method

The gravity wave–related terms are treated with implicit time stepping to alleviate the restriction on time step size, which yields the three-dimensional elliptic equation for the divergence:

where 𝗕 is a *K* × *K* real-valued matrix, *γ* = Δ*t*^{2} with Δ*t* the time step size, and **D** and **G** are column vectors with *K* elements. As was shown in Hoskins and Simmons (1975), Eq. (6) can be written as matrix equations for each spherical harmonics component in the SHM because the spherical harmonics functions are the eigenfunctions of the Laplacian operator: (𝗜 − *γc _{l}*𝗕)

**D**

*,*

_{l}*=*

_{m}**G**

*,*

_{l}*, where 𝗜 is unit matrix,*

_{m}*c*= −

_{l}*l*(

*l*+ 1), and

**D**

*,*

_{l}*and*

_{m}**G**

*,*

_{l}*are column vectors consisting of spectral components, for example,*

_{m}

*D*_{l,m}= {

*D*

^{1}

_{l,m}, . . . ,

*D*

^{K}

_{l,m}}

^{T}. The matrix equations are solved with ease by direct inversion.

The direct inversion for each spectral component, however, is not feasible in the DFS method. Therefore, a different approach is needed in this case. The matrix 𝗕, real-valued and asymmetric, is decomposed into eigenvectors and eigenvalues (Golub and Van Loan 1996):

where 𝗨 is the eigenvector matrix of which each column is an eigenvector and **Γ** is a diagonal matrix of the eigenvalues. The eigenvectors are not orthogonal to one another but all eigenvalues are positive. The first mode is barotropic mode and the others are baroclinic modes that are wavier with the mode number (cf. Hoskins and Simmons 1975). In Fig. 1 the largest five eigenvalues and corresponding eigenvectors are illustrated. Except for the largest two, the eigenvalues are less than 0.1 and decrease sharply as the mode number increases. For example, the largest five are 0.4038, 0.1091, 0.0257, 0.007 470, and 0.002 842, while the smallest five are 0.4151 × 0^{−6}, 0.2371 × 10^{−6}, 0.1194 × 10^{−6}, 0.463 × 10^{−8}, and 0.9385 × 10^{−9}, respectively. Let **D** and **G** be presented with a linear combination of the eigenvectors,

and

where **U*** _{k}* is the

*k*th eigenvector and

*D*(

_{k}*G*) denotes the

_{k}*k*th vertical mode of the divergence (the forcing function);

*G*can be obtained by the use of the inverse matrix 𝗨

_{k}^{−1}. Multiplying the

*k*th column of 𝗨

^{−1}to (7) yields the following elliptic equation for the

*k*th vertical mode:

where *b _{k}* is the

*k*th eigenvalue. This is the same Helmholtz equation used for the semi-implicit time stepping of shallow-water equations but with different coefficient, which can be solved with efficiency using tridiagonal matrix inversion. If

*D*is inverted from (9), the divergence in (7) is obtained from (8a). From (9) the amplitude of vertical mode is reduced more effectively for the mode with larger eigenvalue.

_{k}To investigate the accuracy of DFS method, the elliptic equation is solved with an arbitrarily chosen forcing function *G*(*λ*, *θ*, *k*) = [cos*θ*/cos(*π*/6)] × exp(−*F*^{2}_{1} − *F*^{2}_{2} − *F*^{2}_{3}), where *F*_{1} = (*λ* − *π*)/Δ* _{λ}*,

*F*

_{2}= (

*θ*−

*π*/6)/(3

*π*/16), and

*F*

_{3}= (

*k*− 13)/3 with

*k*= 1, 2, . . . ,

*K*. The coefficient

*γ*is given as 0.01, which is nearly one order larger than the typical value used in the simulations. Such a large value was chosen to better demonstrate the effect of the inversion on the forcing function in (6). Matrix 𝗕 is calculated using the horizontal mean of the relaxation temperature of the Held–Suarez–Williamson simulation (see Fig. 11 later). Figure 2 shows the longitude–height cross section of the forcing function and the divergence field at

*θ*= 28.83° N for T84 (256 × 128 grids). The patterns are so close to each other that they cannot be visibly distinguishable. The vertical extent of the forcing function after inversion is reduced significantly, which indicates the significant reduction of amplitude of the barotropic component. In Fig. 3, the log–log plot of the

*L*

_{2}error norm of DFS, which is measured as the difference from the SHM-T106 over the 3D grid space, is shown for two cases of Δ

*= 3*

_{λ}*π*/16 (circles) and Δ

*= 3*

_{λ}*π*/8 (squares). The plots with squares, representing the case with larger zonal scale than the circles, show smaller error than the case of smaller zonal-scale forcing function. The error remains around

*O*(10

^{−4}), but it decreases with the resolution. The nearly straight lines indicate the exponential convergence of the solutions: The error is reduced by a factor of about 4.7 in case of Δ

*= 3*

_{λ}*π*/16 as the resolution increases from T42 to T106.

### b. Case II: Rossby–Haurwitz wave of zonal wavenumber 4

In this section, the dynamical core is tested with the Rossby–Haurwitz wave of zonal wavenumber 4 (Monaco and Williams 1975; Giraldo and Rosmond 2004). The model is run without the relaxation forcing as well as the Rayleigh damping, that is, *α* = *ν* = 0. The initial conditions for the divergence and vorticity are given as

where *c _{m}* = (

*m*+ 1)(

*m*+ 2),

*a*is the earth’s radius,

*m*= 4, and the first term (second term) in the parenthesis is the wave (zonal mean) component. As in Giraldo and Rosmond (2004), we used

*A*= 50 × (

*a*/

*m*) and

*B*= 0 in the above and following equations in this test. The temperature field is defined as

where

The surface pressure is initialized by

where

with

and with *S* = cos^{2}*θ*, *A _{R}* =

*Aa*

^{−2}, and

*m*=

_{s}*m*

^{2}+ 2

*m*+ 2.

In Fig. 4 the surface pressure, temperature, and meridional velocity in the lowest sigma level are compared for the 5-day integration with T106 resolution. The amplitude and zonal phase of the variables do not exhibit a significant difference between two dynamical cores. Detailed comparison of other variables indicated that both models have nearly identical results to each other. The results look very similar to those produced by other dynamical cores [e.g., see Figs. 12 and 13 of Giraldo and Rosmond (2004)]. The *L*_{2} norm error of the DFS model, which is measured from the SHM-T106 result, is presented in Fig. 5 for the three field variables shown in Fig. 4. The magnitude of the error varies with the variable: The temperature has the smallest error of the variables, which is smaller than the momentum (meridional velocity or equivalently the vorticity) by about three orders. It is found that the solutions converge with increasing model resolution, as is evident from the decrease of error by nearly two orders between T42 and T106. In this test case, the error decreases more sharply between T84 and T106 than other resolutions.

The initial vorticity field is given to be independent of the vertical level, but this barotropic structure may change with time because of the imbalance between the model variables. To see the vertical variation, we show the baroclinic component of the vorticity by taking the difference between *σ* = 0.3 and *σ* = 0.7 (Fig. 6) for T106 resolution. The baroclinic component, which is absent in the beginning, increased during the initial 5 days to reach the amplitude of about 10% of the initial vorticity (not shown).

### c. Case III: Jablonowski–Williamson (JW) test

This test consists of two cases: One is a test on the balanced state and the other is the life cycle of localized baroclinic waves that develop in the balanced initial state superimposed by a small-amplitude perturbation (Jablonowski and Williamson 2002). The model atmosphere is initially balanced by the following set of variables: The surface pressure is given as 1000 hPa in the global domain, and the velocity components are specified as

where *u*_{0} = 35 m s^{−1}, *σ _{υ}* = (

*σ*−

*σ*

_{0})

*π*/2, and

*σ*

_{0}= 0.252. The structure of the temperature is given by

where

where Γ now equals 0.005 K m^{−1} and *g* = 9.806 m s^{−2}. The surface geopotential is defined as

where

The model is time integrated for 30 days. The accuracy and stability of the model is evaluated by the surface pressure deviation from the initial value of 1000 hPa. In Fig. 7a the time variation of the *L*_{2} norm error of the surface pressure is presented for DFS-T62 and SHM-62. The error increases sharply during the first 2 or 3 days, but remains almost unchanged after that. The magnitudes of the errors for DFS and SHM are almost the same, showing only a few percent difference. Unlike the result of the explicit temporal scheme (e.g., Giraldo and Rosmond 2004), the surface pressure in the semi-implicit scheme, as in this study, exhibits nonoscillatory behavior. To see this behavior in more detail, we sampled the surface pressure at a grid point nearest to (*λ*, *θ*) = (0, 0). Time variation of the error (deviation from 1000 hPa) is illustrated in Fig. 7b and looks very similar to Fig. 7a, except for a small-amplitude oscillation that almost disappears around day 10. The *L*_{2}-norm surface pressure error of DFS model, which is measured as the difference from the SHM-T106 model, is shown in Fig. 7c and exhibits an exponential decrease with the resolution.

To test the life cycle of the baroclinic waves with the new dynamical core, the initial balanced state is added by the small-amplitude localized perturbation:

and

where *r _{p}* =

*a*/10 and (

*λ*,

_{c}*θ*) = (

_{c}*π*/9, 2

*π*/9). Since the basic state, which has vertical shear as well as a meridional temperature gradient, is baroclinically unstable, the perturbation grows exponentially with time until it is saturated and then finally breaks down (Jablonowski and Williamson 2002; Majewski et al. 2002; Giraldo and Rosmond 2004). Figure 8 shows the time evolution of the minimum surface pressure, an indicator of the wave development, for four resolutions of DFS and SHM. Note that the variation between the models is so small that the curves overlap to a large extent. Consistent with other dynamical cores (e.g., Giraldo and Rosmond 2004), the difference among simulation results grows with time but it becomes smaller again at days 13–14. A closer look at the solutions (lower panel of Fig. 8) indicates that there is a systematic difference among the results: The minimum surface pressure becomes large with the resolution and the SHM gives slightly larger values than DFS. The difference between DFS and SHM becomes smaller as the resolution increases. Figure 9 illustrates the vorticity field at

*σ*= 0.5 (the vertical layer

*k*= 13) by day 9. The developing baroclinic waves have finite zonal and meridional extent because the initial perturbation is localized. The approximate center of the baroclinic waves activity, which is initially given to be 20°E as in (17), has traveled to the local area around 120°W as a direct consequence of steering by the basic zonal-mean zonal flow (e.g., Simmons and Hoskins 1978). Four contour plots show very similar distributions of the vorticity, and it appears that the difference between different resolutions (e.g., DFS-T42 and DFS-T106) is more conspicuous than the difference between different dynamics (e.g., DFS-T106 and SHM-T106). Such a convergence property of the solution with the resolution is presented in Fig. 10, where the

*L*

_{2}norm error between DFS-model and SHM-T106 was evaluated in terms of the vorticity field as in Fig. 9. From T42 to T106 the error has decreased by nearly one order.

### d. Case IV: Short-term Held–Suarez–Williamson (HSW) simulation

The test cases in this and the next (case V) subsection consist of short- and long-term time integration with the Held–Suarez–Williamson forcing (Held and Suarez 1994; Williamson et al. 1998). The forcing is the Held–Suarez relaxation forcing modified by Williamson et al. (1998) to give a realistic temperature structure in the stratosphere and higher levels. These were designed to investigate the response of the dynamical cores to a realistic climatological forcing; hence, the model response is characterized by a more complicated disturbance field, unlike the previous test cases. Figure 11 presents the reference temperature and relaxation coefficient. The relaxation time scale is small in the lower Tropics and increases with both latitude and height. The Rayleigh damping for the vorticity and divergence of Held and Suarez (1994), a simple realization of momentum dissipation in the planetary boundary layer, is given only for layers *σ* ≥ 0.7 as (*σ* − 0.7)/0.3 per day.

To compare the results of short-term time integration, the initial condition is generated by the use of SHM-T42 as follows. The model is time integrated from the rest state with the temperature of the reference state in Fig. 11 and the surface pressure of 1013 hPa. By day 30 a small-amplitude random disturbance is added to the temperature, vorticity, and divergence fields to break the zonal-symmetric state. Since the zonal-symmetric state established by the symmetric relaxation forcing in the left panel of Fig. 11 is similar to the climatological state, and thus baroclinically unstable to the perturbation, the unstable baroclinic waves grow almost exponentially until their amplitude is saturated (e.g., Simmons and Hoskins 1978; Valdes and Hoskins 1988). The time needed for saturation depends on the initial amplitude of the perturbation because the time scale of the fast-growing mode is basic-flow dependent. The amplitude of the wave components is given between 10^{−8} and 10^{−9} in nondimensional units, which gives about 0.1 K for temperature and 10^{−6} Ω for vorticity and divergence. Beyond day 60, the model variables are well balanced among one another, and the global-averaged eddy kinetic energy settles down to a nearly steady state although it exhibits fluctuation with time. The wave-component data by day 200 are stored to be used as initial conditions for the short-term integrations. The SHM-T42 initial wave-component data can be used directly in the SHM with higher resolution. However, we need three further steps to use the SHM-T42 initial wave-component data in DFS models: The SHM-T42 data are transformed to the gridpoint data with a uniform longitude–latitude interval using the inverse Gauss–Legendre transform. Then, the gridpoint data are transformed to the wave components via double Fourier transform as in Cheong (2000a). Finally, the gridpoint data are transformed to the wave-component data for DFS-T42. This DFS-T42 initial wave data are used in the higher-resolution DFS models.

Figure 12 presents the initial fields of the surface pressure, geopotential height, and the vorticity at *σ* = 0.5 (left column) and the results of 9-day integration with SHM-T106 (right column). Unlike the test cases II and III, the field is characterized by globally distributed disturbances, but the activity is much more concentrated in the middle latitudes because of the strong temperature gradient therein. In Fig. 13 the same variables by day 9 as in the right column of Fig. 12 are shown for DFS-T42 and DFS-T106. Note that the fields are very smooth over the global domain without small-scale numerical noises. Basically the difference between DFS-T106 and SHM-T106 is not determinable from the contour plots, but it seems that the difference of DFS-T42 from both SHM-T106 and DFS-T106 is slightly larger than the difference between SHM-T106 and DFS-T106. Figure 14 provides the convergence of the results with the resolution, where the rmsd of the DFS model was measured from SHM-T106 result. The difference decreases with the resolution for both 5- and 9-day simulations, nearly by more than one order for the resolution increase from T42 to T106. For example, by day 9 the surface pressure (geopotential height) error is 1.5 hPa (12 m) for T42, but it is reduced to 0.06 hPa (0.5 m) for T106. The convergence of the solution with resolution and the large error for low resolution, as shown in this figure, is typical for both two- or three-dimensional models (e.g., Boer and Denis 1997; Williamson et al. 1998; Williamson 1999; Cheong 2000a).

### e. Case V: Long-term HSW simulation

In this section, the results of long-term simulations with HSW forcing are compared in terms of time-zonal- or time-wavenumber-averaged quantities. The model is integrated for 1400 days, of which the last 1200 days are used for time averaging. Figure 15 shows the time-zonal mean-zonal flow (top row) and temperature (bottom row) for DFS-T42, DFS-T106, and SHM-T106. The overall features such as middle-latitude tropospheric jets, equatorial surface easterlies, and stratospheric jets are well reproduced in all simulations although a detailed structure varies from one simulation to another. Common to other studies with HSW forcing, the equatorial easterlies are formed in all levels that have the smallest meridional extent in the middle troposphere. The meridional extent increases both upward and downward from that level and extends to about 45° in the highest sigma level of 0.02 (approximately 20 hPa) in both hemispheres. It is worthy to note that the maximum values for the tropospheric jets are comparable to each other in these simulations. However, the stratospheric easterly jet for the lowest resolution T42 is slightly different from (i.e., stronger than) those of T106 because of the sharp gradient and poor resolvability of it in the low-resolution model. Similar discrepancy between low- and high-resolution simulations is found in the stratospheric temperature field. If we consider the large temperature gradient and thermal wind balance in the upper levels, it is likely that the stratospheric jets will become stronger in the case where the uppermost level of the model is kept higher than the present one.

Figure 16 presents the standard deviation of temporal fluctuation of zonal velocity as a function of latitude and zonal wavenumber. Though not exactly the same, both results of the SHM and DFS models show very similar distribution to each other: The largest amplitude is found in middle latitudes while the smallest is at low latitudes. Note that the second peak in amplitude appears near the poles, as is consistent with the results of Held and Suarez (1994). As the zonal wavenumber increases, the midlatitude maximum is split into two local maxima, one poleward and the other equatorward. Common to both models of DFS and SHM, the hemispheric asymmetry is found for both T42 and T106 resolution even though the time averaging over 1200 days is carried out. This suggests to us that the temporal variation in this simulation with symmetric forcing includes a long-period component that is almost equal to the whole data span of 1200 days.

Time variation of the global-mean surface pressure is important in the context of the model’s performance on the mass conservation. If the logarithmic surface pressure equation is integrated in the vertical direction and is followed by multiplication of *p _{s}*, we have the surface pressure equation in the flux form

Therefore, the total mass (global integration of surface pressure) should be conserved in the dynamical cores:

When the surface pressure Eq. (18) is time integrated, the total mass conservation is guaranteed, as was illustrated in Rivier et al. (2002) for SHM and in Cheong (2000a) for DFS in terms of spectral components, whereas it is not the case for (1e). Therefore, it is worthy to check the total mass conservation in the dynamical cores that incorporate the logarithmic surface pressure instead of surface pressure as the prognostic mass-related variable to facilitate the semi-implicit time stepping.

Figure 17 shows the normalized error of total mass conservation for the 1200-day integration for the T42, T62, and T106 resolutions, where, for a clear comparison, the scale of the ordinate is given differently for three panels. The error grows almost linearly with time. Except for the lowest resolution T42, the SHM gives smaller errors than DFS, but the difference between SHM and DFS is not significantly large. It is interesting to note that the errors that appear in either positive or negative sign in each simulation do not change sign during the 1200-day integration, except for the initial period of a couple of months. In contrast to the cases with T42 and T62, the case of T106 maintained the negative sign on average. This behavior indicates the error is not of systematic bias as it should be. To see the convergence of the error with the resolution, a log–log plot of the absolute values of the error by day 1200 is presented in Fig. 18 for the DFS case results. The error decreases by more than one order between T42 and T106, implying a desirable spectral convergence. Not shown here, for a 10-yr integration with T42 the error of total mass is about 3.0 × 10^{−4} (approximately 0.3 hPa) by year 10 for both SHM and DFS models. This is smaller than the total mass error of an operational climate model nearly by one order [e.g., 2.5 hPa for 10-yr integration with T42 in Van den Dool and Saha (1993)]. The difference may result from the fact that the dynamical cores are so simple in comparison with the operational model, which does not include full physics, radiational forcing, and others.

### f. Computational efficiency of DFS

One of advantages of the DFS spectral method over the SHM is the computational efficiency, particularly for high-resolution models. The computational efficiency of the DFS spectral method can be found in Cheong (2000b) for the shallow-water model in terms of operation counts per time step. For the 3D primitive equation model as is used in this study, the sequential-computing time of DFS is reduced by 1.19, 1.28, 1.31, and 1.66 times for T42, T62, T84, and T106, respectively, as compared with the SHM model of the same resolution. This is a better performance estimated from the operation counts associated with the shallow-water equations without taking the data structure into consideration.

To measure more practical and rigorous computational efficiency of DFS, we used the high-performance computer NEC SX-5, a parallel-vector architecture, in the Korea Institute of Science and Technology Information (KISTI). The NEC SX-5 supercomputer in KISTI has eight CPUs with a theoretical peak performance of 80 gigaflops, and has the vector register length of 512 (see Kwon et al. 2004). For the model runs on this vector architecture, the source codes for DFS and SHM are optimally tuned to obtain maximized performance. More specifically, they include three important aspects: vectorization of do-loops, do-loop splitting for the matrix computation associated with the semi-implicit procedure, and two- or three-dimensional data movement for better vectorized computation. In Fig. 19 the computing time used to iterate 10 time steps of DFS and SHM models are compared for the horizontal resolutions from T106 to T1022. The relative efficiency of DFS over SHM (the ratio of computing time) is smaller than 2 for the lowest resolution in this figure, but increases significantly with the resolution. It needs less computing time than SHM by factors of about 2.09, 3.24, and 7.68 for T212, T682, and T1022, respectively. This efficiency is slightly lower than that evaluated by the operation count for the shallow-water model as in Cheong (2000b), which is attributed to the matrix computation associated with the semi-implicit time stepping.

Another advantage of DFS dynamical core over the SHM is that it does not need memory space for basis functions such as Legendre functions. The total memory required for the DFS dynamical core used in this study is 0.27, 0.99, 2.5, 4.1, 11.4, and 28.2 gigabytes (GB) for T106, T212, T340, T426, T682, and T1022, respectively. On the other hand, the required memory for SHM is 0.37, 1.42, 3.8, 8.1, 24.3, and 68.8 GB for the resolutions stated above. More than double the memory space is necessary for SHM in comparison with DFS in the case of the resolution of T212 or higher.

## 4. Discussion and conclusions

In this study a dynamical core for the atmospheric general circulation is developed using the DFS as basis functions. The model incorporates the hydrostatic approximation and sigma coordinate in the vertical direction without bottom topography. The model atmosphere is divided into 25 layers with uniform sigma depth. The DFS dynamical core is tested in terms of the semi-implicit method (i.e., 3D Helmholtz equation solver), balanced initial state, baroclinic waves, and short- and long-term simulations for four horizontal resolutions of T42, T62, T84, and T106. The DFS model has produced stable and accurate simulation results along with a desirable convergence with the resolution. The error that is measured from the highest-resolution SHM-T106 was found to depend on both the variables and test cases.

For all resolutions shown in this study, the DFS gives better computational efficiency than the SHM model, which reaches about 1.66 times at T106 resolution. As was revealed through the computation of highly optimized models on the high-performance computer NEC SX-5, this increased further as the resolution becomes better due to the availability of FFT for both longitudinal and latitudinal directions, increasing the efficiency by a factor of 7.68 for T1022 resolution. Along with the computational efficiency, the DFS needs less memory space in comparison with SHM because the basis functions are not necessary. The memory space for the DFS model is reduced to less than half of the memory required for SHM when the resolution is T212 or higher.

The advantages of the DFS dynamical core, such as a smaller memory requirement and better computational efficiency, are dependent on the horizontal resolution as well as the vertical resolution. As the vertical layers increase with the horizontal resolution fixed, the storage saving becomes less significant. At the same time, the computing time needed for the matrix operation associated with the semi-implicit procedure is proportional to the square of the number of vertical layers. Therefore, the computational efficiency stated above may be diminished slightly when the number of vertical layers increases. Further computations on the same high-performance computer with double the number of vertical layers (i.e., from 25 to 50 layers), in order to make this point clearer, demonstrated that the results were not significantly changed from those presented above. The computational efficiency of DFS relative to SHM was lowered by only about 2%–3% for the resolutions T212, T681, and T1022. The low sensitivity to the number of vertical layers is because the horizontal resolution is, in general, greater than the vertical resolution by one or two orders.

The DFS model uses the spherical (longitude–latitude) coordinate system without inviting any pole problems. Unlike the nonspherical coordinate systems, which are incorporated in new dynamical cores in an attempt to keep almost-even resolution on the sphere and avoid the pole problem, the spherical coordinate system makes the parallel (and also sequential) implementation of the numerical code easy and efficient because the data for the horizontal grid points can be processed with two-dimensional array in the model.

## Acknowledgments

This study was supported by the research and development project of Korea Meteorological Administration in the area of improvement of medium-range weather forecasts. This work was done while the author visited the department of Mechanical and Aerospace Engineering, University of California, San Diego, in 2004. Professor Paul F. Linden is acknowledged for generously providing the research facilities during his stay. Dr. Masao Kanamitsu at Scripps Institution of Oceanography is appreciated for constructive suggestions and useful comments for this study. The author thanks Dr. Francis Giraldo in the Naval Research Laboratory for helpful discussions on the test results in section 3. Anonymous reviewers are acknowledged for helping the author improve the manuscript. The author was partly supported by the financial aid from Pukyong National University, Korea, during the research year in 2004.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

^{−3}and k

^{−5/3}energy spectrum of atmospheric turbulence: Quasigeostrophic two-level model simulation.

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Hyeong-Bin Cheong, Dept. of Environmental Atmospheric Sciences, Pukyong National University, 599-1 Daeyeon-3-dong, Namgu, Pusan 608-737, South Korea. Email: hbcheong@pknu.ac.kr