. 84
( 136 .)


Schuster(1903) Triangular matrix multiplication transform
Dilts (1985) Modern reinvention of Schuster™s scheme
Brown (1985) Fast spherical harmonic transform
Orszag (1986) Fast transform for any basis satisfying
3-term recurrence relation
Elowitz&Hill Compares fast transforms of Dilts (1985),Brown (1985)
&Duvall (1989) but conclusion is false
Alpert&Rokhlin(1991) Fast Legendre-to-Chebyshev transform;
does not generalize to associated Legendre
Boyd (1992c) Showed FMM could be applied to non-Chebyshev basis
Dutt&Rokhlin (1993,1995) FMM nonequispaced FFT
Jakob-Chien&Alpert(1997) Grid-to-spherical-harmonics-to-grid
Yarvin&Rokhlin(1998) projective ¬lter for longitude/latitude
time-dependent double Fourier series
Foster&Worley(1997) Comparison of parallel spherical harmonic transforms
Healy&Rockmore& Fast spherical harmonic transform;
Kostelec&Moore(1999) freeware at www.cs.dartmouth.edu/ geelong/sphere
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

18.13 Equiareal Resolution and the Addition Theorem
De¬nition 40 (Equiareal Resolution) A numerical algorithm which has the property that its
numerical characteristics are invariant to a rotation of the north pole of the coordinate system so
that features of a given size are resolved equally well or badly regardless of whether they are located
at the poles, equator, or anywhere in between.

A so-called “triangular truncation” of a spherical harmonic basis has this property be-
cause of the following theorem and its corollaries.

Theorem 36 (ADDITION THEOREM:) Let (» , θ ) denote longitude and latitude as measured
relative to a set of coordinate axes rotated with respect to the original axes. Let (», θ) denote the
angles measured relative to the original, unrotated coordinate system. Then
m m
Yn (» ,θ)= Dmm (R) Yn (», θ)
m =’n

where the coef¬cients Dmm are functions of the rotation angle.
COROLLARY 1: The spherical harmonics of degree n form a (2n+1)-dimensional representation of
the continuous rotation group.
COROLLARY 2: A triangular truncation of spherical harmonics, that is, keeping only those har-
monics such that

n¤N all m (18.46)

gives equal resolution to equal areas on the globe, regardless of where those areas are located.

PROOF: Eq. (18.45) is a classical result discussed in most quantum mechanics texts such as
Merzbacher (1970).

The ¬rst corollary is merely a way of restating the theorem in the jargon of group theory.
When we rotate the pattern that is a particular spherical harmonic Yn , we create a new
function of latitude and longitude. Like all such functions, this may be expanded in a
spherical harmonics series, but there is no obvious reason why these series should not
contain an in¬nite number of terms. The theorem shows, however, that the series contains
at most (2n+1) terms, and all the non-zero harmonics have the same degree n as the function
that is rotated. Thus, the spherical harmonics of degree n form a closed subset under rotation
through arbitrary angles ” and this is what is required to form a “representation of the
rotation group”.
The collection of spherical harmonics that are retained in a “triangular truncation”
therefore are closed under rotation, too. This implies the property of equiareal resolution.
In the words of Orszag (1974), “The basic mathematical reason is that, under arbitrary ro-
tations, expansions truncated at harmonics of degree N remain truncated at degree N so
that the resolution of such series must be uniform over the sphere.”

