. 82
( 136 .)


or ¬nite elements grids are popular is that x is roughly the same everywhere so that the
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

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

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

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
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.
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

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

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.

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
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
Ck (µ) Cj (µ) (1 ’ µ2 )m dµ
m m
if j = k (18.37)
= 0

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

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

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
max |Ck (x)| ∼ O k m+1/2 (18.38)
x∈[’1, 1]
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

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

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 )].

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
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


. 82
( 136 .)