time-step limit is not unduly expensive.

On a sphere, the simplest tactic is to use latitude and longitude as coordinates and apply

an even grid spacing in » and θ. Denoting the radius of the sphere by a, the rub is that the

distance x between two grid points on a circle of colatitude θ is

(18.32)

x = a sin(θ) »

which tends to zero at the poles as shown graphically in Fig. 18.7.

Figure 18.7: A grid with uniform spacing in latitude and longitude as viewed from Polaris,

the Pole Star. The meridians (lines of constant longitude) all converge at the pole. Conse-

quently, x ’ 0 as the pole is approached, even though » is constant. The very small

grid spacing near the pole is the “pole problem”: One must use a very, very short time

step, or the numerical ¬‚ow will violate the Courant-Friedrichs-Levy criterion and become

unstable.

The result is that a time step which is stable and reasonable for tropical and temper-

ate latitudes ” typically 10 minutes for a fully explicit, “primitive equations” numerical

weather prediction model ” will give instability near the poles. As the run continues, this

localized instability will gradually spread equatorward like a cancer until the whole globe

is consumed with numerical error.

One remedy is an arti¬cially high viscosity near the poles. Heavy damping destroys

the accuracy of the model at high latitudes, but it is known that high latitudes have little

physical in¬‚uence on lower latitudes, and few complaints about inaccurate local forecasts

have come from polar bears.

The alternative is to attack the root of the problem, which is that, close to the poles, the

resolution in longitude is needlessly high. Why use very small x near the poles when the

accuracy of the model as a whole is constrained by the much poorer resolution everywhere

18.11. SPHERICAL HARMONICS: INTRODUCTION 399

else? The high zonal wavenumber components of the solution are the only ones that be-

come unstable. Therefore, the preferred methods for solving the “pole problem” are based

on lowering zonal resolution near the poles.

One approach, used in the GFDL-Princeton General Circulation Mode, is to selectively

delete grid points from the latitude-longitude grid at high latitudes. This is messy to pro-

gram. Worse still, the grid becomes irregular and the resulting uncentered ¬nite difference

approximations are not very accurate.

A second alternative is a zonal-scale-selective ¬ltering at high latitudes. This is merely a

variant of the arti¬cial polar damping technique, of course, but a scale-selective dissipation

is physically reasonable in a way that a scale-independent ¬ltering is not. The reason is

that the model cannot resolve the very short longitudinal scales of the damped components

except near the poles. There is no reason to mourn the damping of zonal wavenumber 30 at

85 degrees of latitude where its associated zonal scale is only 20 km when the longitudinal

resolution at the equator is as coarse as 200 km.

The third alternative is to use a spectral or pseudospectral method with spherical har-

monics as the basis set. Because of the close connection between the spherical harmonics

and the geometry of the sphere, this basis set plays no favorites with respect to latitude: the

spherical harmonics have the property of “equiareal resolution” on the sphere as explained

in Sec. 18.13 below.

To be sure, the spherical harmonics are more complicated to program than simple ¬nite

differences on a uniform latitude-longitude grid. But if one wants to survive the cancerous

instabilities of the “pole problem”, one is forced to use something more complicated than a

uniform » - θ grid.

18.11 Spherical Harmonics: Introduction

Spherical harmonics ¬x the pole problem, contain the proper “parity factors”, and offer

exponential convergence for functions that are in¬nitely differentiable on the sphere. Un-

fortunately, spherical harmonics are more complicated than any of the other basis sets in

common use because they are two-dimensional. We will later discuss icosahedral grids and

near-uniform triangularizations to emphasize that this dif¬culty is unavoidable: there is no

way to solve problems on the sphere that does not require some complicated programming

and careful attention to detail.

The spherical harmonics are

m, n non-negative integers

Yn (», θ) ≡ eim» Pn (θ)

m m

(18.33)

such that n ≥ m

Yn (», θ) ≡ e’im» Pn (θ) (’1)m

’m m

Thus, in longitude, the spherical harmonics are merely a Fourier series. Eq. (18.33) dis-

plays the Fourier series written in the form of complex exponentials; some authors prefer

ordinary sines and cosines in », but the difference is only one of convention.

m

The functions Pn (θ), the “associated Legendre functions” of order [zonal wavenum-

ber] m and degree n, are where all the complications lie. They are de¬ned by

Pn (θ) ≡ sinm (θ) Cn’m (cos[θ])

m m

(18.34)

where the Cn (cos[θ]) are the “Gegenbauer polynomials”3 . The subscript denotes the de-

m

gree of the polynomial; the coef¬cients of the polynomials are different for each different

zonal wavenumber m.

