## Abstract

The authors describe a recent development and some applications of a spectral element dynamical core. The improvements and development include the following: (i) the code was converted from FORTRAN 77 to FORTRAN 90; (ii) the dynamical core was extended to the generalized terrain-following, or hybrid *η*, vertical coordinates; (iii) a fourth-order Runge–Kutta (RK4) method for time integration was implemented; (iv) moisture effects were added in the dynamical system and a semi-Lagrangian method for moisture transport was implemented; and (v) the improved dynamical core was coupled with the Community Atmosphere Model version 2 (CAM2) physical parameterizations and Community Land Model version 2 (CLM2) in such a way that it can be used as an alternative dynamical core in CAM2. This spectral element version of CAM2 is denoted as CAM-SEM. A mass fixer as used in the Eulerian version of CAM2 (CAM-EUL) is also implemented in CAM-SEM. Results from multiyear simulations with CAM-SEM (coupled with CLM2) with climatology SST are also presented and compared with simulations from CAM-EUL. Close resemblances are shown in simulations from CAM-SEM and CAM-EUL. The authors found that contrary to what is suggested by some other studies, the high-order Lagrangian interpolation (with a limiter) using the spectral element basis functions may not be suitable for moisture and other strongly varying fields such as cloud and precipitation.

## 1. Introduction

The spectral element method (SEM) combines the high-order accuracy of spectral methods and the gridding flexibility of the finite element method. Although SEM has been used in the ocean and atmospheric modeling by several research groups (e.g., Ma 1993; Iskandarani et al. 1995; Haidvogel et al. 1997; Taylor et al. 1997; Giraldo and Rosemond 2004; S. Thomas et al. 2006, personal communication), the atmospheric modeling based on SEM coupled with full physics (including topography and the land surface model) is just beginning to be explored (cf. Wang et al. 2004). This article presents some details of the recent development of a spectral element dynamical core. We also present simulations from the spectral element version of the Community Atmosphere Model version 2 (CAM2) (CAM-SEM) and compare them with simulations from the Eulerian version of CAM2 (CAM-EUL).

In atmospheric modeling, SEM is usually used for the horizontal discretization, while the conventional finite difference method is used in the vertical discretization. Fournier et al. (2004) presented the idealized Held–Suarez forcing experiment with a spectral element dynamical core with a *σ*-coordinate in the vertical. This spectral element dynamical core was based on the previous work of Taylor et al. (1997) on global shallow water modeling with SEM. In the current study, the dynamical core is extended to the generalized terrain-following, or hybrid *η*, vertical coordinate. The improvements and further development of the spectral element dynamical core include the following: (i) the code was converted from FORTRAN 77 to FORTRAN 90; (ii) the dynamical core was extended to the generalized terrain-following, or hybrid *η*, vertical coordinates; (iii) a fourth-order Runge–Kutta (RK4) method for time integration was implemented; (iv) the moisture effects were added in the dynamical system and a semi-Lagrangian method for moisture transport was implemented; and finally (v) the spectral element dynamical core was coupled with the CAM2 physics and Community Land Model version 2 (CLM2) in such a way that it can be used as an alternative dynamical core in CAM2. This spectral element version of CAM2 will be denoted as CAM-SEM. A mass fixer as used in CAM-EUL is also implemented in CAM-SEM. Some details of these developments will be presented in the following section.

The purposes of this paper are 1) to describe the recent developments of the spectral element dynamical core and basic features of CAM-SEM and 2) to present the simulated states from multiyear simulations using CAM-SEM and compare them with simulations using CAM-EUL. This paper will focus on the analysis of multiyear simulations using CAM-SEM (including CLM2) with climatological sea surface temperature (SST). Other simulations, such as aqua-planet experiments, (Atmospheric Model Intercomparison Project) AMIP-type simulations, and simulations with local mesh refinement (LMR) will be represented in separated papers. Some preliminary results were presented in Wang et al. (2004) and Baer et al. (2006).

## 2. The spectral element dynamical core

In this section, we present some details on the recent development of the spectral element dynamical core: (i) the *η*-version of the spectral element dynamical core, (ii) the RK4 method for time integration, (iii) the semi-Lagrangian tracer transport, and (iv) coupling with the CAM2 physical parameterizations and CLM2.

