points is large. The rub is that the grid nevertheless is non-uniform and this leaves two

problems: (i) programming is much more dif¬cult than for a uniform grid and (ii) the ac-

curacy of ¬nite differences deteriorates because the grid is not entirely regular. In addition,

recent models have been bedeviled by stability problems and spurious excitation of lon-

gitudinal wavenumber ¬ve disturbances. Heikes and Randall(1995a,b) used a “twisted”

icosahedral grid which respects the parity symmetry with respect to the equator, but at

the expense of more complication. Algorithms that replace the icosahedron with another

Platonic solid, such as a cube or octagon, have been tried, too (Table 18.8).

All of these variants work, but all are messy and most have required some tweaking.

These faults of almost-uniform grid methods are some consolation for the intricacies of

spherical harmonics.

18.26 Spectral Alternatives to Spherical Harmonics, I: Robert

Functions

The Fast Fourier Transform (FFT) converts grid point values of a function to its Fourier

coef¬cients in only O(N log2 N ) operations. Since the spherical harmonics are an ordinary

Fourier series in longitude, we can apply the FFT in », but not to the Gegenbauer polyno-

mials in θ. This inspired serious efforts in the late 60™s and early 70™s to ¬nd substitutes for

the latitudinal basis functions.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

434

Andr´ Robert suggested the basis functions

e

sin|m| (θ) cos(nθ) [“modi¬ed Robert functions”] (18.89)

These basis functions reproduce the polar behavior of the spherical harmonics. The high

order zeros for large wavenumber m eliminate the “pole problem”; one can integrate using

(18.89) with a fairly large time step and no special tricks. Merilees (1973a) pointed out

the fatal ¬‚aw of the “Robert” functions: they are horribly ill-conditioned. The problem is

that for large m, the sin|m| (θ) factor is zero except in a narrow boundary layer around the

equator. Within this layer

n2 θ2

cos(nθ) ≈ 1 ’ ≈1 if |θ| (18.90)

1

2

so that the lowest few Robert functions are graphically indistinguishable for large m. This

translates into large round-off errors when a function is expanded in terms of these (al-

most) linearly dependent functions. The result is that the Robert functions have never

been employed for large-scale numerical models.

18.27 Spectral Alternatives to Spherical Harmonics, II: Parity-

Modi¬ed Fourier Series

Merilees (1973b) and Orszag (1974) have discussed parity-modi¬ed Fourier series in which

the latitudinal dependence of the zonal wavenumber m component of a scalar function is

expanded in the form

∞

s

(18.91)

fm (θ) = sin (θ) an cos(nθ)

n=0

where

m even

0

s≡ (18.92)

n odd

1

The pseudospectral grid points are

π (2i ’ 1)

θi ≡ (18.93)

i = 1, . . . , N,

2N

that is, an evenly spaced grid with all points on the interior of θ ∈ [0, π]. This implies

that we may freely divide out sins (θ) and then apply the FFT (a Fast Cosine Transform,

actually). The coordinate singularities of the differential operators at θ = 0 and π cause no

major problems because we do not need to evaluate any quantities at the poles.

One obvious concern is that the spherical harmonics (and all well-behaved functions on

the sphere) have m-th order zeros at the poles whereas the modi¬ed Fourier basis functions

in (18.91) have at most a ¬rst order root. However, to quote Orszag (1974): “The detailed

behavior expressed by the sin|m| (θ) factor is usually of little direct interest in a numerical

simulation . . . The above argument suggests that surface harmonic-series contain much

information on the behavior of functions near the poles that is not of primary interest.”

Boyd (1978b) offers an amusing illustration of Orszag™s point by applying the pseudospec-

tral method with 60 basis functions to calculate the eigenvalues of the two-dimensional

18.28. PROJECTIVE FILTERING FOR LATITUDINAL FOURIER SERIES 435

Table 18.9: Spectral Alternatives to Spherical Harmonics: Bibliography

Note: “Double Fourier” denotes a Fourier basis in both latitude and longitude.

References Comments

Robert (1966) “Robert” function basis, but these are very ill-conditioned

Ill-conditioning of “Robert” basis sinm (θ) cos(nθ)

Merilees (1973a)

Merilees (1974) Double Fourier series with high latitude ¬ltering

Orszag (1974) Comparison of spherical harmonics with double Fourier series

Boyd (1978b) Comparison of Chebyshev, Fourier & Associated Legendre

for boundary value and eigenvalue problems

Yee(1981) Poisson equation on sphere; double Fourier basis

Spotz&Taylor high-latitude-¬ltered double Fourier series

& Swarztrauber(1998)

Shen(1999) double Fourier basis

Cheong(2000, 2001) double Fourier series

Laplacian, i. e. the spherical harmonics themselves, for m = 49. Although the 60 term ap-