3 Warning: The superscript m used in this section is related to the superscript m of Appendix A via m=m -1/2.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

400

The factor sinm (θ) in (18.34) has a three-fold signi¬cance. First, recall that sin2k (θ) ≡

(1’cos2 [θ])k = a ¬nite Fourier cosine series with 2k terms. Thus, all the associated Legendre

functions with m even may be written in terms of a Fourier cosine series. However, when

m is odd, one factor of sin(θ) cannot be converted into a cosine polynomial, and [combining

m

sin(θ) with the cosine polynomial via the usual trigonometric identities] Pn (θ) is a (¬nite)

Fourier sine series. Thus, sinm (θ) is the “parity factor” discussed earlier in the chapter; it

is automatically built into the spherical harmonics.

The second signi¬cance of the sine factor is that it implies that the spherical harmonics

have an m-th order zero at the poles ” a multiple root in colatitude whose degree is de-

termined by the wavenumber in longitude. It is this high-order zero which enables spher-

ical harmonics to avoid the pole problem: the basis functions of high zonal wavenumber

m would have wastefully high longitudinal resolution at high latitudes except that the

sinm (θ) factor forces these spherical harmonics to have negligible amplitude near the poles.

The third consequence of the sinm (θ) factor is that modes with large m and (n’m) m

have little amplitude outside a narrow band centered on the equator. This justi¬es the

Hermite function asymptotic approximation discussed later.

The de¬nition of the “degree” n in (18.34) seems rather peculiar in that the degree of the

Gegenbauer polynomial is (n ’ m); it would seem much more natural to de¬ne the degree

of the spherical harmonic to equal that of the polynomial. However, the convention shown

in (18.34) has become universally adopted because this de¬nition of degree-of-harmonic

simpli¬es almost all other descriptions of the associated Legendre functions.

The most important of these is truncation: what cutoffs on m and n are most ef¬cient?

As will be explained below, one can prove that the spherical harmonics give equal resolu-

tion of all areas of the globe if and only if the truncation is triangular, that is, the retained

basis functions satisfy the inequalities

|m| ¤ N n¤N

and “triangular truncation” (18.35)

for some integer N as illustrated in Table 18.4. A triangular truncation up to and including

wavenumber N is usually denoted by the shorthand of “T” plus the numerical value of N .

A “TN ” truncation retains a total of (N + 1)2 spherical harmonics.

The Gegenbauer polynomials satisfy an orthogonality relation that can be expressed in

two illuminating ways:

π

Ck (cos[θ]) Cj (cos[θ]) sin2m (θ) dθ

m m

if j = k (18.36)

= 0

0

1

Ck (µ) Cj (µ) (1 ’ µ2 )m dµ

m m

if j = k (18.37)

= 0

’1

First, note that both forms of orthogonality apply only to different polynomials of the same

zonal wavenumber m. With the spherical harmonics, one is not dealing with a single set of

orthogonal polynomials, but rather with a countable in¬nity of sets.

Second, (18.36) shows a fourth role for the factor of sinm (θ) in (18.34): when two asso-

ciated Legendre functions of the same m are multiplied together, they collectively supply

the sin2m (θ) weight factor in (18.36). Thus, the associated Legendre functions of the same

order m are orthogonal on θ ∈ [0, π] with a weight function of unity.

When written in terms of µ = cos(θ), the Gegenbauer polynomials are ordinary (rather

than trigonometric) polynomials. Eq. (18.37) shows that the Gegenbauer family includes

the Chebyshev and Legendre polynomials as the special cases m = ’1/2 and m = 0. Like

other orthogonal polynomials, the Gegenbauer polynomials and their derivatives may be

18.11. SPHERICAL HARMONICS: INTRODUCTION 401

m

Table 18.4: An illustration of allowed combinations (m, n) for the spherical harmonics Yn .

The X™s show the modes that retained in a triangular truncation with a cutoff at N = 3

(“T3”); the O™s show the modes that are added when the cutoff is increased to N = 5

(“T5”). The total number of modes is (N + 1)2 in a TN truncation.

m’ -5 -4 -3 -2 -1 0 1 2 3 4 5

n

“

0 X

1 X X X

2 X X X X X

3 X X X X X X X

4 O O O O O O O O O

5 O O O O O O O O O O O

evaluated for arbitrary m and n via three-term recurrence relations. Handbooks on mathe-

matical functions are stuffed with such identities, so we shall here focus on the qualitative

properties of the basis functions.

One important property is that as m increases, the roots of the n-th degree polynomial

move closer and closer to the origin (equator). Since basis functions give good resolution

where they are oscillating and poor resolution where they vary slowly and monotonically,

it follows that the high order Gegenbauer polynomials give good resolution near the equa-