18.14 Variable Resolution Spherical Harmonics Models
“Limited-area” weather forecasting models offer the advantages of very high resolution
over a small portion of the earth™s surface ” higher than would be affordable if this resolu-
tion were extended over the entire globe. Most national weather services run limited-area
models targeted at the country that pays for them. Limited-area models are also used to
track tropical hurricances, which have such small scales that it is dif¬cult for uniform res-
olution global models to track them accurately.
One obvious strategy is to employ a dense but uniform grid over a small portion of the
globe and specify in¬‚ow-out¬‚ow conditions at the sides, usually from climate data or from
interpolation of a global model. Such non-global limited-area codes are in wide use, but
there are dif¬culties. One is that specifying lateral boundary conditions turns out be very
hard; rather elaborate blending procedures are necessary to avoid corruption of the high-
resolution data within the domain by the low-resolution data at the boundaries. Moreover,
as models incorporate more and more physics and chemistry, it has become increasingly
painful to maintain two completely separate models, one global and one limited-area, us-
ing different numerical schemes, tuning parameters and so on.
An alternative is to use the global model as the limited-area model, too. This can be
done by a smooth change of coordinates that maps the surface of the sphere into itself. In
physical space, the transformed grid has a high density of points over the region of interest,
but decreases to lower and lower density as one moves away from the target region. No
arti¬cial sidewall boundary conditions are needed because there are no sidewalls. A single
model can be used for both global and regional forecasting by switching on or off a few
metric factors in the evaluation of derivatives.
Schmidt(1977, 1982) proposed a conformal sphere-to-sphere mapping which he tested
successfully in a simple code. With re¬nements, this has been adopted by M´ t´ o-France
for its primary operational weather forecasting code (“Arpege”) as described by Courtier
& Geleyn(1988), Courtier et al.(1991), and D´ qu´ &Piedelievre(1995). Hardiker(1997) has
shown that such mappings are equally effective for tracking hurricanes.
These variable-resolution models have been suf¬ciently successful that the desire to
combine global and regional models into a single code is not a major threat to continued use
of spherical harmonics for weather forecasting. Rather, the big Thing-That-Goes-Bump-in-
the-Night is the switch to massively parallel machines, which may or may not be happy
doing PMMTs for very large N .

Table 18.6: Variable Resolution Spherical Harmonic Models

References Comments
Schmidt(1977,1982) Conformal mapping to give high local resolution
to track cyclones
Courtier&Geleyn(1988) Variable resolution (mapped) weather prediction
Courtier et al.(1991) Arpege project at M´ t´ o-France:
global weather model with high resolution in Europe
D´ qu´ &Piedelievre(1995)
ee Experience with Arpege variable resolution weather model
Hardiker(1997) Conformal mapping for variable resolution

18.15 Spherical Harmonics and Physics
The spherical harmonics are the eigenfunctions of the two-dimensional Laplacian operator
Yn = ’n(n + 1) Yn
2 m m
‚2 ‚2
‚ 1
≡ 2 + cot(θ)
sin2 (θ) ‚»2
‚θ ‚θ
Because of this, the spherical harmonics are fundamental solutions to many problems in
In geophysics, for example, Haurwitz (1940) showed that the streamfunction for lin-
ear, barotropic Rossby waves was proportional to a spherical harmonic, that is to say, the
spherical harmonics are the quasi-geostrophic normal modes of the earth™s atmosphere.
Similarly, Longuet-Higgins (1968) has shown that for gravity waves in the barotropic limit,
the velocity potential is proportional to a spherical harmonic. In both cases, elementary
identities show that the velocities and other quantities may be written as pairs of spherical
Barrett (1958) observed that since the (purely westward) phase velocity is
cphase = ’ [Rossby-Haurwitz waves] (18.49)
n (n + 1)
one may synthesize a uniformly propagating disturbance from an arbitrary sum of spheri-
cal harmonics of the same degree n. In particular, an “eccentric” spherical harmonic, that is,
a spherical harmonic rotated so that its pole at some arbitrary latitude θp , is a steadily prop-
agating Rossby wave because of the addition theorem of Sec. 18.13. In more recent times,
such rotated spherical harmonics have been the basis for constructing nonlinear modons
in spherical geometry as explained by Tribbia(1984b).
One could multiply these examples with dozens from other ¬elds (Morse and Fesh-
bach, 1953). This close and intimate connection between the spherical harmonics and the
physics, as well as their good numerical properties, have helped to make spherical har-
monics popular.