proximation may have at most 61 zeros (counting multiplicity and the two roots of sin(θ))

and the spherical harmonics have 98 + (n ’ m) roots, the lowest twenty eigenvalues are

calculated to within an accuracy of 1 part in 70 million.

For so-called “jury” problems, that is to say, for boundary value and eigenvalue prob-

lems, Boyd (1978b) argues that the parity-modi¬ed Fourier series is the best basis set. It is

much easier to program than spherical harmonics and the accuracy does not seem poorer.

For time-dependent or “marching” problems, however, the “pole” problem returns in

full force. Orszag (1974) has shown that the time-step limits for realistic calculations will

be quite severe. He suggested several remedies. First, one may apply a polar ¬ltering ev-

ery ten or twenty time steps. Second, one may split the convective terms into (i) a surface

harmonic of degree 2 that gives the ¬‚ow over the poles and (ii) the rest of the convective

terms and then treat (i) implicitly; the bulk of the terms, i. e. (ii), may be treated explic-

itly without instability because (ii) does not involve ¬‚ow over the pole. Third, treating

viscous terms implicitly ” almost always necessary with spectral and pseudospectral al-

gorithms ” tends to stabilize the convective terms, too, at least for a Backwards-Euler or

other implicitly-dissipative scheme.

Even so, longitude/latitude double Fourier series have been little used for time-dependent

calculations in the quarter-century after Orszag™s article. Part of the reason is inertia; as ex-

plained earlier, the relatively slow cost of associated Legendre transforms has not been

a major liability for spectral models even at the highest resolution employed in weather

forecasting at the dawn of the third millenium, T 639.

In recent years, however, double Fourier series have made a comeback in the articles

of Shen(1999), Cheong(2000, 2001) and Spotz, Taylor and Swarztrauber(1998). Spotz et al.,

for example, revived an idea of Merilees(1974) with an improved ¬lter, and showed that

latitudinal Fourier series can be both fast and accurate when the polar ¬ltering is strong

enough.

18.28 Grid-to-Spherical-Harmonics-to-Grid Projective Filter

for Time-Dependent Latitudinal Fourier Series

The ultimate polar ¬lter is to transform grid point values on a tensor product latitude/longitude

grid to a triangular truncation of spherical harmonics and then sum the spectral series to

obtain a different set of grid point values. The output of the fore-and-back transform will

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

436

inherit the equiareal resolution property of spherical harmonics. The pole problem is com-

pletely suppressed.

This combination of transforms has been described as both a “spherical ¬lter” and a

“projection” in the literature; we shall call it a “projective ¬lter”. The reason for the termi-

nological confusion is that the projective ¬lter is different from other ¬lters in that its pur-

pose is not to provide arti¬cial viscosity or suppress small-scale noise and aliasing instabil-

ity. Rather, the grid-to-spectral half of the projective ¬lter is a non-dissipative transform. It

is, in the usual parlance of functional analysis, a “projection” of the function f (», θ) onto the

subspace spanned by the spherical harmonics. The spectral-to-grid half of the transform

is merely a summation of the truncated series. The grid-to-grid transform is more than a

projection (because it is a summation, too), but in some sense less than a ¬lter (because the

only ¬ltering is a truncation).

The weakness of a latitude/longitude Fourier method is that projective ¬ltering has a

cost which is the same order of magnitude as an Associated Legendre transform. However,

because only the unknowns themselves must be projectively ¬ltered ” no derivatives need

be transformed ” double Fourier algorithms have the potential of extending the lifetime

of global spectral methods on the sphere for at least another decade.

The method which Spotz and Swarztrauber (2000) label both “naive” and “standard”

is to compute both halves of the projective ¬lter through separate Parity Matrix Multipli-

cation Transforms (PMMT) in latitude. (All transformations in longitude are performed

by the Fast Fourier Transform (FFT), which convert grid point values into functions fm (θ)

which describe how the coef¬cient of cos(m») or sin(m») varies with colatitude.) The ad-

jective “Parity” means that grid point values of f (θ) are ¬rst split into components that

are symmetric and antisymmetric with respect to the equator θ = π/2 through f S ≡

{ f (θ) + f (π ’ θ) } and f A ≡ { f (θ) ’ f (π ’ θ) }. Each component is then projectively

¬ltered separately; the even and odd parity ¬ltrates are recombined by adding their grid

point values in the northern hemisphere and subtracting the ¬ltered values of f A from the

¬ltrate of f S in the southern hemisphere.

Yarvin and Rokhlin(1998) noted that one could reduce the cost by 25 % by combining

the forward and backward transform matrices into a single N — N matrix for small m

while performing the forward and backward transforms separately for large m. (Note

that the forward and backward transform matrices are rectangular with smaller dimension

(N ’ m); for the extreme of m = N , the transform matrices are only vectors; two vector-

vector multiplies are obviously cheaper than an N — N matrix-vector multiply.)