### a. The η version of the spectral element dynamical core

To facilitate the coupling with CAM2 physical parameterizations, the code is first converted from FORTRAN 77 to FORTRAN 90, and FORTRAN 90 modules are used. It is also necessary to extend the vertical coordinate from *σ*-coordinate to the generalized terrain-following vertical coordinate, or *η*-coordinate, as used in CAM2. The vertical discretization is the same as the one used in CAM2 (Collins et al. 2003). Including the effects of moisture, the hydrostatic atmospheric equations for moist air in the hybrid *η*-coordinates can be written in vector form as follows:

where *σ* = *p*/*π*, *π* denotes the surface pressure, *E* = ½(**v** · **v**), *c**_{p} = {1 + [(*c _{pυ}*/

*c*) − 1]

_{p}*q*}

*c*, and

_{p}*T*= {1 + [(

_{υ}*R*/

_{υ}*R*) − 1]

*q*}

*T*. It should be noted that the wind vector

**v**refers to the horizontal wind components, while the gradient and divergence operators are computed along the constant

*η*surfaces.

Note that to include the effects of moisture, the virtual temperature *T _{υ}* is used in place of temperature

*T*,

*only*in the pressure gradient term −

*RT*

_{υ}**∇**ln

*p*in the momentum Eq. (2.1), in the

*ω*term (

*R*/

*c**

_{p})

*T*(

_{υ}*ω*/

*p*) in the thermodynamic equation [(2.2)], and in the hydrostatic balance equation [(2.5)].

Calculation of the vertical advection terms may be written as η̇(∂*ψ*/∂*η*) = η̇(∂*p*/∂*η*)(∂*ψ*/∂*p*). The vertical discretization scheme is the same as used in CAM2 (Collins et al. 2003). The calculations of the horizontal (or lateral) pressure gradient terms on the constant *η*-surfaces are done as follows:

With the definition of *p*(*η*) = *A*(*η*)*p*_{0} + *B*(*η*)*π*, we have (1/*p*)∂*p*/∂*π* = *B*/*p*, where A(*η*) and B(*η*) are two constants defined in such a way that a set of levels *η _{k}* may be defined such that

*η*=

_{k}*A*+

_{k}*B*, and

_{k}*p*

_{0}= 10

^{5}Pa,

*p*(1) = 354.46 Pa. The upper and lower boundary conditions are η̇(

*p*

_{t},

*π*) =

*ω*(

*p*) = 0 and η̇(

_{t}*π*,

*π*) = 0, respectively [Collins et al. 2003, his Eqs. (3.2) and (3.1)].

The basic validation simulation with the Held–Suarez idealized forcing (Held and Suarez 1994) using the *η* version of the spectral element dynamical core was performed and compared with the simulation using CAM-EUL (Wang et al. 2004; Baer et al. 2006). In the simulations presented in this paper, the RK4 method was used for time integration. The RK4 method is described in the following subsection.

### b. The RK4 method for time integration

To provide efficient time integration as well as flexible coupling between dynamics and physical parameterizations, a fourth-order Runge–Kutta method is implemented in the spectral element dynamical core, replacing the previous leapfrog scheme. For a system of equations of the form

where *ϕ* and *ψ* are prescribed at *t* = *t*_{0}, the RK4 scheme can be written as follows (Hildebrand 1987, p. 291):

with

and

where *h* is the time step size used in dynamics. Because of the explicit and one-step fashion of the RK4 method, subcycling fractional steps can be conveniently used for dynamics so that a large time step can be used for physics. For example, in the simulations presented in this paper, one large time step Δ*t* = 2400 s is used for physics while 8 substeps with time step *h* = Δ*t*/8 are used for dynamics. Note that the tendencies of *u*, *υ*, and *T* from physical parameterizations are calculated/updated only at large time steps and are kept unchanged during the dynamical substeps. The specific humidity *q* is updated from physical parameterizations at large time steps and is also kept unchanged during the dynamical substeps. The transport of *q* is done through the semi-Lagrangian scheme as described in the next subsection.

### c. The semi-Lagrangian tracer transport

Just as in the global spectral method, the spectral element method suffers from its difficulty in maintaining monotonicity and its tendency to produce negative moisture. One way to alleviate the problem is to use the semi-Lagrangian transport scheme with monotonicity or maximum/minimum limiter imposed.

The equation for the specific humidity *q* of water vapor and mixing ratio of other constituents is written in the *η* coordinate, excluding sources and sinks:

The contributions from source terms (i.e., physical parameterizations) have already been added to the moisture fields at the time level (*n*). As in CAM-EUL, the transport is divided into horizontal and vertical transport substeps (Collins et al. 2003). In the semi-Lagrangian form, this can be written in general as

The operator *L _{λ}*,

*in (2.13) represents the horizontal interpolation (of*

_{ϕ}*q*) to the departure point assuming the vertical velocity η̇ = 0, while the operator

^{n}*L*in (2.14) represents the vertical interpolation (of

_{η}*q*

^{n}^{+1}

*) to the departure point assuming the horizontal velocity*

_{h}**v**= 0. It should be noted that in the semi-Lagrangian schemes commonly used in global models, the fully three-dimensional trajectory calculations are usually used, and the source terms are usually extrapolated and integrated along the trajectories (e.g., Collins et al. 2003; Hortal 1999).

The trajectories are computed in the spherical coordinates except near the polar regions where, for |*ϕ*| > 70°, the trajectories are computed on the cube surfaces (for mapping onto a cube, see Taylor et al. 1997). A second-order Runge–Kutta (RK2) method, referred to as the *modified Euler method* or the *improved polygon method* (Lambert 1973, p. 118), is used for trajectory calculations, that is,

where *u* and *υ* are normalized by the radius of the earth. Therefore, only *one interpolation* for the median value of velocity at time level (*n* + ½), *u ^{n}*

^{+1/2}(

*λ*,

_{M}*ϕ*) and

_{M}*υ*

^{n}^{+1/2}(

*λ*,

_{M}*ϕ*) is needed in evaluating the trajectory and calculating the departure points (

_{M}*λ*,

_{d}*ϕ*), and only

_{d}*one search and one interpolation*are needed in determining the specific humidity at time level (

*n*+ 1), or

*q*

^{n}^{+1}

*=*

_{h}*q*(

^{n}*λ*,

_{D}*ϕ*). Because of the fractional time steps used in the dynamical module, the velocity fields at time level

_{D}*n*are available; and the above algorithms provide easy and convenient calculations for the semi-Lagrangian transport scheme.

The vertical transport is the same as used in CAM-EUL (i.e., the Hermitian cubic interpolation is used with a sufficient condition for monotonicity imposed) except that the RK2 scheme is used for the trajectory calculations as in the horizontal trajectory calculations.

In CAM-SEM, the horizontal transport is done on the spectral element nodes. We initially used the spectral element basis functions for horizontal interpolations. That is, a seventh-order Lagrangian interpolation on the Gauss–Lobatto–Legendre (GLL) nodes (see appendix for the details) was used. In addition, a maximum/minimum limiter is imposed on the moisture field such that the interpolated value of moisture at the departure point is within the minimum and maximum of the moisture field in the element of the departure point. This limiter is the same as that used in Giraldo et al. (2003).

However, we found that when the spectral element basis functions are used as the interpolants (with a limiter), the moisture field produced may be too noisy and that many negative moisture points are produced out of the moist convection parameterization scheme. We therefore tested the bilinear interpolation and monotonic Hermite cubic interpolation for moisture fields, while the spectral element basis functions are still used for interpolations of wind fields. We found that the negative moisture from the moist convection parameterization scheme never occurred when using bilinear interpolation and rarely occurred when Hermite cubic interpolation was used.

The monotonic Hermite cubic interpolation is similar to the one used in CAM-EUL (Collins et al. 2003; Williamson and Rasch 1989). The sufficient condition for monotonicity with degrees Celsius continuity is imposed as in CAM-EUL. The one-sided derivatives are used near the element boundaries. The pure advection test of a cosine bell around the globe (e.g., section 4 in Williamson and Rosinski 2000) is used to test the above schemes. Figure 1 presents normalized *l*_{2} and *l*_{∞} errors of the schemes using spectral element interpolation and monotonic Hermite cubic interpolation. With the latter interpolation scheme, the normalized *l*_{2} error is just slightly above 0.5 and the normalized *l*_{∞} error is about 0.4–0.5, similar to the results of the monotonic Hermite cubic interpolations of Williamson and Rosinski (2000, their Figs. 5a, b). The results are slightly worse than that presented in Williamson and Rosinski (2000). One reason may be that the one-sided derivatives are used along the elemental boundaries here. With the spectral element interpolation scheme, the normalized *l*_{2} error is between 0.1 and 0.2 and the normalized *l*_{∞} error is about 0.2. It should be noted that although the tests of pure advection of a cosine bell around the globe showed reasonable results when the spectral element interpolation scheme was used, tests of advection of the moisture field showed that the monotonic Hermite cubic interpolation scheme is the only acceptable scheme among the schemes tested. This also indicates the limited usefulness of the idealized advection of the cosine bell, at least in the context of numerical experiments with this version of the model. All the simulations presented in this paper use the monotonic Hermite cubic interpolation.

Finally, a mass fixer as used in CAM-EUL (Collins et al. 2003) is also implemented in CAM-SEM to prevent the drift of total dry mass and to remedy the nonconservative property of the semi-Lagrangian moisture transport scheme. As in CAM-EUL, the maximum 10% threshold of correction is imposed, though it is never materialized. The effects of the mass fixer can be seen in the comparison of surface pressure distributions with or without the mass fixer (Fig. 3 in section 3b below).

### d. Coupling with the CAM2 physical parameterizations and CLM2

The spectral element dynamical core with features described above is coupled with the CAM2 physical parameterizations and CLM2 in such a way that it can be used as an alternative dynamical core in CAM2. CAM2 version cam2_0_dev15 and CLM2 version clm2_deva_08 were used for this purpose.

To prepare the initial and boundary conditions for CAM-SEM, the initial conditions and boundary conditions such as sea surface temperatures and sea ice fractions of T42 CAM-EUL were used to interpolate bilinearly to the CAM-SEM grids. The same 26 levels of hybrid *η*-coordinates are used in the vertical.

The processing of surface data from the so-called *raw data* on latitude–longitude grids to CLM2 surface data on spectral element grids follows from the procedures used in CLM2, except that bilinear or *distance-weighted* interpolation is used in place of *area-weighted* interpolation. The surface data on CAM-SEM grids obtained in this way were compared with the T42 CLM2 surface data to check the validity of the data. In Table 1, we also present the global average of a few variables from surface data on CAM-SEM grids and on CAM-EUL Gaussian grids.

Before we look at the long-term simulations, we examine the timing of the simulation using CAM-SEM and compare it with the timing of CAM-EUL. First, we estimate the model resolution for CAM-SEM as follows. We calculate the area of each element *A _{e}.* The average element length

*L*is estimated as the square root of

_{e}*A*, or

_{e}*L*

_{e}=

*A*

_{e}. The estimates are given in Table 2. We use the notation E96N64 to represent CAM-SEM grids with total 96 elements (4 × 4 = 16 elements on each cube face) and 64 nodes in each element. The average node spacing is about 320 km with a minimum of about 150 km for CAM-SEM grids E96N64. The standard T42 Eulerian version of CAM2 has an estimated resolution of about 320 km or coarser (see discussion in Laprise 1992).

To assess the timing, we run CAM-SEM (E96N64) and CAM-EUL (T42) for 10 days with various numbers of processors on IBM machine bluesky at the National Center for Atmospheric Research (NCAR). The time metrics are the wall-clock time used in stepon (the time stepping subroutine of CAM2) for a 10-day simulation, according to the standard timing library diagnostic in CAM2. These tests are for pure Message Passing Interface (MPI), and the tests where the number of tasks is one were performed with Single Program, Multiple Data (SPMD) (MPI) turned off. The total wall-clock time output from stepon is shown in Table 3. The same data are also used to calculate the scaling performance shown in Fig. 2. The timing performance of CAM-SEM, as it stands now, is comparable with CAM-EUL. It is expected that with further improvement on implementation and performance tuning, its performance can be improved further. This is only one set of tests run on the IBM platform and run with a pure MPI configuration. On a different platform, the results may be different. But it provides a basic benchmark comparison. It should also be noted that although the pure spectral element dynamical core scales almost linearly (cf. Fournier et al. 2004), CAM-SEM scales only asymptotically. One of the reasons may be due to the performance of coupling procedures between dynamical core and physics and between physics and CLM2. It would be worthwhile to examine the coupling performance in detail in the future.

## 3. Numerical simulations with CAM-SEM

In this section, we analyze CAM-SEM simulations with climatological SSTs and compare with CAM-EUL simulations as well as observations. Both CAM-SEM (E150N64) and CAM-EUL (T42) were run for more than 5 yr, initialized from 1 September of model year 0. The 3-yr simulations from model years 2–4 are taken for the analysis presented in this paper.

For experiments with CAM-SEM, there is no adjustment or tuning of parameters in the physical parameterizations, except that the threshold relative humidity value for high clouds is increased from 90% in CAM-EUL to 95% in CAM-SEM. It should be noted that the time step for CAM-SEM (E150N64) is Δ*t* = 2400 s while the time step for CAM-EUL (T42) is *δt* = 1200 s. However, the leapfrog time integration scheme is used in CAM-EUL, so the time step for physical parameterization calculations is 2*δt* = 2400 s, the same as that used in CAM-SEM. It should also be noted that the difference is small between CAM-SEM simulations using Δ*t* = 2400 s and Δ*t* = 1200 s.

To assess the effects of the mass fixer on model simulations, we have also run CAM-SEM with the mass fixer switch off, branched off from the 1-yr simulation of the above CAM-SEM control run. We compare the three–model year simulation from year 2 to 4 with that of the control run.

The standard CAM2 model output to history files in NetCDF format on CAM-SEM (element, node) grids is used in CAM-SEM. It is interpolated from CAM-SEM [element, node] grids to CAM-EUL T42 ([latitude, longitude) Gaussian grids. The high-order (seventh order) polynomials of SEM basis functions are used in the interpolations, except for precipitation and other two-dimensional fields where the bilinear interpolations are used.

### a. Global averages

The global annual average properties of the simulations are given in Table 4. The uncertainty range (error bar) of the satellite data is about ±5 W m^{−}^{2} (Wild and Roeckner 2006). The differences between CAM-SEM and CAM-EUL are all within the range of the certainty error bar.

To assess the effects of interpolation, the global means computed directly on the native SEM grids are also included in Table 4. The difference between the global averages calculated on the native SEM grids and those calculated from the interpolated ones are negligibly small. All the analyses below will use the interpolated fields, unless otherwise noted.

Table 5 compares the global annual average properties of simulations from CAM-SEM with and without mass fixer. All quantities are computed from data on the native SEM grids. The global average surface pressure (PS) and the sea level pressure (PSL) are also shown in the table. It seems that the largest effects of the mass fixer are on the surface energy budgets. This is understandable in terms of the large moisture content of the lower atmosphere. Another factor that may have contributed to the difference is that the initial conditions of the land surface model are not from any previous runs, thus not in a model balanced state. Initial conditions for CLM2 are set to the default, “non-spunup values” in the model. They are identical for both CAM-SEM and CAM-EUL (the control run).

### b. Zonal means

Figure 3 shows the zonally averaged annual mean PSL from CAM-SEM (with or without the mass fixer), CAM-EUL, and the National Centers for Environmental Prediction (NCEP) analysis. The difference between CAM-SEM with the mass fixer and the one without is very small. Large differences exist in the two polar regions between CAM-SEM, CAM-EUL, and the NCEP analysis. One reason for these large differences may be the different topography used in each case. However, a detailed analysis of this important issue of simulation in the polar regions in the context of the spectral element model is deferred for a future study.

The zonally averaged annual mean precipitation is shown in Fig. 4. The difference between CAM-SEM and CAM-EUL is very small. Both show similar biases, in particular the so-called double ITCZ problem.

The zonally averaged annual mean temperature is shown in Fig. 5. Although the cold bias in the tropical upper troposphere is reduced in CAM-SEM, the cold bias in the polar regions from tropopause up is increased in CAM-SEM. The increased latitudinal temperature gradient is also shown in the positive bias of the zonally averaged annual mean zonal wind of the regions (Fig. 6). It should also be noted that it takes longer time (about 10 yr) integrations in order to have a reliable simulation of polar jets (P. Rasch 2006, personal communication).

The cold bias in the polar regions from the tropopause up in CAM-SEM is correlated with the positive bias of high cloud fraction in these regions. The zonally averaged annual mean cloud fraction of CAM-SEM and CAM-EUL is shown in Fig. 7. The largest difference is the high cloud fraction in both polar regions, especially in the North Pole region. This may be due to the cold bias in the two polar regions, especially in North Pole region (Fig. 5).

### c. Geographic distributions

The annual mean sea level pressure from CAM-SEM and CAM-EUL is shown in Fig. 8. The main features are quite similar to each other, with the possible exception of the related surface zonal winds in the Southern Hemisphere, where we found that CAM-SEM winds are actually in better agreement with NCEP reanalyses (not shown).

The annual mean precipitation is shown in Fig. 9. The two model simulations are in good agreement with each other. Large differences exist in the regions of high topography and in some tropical islands, which may by partially due to the somewhat different topography used in each simulation.

## 4. Summary and conclusions

We have successfully coupled a spectral element dynamical core with the CAM2 physical parameterizations, implemented in such a way that it can be used as an alternative dynamical core in CAM2.

Some of the new developments of the spectral element dynamical core include the following: (i) the code was converted from FORTRAN 77 to FORTRAN 90; (ii) the dynamical core was extended to the generalized terrain-following, or hybrid *η*, vertical coordinates; (iii) a fourth-order Runge–Kutta (RK4) method for time integration was implemented; (iv) moisture effects were added in the dynamical system and a semi-Lagrangian method for moisture transport was implemented; and finally (v) the improved dynamical core was coupled with the CAM2 physical parameterizations and CLM2 in such a way that it can be used as an alternative dynamical core in CAM2.

The spectral element version of CAM2 (CAM-SEM) has been integrated for more than 5 yr with climatological SSTs and the results are compared with those of the Eulerian version of CAM2 (CAM-EUL) and observations. The close similarity between CAM-SEM and CAM-EUL is shown in the results.

Like the global spectral method, the spectral element method has serious deficiencies in the advection of moisture fields that usually have large horizontal gradients. A semi-Lagrangian scheme for moisture transport is implemented in CAM-SEM. We found that when high-order spectral element basis functions are used to interpolate moisture in the semi-Lagrangian scheme, negative moisture occurs frequently from the convective parameterization scheme, even though the pure advection test of cosine bell has produced better results (smaller errors). Even though this semi-Lagrangian transport problem may be model specific, monotonic cubic Hermite interpolation is a better choice for moisture interpolation, both from theoretical considerations (shape preserving) and numerical experiments (no negative moisture from convective scheme). Therefore, the pure advection test case should be used to check for gross error only.

The CAM-SEM model timing performance is similar to that of CAM-EUL. Improvement may be expected from refinement of the implementation of the semi-Lagrangian transport scheme in CAM-SEM. Searching in the semi-Lagrangian transport scheme in CAM-SEM starts from the departure point of the previous step (or the base element) and then the elements connecting this base element until the departure point is found. It may be desirable to have a preferential searching *direction*, but this is not currently implemented. Another area that may affect the model timing performance of CAM-SEM is the time integration schemes. It would be worthwhile to test different time integration schemes, such as the semi-implicit scheme.

## Acknowledgments

We would like to extend our appreciation to Dr. David Williamson (NCAR) for several beneficial discussions on model development and implementations. Thanks also go to Jerry Olson of CGD/NCAR for discussion on the implementation of the semi-Lagrangian scheme and mass fixer in CAM and many staff members of CGD/NCAR for sharing their knowledge of CAM2 and CLM2. Comments on the manuscript by Dr. Phil Rasch (NCAR) are also greatly appreciated. The research described herein was supported in part by Grant DEFG0201ER63247 to the University of Maryland from the DOE/BER/CCPP program. Computations were performed on computers at NCAR. Some computations during early model development and tests were also done on computers at NERSC and ORNL.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

### APPENDIX

#### The Spectral Element Methods

The spectral element method combines the grid flexibility of the finite element method and the high-order accuracy and efficiency of the spectral method. A nice introduction to SEM can be found in Deville et al. (2002, their chapter 2, section 2.4). The extension to the multidimensional problems using tensor-product forms is discussed in Deville et al. (2002, their chapter 4). A brief introduction to the basics of SEM based on the Gauss–Lobatto–Legendre integration is given in this appendix.

To apply the spectral element method to the global atmospheric modeling, it is convenient to map a sphere to a cube (Taylor et al. 1997). The standard spectral element method is applied on atmospheric equations on the cubic surface for the horizontal discretization. In the vertical discretization, conventional finite difference is usually used.

The interpolation in SEM is based on the Lagrangian interpolation. We will review the Lagrangian interpolation first. The GLL quadrature will then be reviewed, followed by a brief description of the Boyd–Vandeven (BV) filter, which are usually used in spectral element modeling.

##### The Lagrangian interpolation

The interpolation in SEM is based on the Lagrangian interpolation. Following Lanczos (1997, 5–6), the Lagrangian interpolation by polynomials can be briefly summarized as follows: To construct a polynomial of the order *N* that coincides with *u*(*x _{k}*) at the

*N arbitrarily given*points

*x*for

_{k}*k*= 0, 1, . . . ,

*N*, we first construct the

*fundamental polynomial F*

_{N}_{+1}(

*x*) by multiplying all the root factors,

Dividing by the root factors *x* − *x _{k}*, we construct the

*N*+ 1

*auxiliary polynomials*,

and

with *k* = 0, 1, . . . , *N*. The auxiliary polynomials *π _{k}*(

*x*) have the property that they are zero at all the root points

*x*, except at

_{i}*x*where they have the value 1,

_{k}where *δ _{k}*,

*is the*

_{i}*Kronecker delta function*. If we now form the sum

we obtain a polynomial that has the following properties: its order is *N* and it assumes the values at the prescribed points *x _{k}* for

*k*= 0, 1, . . . ,

*N*. This is called

*Lagrange’s interpolation formula*.

##### The GLL quadrature

The spectral element method can be regarded as a special type of high-order, collocation finite-element method. SEM based on Legendre polynomials makes use of the Gauss-Lobatto-Legendre quadrature. The GLL quadrature nodes are also the collocation points.

For completeness, we derive the GLL quadrature as follows: for the 1D problems, the (*N* + 1) GLL points Ξ_{N}_{+1} = (*ξ*_{0}, *ξ*_{1,} . . . , *ξ _{N}*) are solutions to the equation (Deville et al. 2002, p. 63; Davis and Rabinowitz 1984, p. 104):

(which includes two end points ±1) where the *L*′* _{N}*(

*ξ*) is the derivative of the Legendre polynomial of degree

*N*,

*L*(

_{N}*ξ*), which is the eigensolution of the differential equation

The Lagrange polynomial with nodal points at the GLL points can be derived as follows: Let the fundamental polynomial *F _{N}*

_{+1}(

*x*) be

The Lagrange polynomials at the GLL quadrature points are

Therefore, the Lagrange polynomial with nodal points at the GLL points can be written as

These are the basis functions used in the spectral element method. The Lagrangian interpolation is given by

Obviously,

(i.e., the GLL quadrature points are also the collocation points). The GLL quadrature formula is given by

and the GLL quadrature weights are given by (Deville et al. 2002, p. 451)

Therefore, the spectral element method can be considered a high-order collocation finite element method.

##### The spectral filter

The Boyd–Vandeven filter of order *N* has the form (Taylor et al. 1997)

where erfc(*x*) = 1 − erf(*x*) is the complementary error function. The BV filter was also used in a spectral element ocean model (Levin et al. 1997).

## Footnotes

*Corresponding author address:* Dr. Houjun Wang, NCAR, P.O. Box 3000, Boulder, CO 80307. Email: houjun.wang@noaa.gov