18.16 Asymptotic Approximations I: Polar-Cap and Bessel
At high latitudes,
cos(θ) ≈ 1 sin(θ) ≈ θ (18.50)
& θ 1

and the horizontal Laplacian operator becomes
‚2 1 ‚2

+ +2 2
‚θ2 θ ‚θ θ ‚»
which is identical in form with the Laplacian in plane polar coordinates if we identify θ
with radius r and » with the polar angle. For a given zonal wavenumber m, the second
»-derivative in (18.51) becomes multiplication by (’m2 ) and the eigenequation 2 Yn =

’n(n + 1) becomes Bessel™s equation:
d2 m2
+ k’ 2
+ Jm (kθ) = 0
dθ2 θ dθ θ
k≡ (18.53)
n (n + 1)
A more heuristic (but equally correct) way of justifying this planar, polar coordinate
approximation is to simply look down on a globe from above one of the poles. As one
moves closer and closer to the poles, the sphere ¬‚attens into a plane, and the meridians
form a network of radial lines.
Since this “polar-cap” approximation becomes exact near the pole, the Bessel functions
must be consistent with the known behavior of the spherical harmonics at the pole. We
note that Jm (r) has an m-th order zero at r = 0 ” mimicing the m-th order zero of Yn at
the pole. If we expand Yn as a power series in θ, we ¬nd that the expansion contains only
every other power of θ; even powers of θ when m is even and odd powers of θ (because of
the “parity factor”, sin(θ)) when m is odd. Similarly, the expansion of Jm (r) is in alternating
powers of r; since the expansion begins with rm , the powers are the odd powers of r when
m is odd and the even powers of r when m is even.
The Bessel function (for m > 0) rises monotonically to a turning point at r ≈ m and
then oscillates for larger values of its argument. Fig. 18.9 illustrates three typical Bessel
functions. The turning colatitude is
θt ≈ [“turning colatitude”] (18.54)
n (n + 1)
When m is roughly equal to n (recall that n cannot be smaller than m), then the predicted
turning latitude is θt ≈ 1, that is, only about 30 degrees from the equator. The polar cap
approximation is not accurate that close to the equator, so what the Bessel approximation
tells us in this case is simply that the spherical harmonic has most of its amplitude at low
latitudes (where we shall use a different approximation given in the next section). Near the
pole, the Bessel approximation is still valid ” but the behavior of both Jm (kθ) and Yn is
dominated by the m-th order zero at the pole, so the polar cap approximation does not tell
us anything we did not already know.
When n m, however, the predicted turning point is close to the poles and the Bessel
approximation is more useful. There is some region equatorward of θt where the Bessel
approximation is legitimate, and we can use the asymptotic approximation
4 (2m+1) π
1 1
[n(n+1)] θ ∼ cos [n(n+1)] 2 θ ’ (18.55)
Jm θ θt
2 n(n+1) θ 2
π 4
Eq. (18.55) shows explicitly that a spherical harmonic oscillates equatorward of its turning
latitude. It also shows that the amplitude of the oscillation decreases in the direction of

the equator (as 1/ θ) ” the maximum of the harmonic is just equatorward of the turning
The “polar-plane” approximation, which invented by B. Haurwitz, has been systemat-
ically developed for geophysical applications (Bridger and Stevens, 1980).

Figure 18.9: Three Bessel functions, illustrating the behavior of spherical harmonics near
the pole, r = 0. Note that J10 (r), which is representative of high degree harmonics, has a
“turning latitude” at about r = 10: the function is exponentially decaying towards the pole
for small r, and oscillatory for larger r. (r ≡ n(n + 1) θ increases toward the equator.)

18.17 Asymptotic Approximations, II: High Zonal Wavenum-
ber & Hermite Functions
Abramowitz & Stegun (1965) give the asymptotic approximation

Pm+n [cos(θ)] ∼ q exp(’ mφ2 ) Hn ( m φ)


. 84
( 136 .)