Swarztrauber and Spotz (2000) developed the “Weighted Orthogonal Complement”

projective ¬lter which has only half the cost of the “naive”/“standard” method and only

two-thirds the expense of Yarvin and Rokhlin™s “direct” method. Better yet, they prove

that the WOC matrices can be written as subsets of a couple of “master” matrices, reduc-

ing storage for transform/projection matrices from O(N 3 ) to only O(N 2 ).

For the shallow water wave equations, nine Legendre transforms are needed at each

time step with a spherical harmonic basis. With a double Fourier expansion and the WOC

projective ¬lter, the transform cost is reduced to the equivalent of only three Legendre

transforms.

Jakob-Chien and Alpert(1997) developed a projective ¬lter based on Fast Multipole

Method (FMM) ideas; this was improved by Yarvin and Rokhlin(1998). As with other

multipole transforms, the FFM projective ¬lter is faster than the WOC or any MMT-type

method in the limit N ’ ∞, but the cross-over point where the FMM scheme is faster is

likely to be at very large N .

Spotz and Swarztrauber(2000) carefully compare all these projective ¬lters for various

N , both with and without the application of the Two-Thirds Rule. (Dealiasing truncation

reduces the advantages of the newer MMT schemes relative to the “naive” method.) They

18.29. SPECTRAL ELEMENTS ON THE SPHERE 437

Table 18.10: Projective Filters Bibliography

References Comments

Jakob-Chien&Alpert(1997) Grid-to-spherical-harmonics-to-grid

Yarvin&Rokhlin(1998) projective ¬lter for longitude/latitude

time-dependent double Fourier series

Swarztrauber&Spotz(2000) “Weighted Orthogonal Complement” algorithm reduces

Spotz&Swarztrauber(2000) storage by O(N ); faster than alternatives

because algorithm mostly stays in the on-chip cache

Cheong(2001) Approximate projective ¬lter (hyperviscosity)

which combines damping with projection

¬nd that their WOC algorithm is faster than the “direct” and “standard” PMMT and the

FMM at least up to N = 200. The WOC method is moderately faster than formal operation

counts would suggest; the very low storage of the WOC leads to more ef¬cient utilization

of the on-chip cache memory of current machines. The WOC method also parallelizes very

well.

Cheong(2001) employs an alternative: Polar ¬ltering by a combination of fourth order

m

and sixth order hyperviscosity. Recall that Yn is an eigenfunction of the two-dimensional

Laplace operator on the surface of the sphere with the eigenvalue »n = ’n(n+1). It follows

that a dissipation which is proportional to a high power of the Laplacian can be tuned to

strongly damp all spherical harmonics with n > N while having only a slight effect on

harmonics retained within a triangular truncation. The Fourier-Galerkin discretization of

the Laplace operator is a tridiagonal matrix, so the ¬lter is very fast.

In contrast to the Spotz-Swarztrauber and FMM methods, Cheong™s hyperviscosity is

not an exact projection onto spherical harmonics. Since most time-dependent models re-

quire a little scale-dependent dissipation for noise-suppression anyway, an approximate

projective ¬lter is satisfactory. However, there has not been a detailed comparison of the

relative speed and accuracy of exact projective ¬lters versus approximate ¬lters that damp

and project simultaneously.

Swarztrauber and Spotz (2000) note, “The goal of an O(N 2 log N ) harmonic spectral

method remains elusive; however, the “projection” method provides a new avenue of re-

search. Perhaps the development of a fast projection [projective ¬lter] will prove to be

easier than the development of a fast harmonic transform.” Suf¬cient progress has al-

ready been made to encourage the development of double Fourier methods on the sphere.

However, as of 2000, latitude/longitude Fourier series have not yet been applied to an

operational forecasting or climate model.

18.29 Spectral Elements on the Sphere

Another alternative to global spectral methods is domain decomposition. A moderately

high order Legendre series is used in each subdomain, and then the pieces are stitched

together as in ¬nite elements to obtain a global approximation. The great virtue of spec-

tral elements is that each subdomain can be assigned to a different processor of a multi-

processor machine. The articles listed in Table 18.11 show that the algorithm parallelizes

very well, and seems promising for both atmospheric and oceanic modelling.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

438

Table 18.11: Spectral Elements on the Sphere Bibliography

References Comments

Ma(1992,1993a,b) Ocean basin models

Iskandarani&Haidvogel global ocean model

&Boyd(1995)

Levin&Iskandarani ocean model

&Haidvogel(1997)

Taylor&Tribbia & atmospheric model;

Iskandarani(1997) shallow water equations

Curchitser&Iskandarani Ocean modelling on a variety of

& Haidvogel(1998), massively parallel computers

Haidvogel&Curchitser

&Iskandarani &Hughes

&Taylor(1997)