tor and poor resolution near the poles ” exactly what is needed to avoid the pole problem

and give uniform resolution over the whole globe. [When we consider a triangular trunca-

tion of harmonics, the small m functions resolve the small polar cap region; the large areas

near the equator are resolved through the combined efforts of all the harmonics.]

Recall that an N -term orthogonal polynomial series minimizes the mean square of the

error where the error is weighted by the weighting function in the orthogonality integral.

The weight function (1 ’ µ2 )m tolerates large error near the poles for the sake of higher

accuracy near the origin.

Another consequence of this stress on the interior is that when the Gegenbauer polyno-

mials are normalized so that the integrals in (18.36) and (18.37) are unity when j = k, one

¬nds that

˜m

max |Ck (x)| ∼ O k m+1/2 (18.38)

x∈[’1, 1]

m

where the tilde over Ck denotes the normalized Gegenbauer polynomials. In other words,

Eq. 18.36 implies that the maximum values of a normalized polynomial are large compared

to its O(1) average value on the interval x ∈ [’1, 1]. This implies that the amplitude of the

polynomial™s oscillations is highly nonuniform (in contrast to the Chebyshev polynomials,

whose maxima and minima are all the same). The amplitude increases rapidly as one ap-

proaches either pole. However, these oscillations of the polynomials are tempered by the

decay of the sinm (θ) factor as we move away from the equator so that the associated Leg-

endre function rather resembles a Hermite function ” that is, to say, there is a band around

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

402

Figure 18.8: Schematic of typical spherical harmonics. The “zonal” harmonics have zero

zonal wavenumber and their latitudinal structure is that of the ordinary Legendre polyno-

mials Pn (cos θ) where θ is colatitude. The “sectoral” harmonics are the modes which have

no latitude circles as nodal lines. The “sectoral” harmonics of high zonal wavenumber

are equatorially trapped and have negligible amplitude outside the tropics. The “tesseral”

harmonics are the general case.

the equator where it oscillates with many (roughly) equal crests and troughs; outside this

band, Pn (cos[θ]) simply decays monotonically to zero. (This band is very narrow if n ≈ m

m

1 but extends almost all the way to poles if (n ’ m) is not small in comparison to

and m

m.) For the special case of the so-called “zonal harmonics”, m = 0, the band of oscillations

is global and the high latitude “decay zones” do not exist. Fig. 18.8 schematically illustrates

typical spherical harmonics.

18.12 Legendre Transforms and Other Sorrows

18.12.1 FFT in Longitude/MMT in Latitude

In longitude, spherical harmonics are just a Fourier series. The optimum longitudinal grid

is therefore evenly spaced. Since two-dimensional transforms are most ef¬ciently per-

formed as a nested sequence of one-dimensional partial sums or transforms, irregardless of

how these one-dimensional sums are evaluated, the longitudinal partial sums are invari-

ably evaluated via the FFT. In latitude, however, it is necessary to use Gaussian quadrature,

which is slower [O(N 2 ) versus O(N log2 N )].

18.12. LEGENDRE TRANSFORMS AND OTHER SORROWS 403

The optimum Gaussian quadrature abscissas for a given wavenumber are the roots of

an m-th order Gegenbauer polynomial. However, because we must calculate transforms

in both » and θ, it is necessary to use the same points in colatitude for all m. Thus, the

spherical grid is the direct product of an evenly spaced Fourier grid in longitude with the

quadrature abscissas for Gauss-Legendre integration. (These grid points in colatitude are

the roots of a Legendre polynomial, and must be recalculated whenever the resolution is

changed.)

Orszag (1974) makes the comment “ . . . the goal of numerical simulations is accurate

reproduction of physical hydrodynamics; it is not clear, at least to the author, that quadratic

conservation properties [and other means of suppressing aliasing at the expense of addi-

tional operation] have very much to do with the attainment of that goal.” Chen(1993) has

performed numerical experiments suggesting that dealiasing is simply a waste of com-

puter resources if the wavenumbers near N are appropriately ¬ltered. However, his con-

clusion is controversial in the sense that models in Eulerian coordinates require such strong

¬ltering near the aliasing limit that the upper third of the wavenumber spectrum is largely

destroyed as it would be in a dealiased scheme. Most spherical harmonics models, espe-

cially those using Eulerian rather than semi-Lagrangian advective schemes, are de-aliased.

However, the semi-Lagrangian treatment of advection has allowed explicit dealiasing

to be dropped from the current ECMWF operational forecasting model. The nonlinear ad-

vective terms are never directly evaluated in a semi-Lagrangian method; instead, ¬‚uid is

advected through the method of characteristics and the resulting ¬elds are interpolated