. 85
( 136 .)



m ’ ∞, n ¬xed, φ ∼ O

where q is a constant, φ is latitude [not colatitude] and Hn is the usual Hermite polynomial.
This approximation is complementary to that of the previous section. First, the “polar cap”
approximation is accurate only for high latitudes while (18.56) is limited to a band around
the equator. Second, the Bessel approximation becomes more useful when n m whereas
the Hermite formula is most accurate when m n , that is, in the limit of large zonal
wavenumber for ¬xed n where n is the number of zeros of the harmonic (excluding those
at the poles).
With a triangular truncation at degree N , there are but two harmonics with zonal
wavenumber N , and their common latitudinal structure is, according to (18.56), given

by the Gaussian, exp(’ N φ2 ). Now the highest zonal wavenumber would impose the

Figure 18.10: A comparison of the exact associated Legendre function, P9 (φ) [solid] with
its asymtotic approximation, the Hermite function exp(’[9/2] φ2 ) [dashed]. φ is colatitude
& φ = π/2 is the equator.

strictest requirement on the time step if we expanded the unknowns in a double cosine
series (instead of spherical harmonics)“ but zonal wavenumber N does not translate into
a very small longitudinal grid interval near the poles for the spherical harmonics because
the Gaussian has an exponentially small amplitude near the pole. The equatorial con¬ne-
ment embodied in (18.56) is how the “pole problem” is avoided when we use spherical
harmonics to integrate time-dependent equations.
Figs. 18.10 and 18.11 show that the Hermite approximation is not uniform in degree
n (or n ); as n increases for ¬xed zonal wavenumber, the harmonic becomes wider and
wider in latitude and (18.56) becomes less and less accurate. However, inspecting Fig. 18.11
more carefully, we note that even as the approximation becomes numerically less accurate,
it remains qualitatively correct in the sense that the harmonic oscillates close to the equator
and then decays monotonically towards the poles on the far side of the turning latitudes.√
[One may in fact obtain a very accurate approximation to P13 merely by replacing m in
(18.56) by 3.33 to narrow the Hermite function.]
One obvious question is: what is the physical reason for this equatorial con¬nement of
spherical harmonics with (n ’ m) m? Fig. 18.12 shows a ray-tracing argument that
makes the trapping at least plausible. Because of the spherical geometry, a ray that leaves
the equator at say a 45 degree angle (that is, moving in a north-east direction) will not reach
the poles if it moves in a straight line on the sphere, that is, if the ray moves as a geodesic.
When the zonal wavenumber is large in comparison to the latitudinal wavenumber, i. e.
(n ’ m), then the equivalent ray is almost parallel to the equator and the turning
latitude is only a short distance from the equator.
To be sure, this argument is only heuristic. We have not attempted to justify ray-tracing.
Indeed, the spherical harmonics often appear in electrostatic and gravitational problems
where no rays are evident. However, this refraction-by-geometry is a real and physically
important effect in many wave problems. Boyd (1985b) gives a geophysical discussion.

Figure 18.11: A comparison of the exact P13 (φ) [solid] with its asymptotic approximation,
exp[’(9/2)φ2 ] H4 (3φ), which is shown as the dashed curve. For a given zonal wavenumber
[m = 9 for both this ¬gure and its predecessor], the approximation worsens with increasing
degree, but remains qualitatively accurate even for fairly large degree.

18.18 Software: Spherical Harmonics
John Adams and Paul N. Swarztrauber of the National Center for Atmospheric Research
have written a very comprehensive software library to compute spherical harmonics, their
derivatives, transforms from grid point values to spectral coef¬cients, interpolation from
an evenly spaced latitudinal grid to the Gaussian grid employed by forecasting models,
and computation and manipulation of vector spherical harmonics. SPHEREPACK 3.0
is now available at http://www.scd.ucar.edu/css/software/spherepack. A shallow water
model, complete with a full suite of test cases, and a three-dimensional global circulation
model are also available from NCAR, complete with documenation.

Figure 18.12: Schematic of the behavior of rays trapped on the surface of a sphere. In
the absence of other refractive mechanisms, the rays follow geodesics, that is, great cir-
cle trajectories. The illustrated ray cannot propagate poleward of the dotted circles: it is
latitudinally con¬ned.

18.19 Semi-Implicit Spherical Harmonic Methods for the Shal-
low Water Wave Equations
Real forecasting models have three spatial dimensions, but since the radial (vertical) coor-
dinate is irrelevant to spherical harmonics, we will discuss the simpler case of the shallow
water wave equations, which are two-dimensional in space, based on the classic paper of
Bourke (1972).
The shallow water equations may be written in the form

= ’(ζ + f ) k — V ’ V·V (18.57a)
Vt φ+
=’ · (φ V ) ’ ¦ D (18.57b)

where ζ is the vorticity [strictly, the vertical component of the vorticity vector], f is the
Coriolis parameter, k is a unit vector in the vertical direction, V is the horizontal velocity
vector, and D is the divergence of V . The height ¬eld has been split into a mean part,
¦, and deviations from the mean, φ , because these parts will be treated differently in the
time-marching algorithm.
The major problem in Eq. (18.57) is that its free oscillations are composed of two sets of
waves with wildly disparate time scales: (i) slow, low-frequency Rossby motions and (ii)
high frequency gravity waves. Existing data networks cannot adequately resolve the grav-
ity waves, and these high frequency motions have little effect on weather or climate. From
the perspective of global weather modelling, the gravity waves are simply noise. Unfor-
tunately, this noise shortens the maximum stable time step for an explicit time differencing
scheme from about 60 minutes to 10 minutes.
For this reason, so-called “semi-implicit” time marching schemes (Chapter 12) have
become very popular in meteorology. A small set of terms which are essential to the (linear)
propagation of gravity waves are treated implicitly, but most of the linear terms and all the
nonlinear terms are handled explicitly. This makes it possible to increase the time step by
a factor of six with no signi¬cant loss of accuracy and without the iterations that would be
needed in a fully implicit code.
One major reason why spherical harmonics have carried the ¬eld in meteorology is that
these expansions allow very simple and ef¬cient semi-implicit time-stepping algorithms.
The two crucial tricks are to (i) decompose the horizontal velocity into the contributions of
a scalar streamfunction ψ and scalar velocity potential χ via

V ≡k— (18.58)
ψ+ χ

and (ii) replace the horizontal momentum equations by a vorticity equation and a diver-
gence equation, which are obtained by taking the curl and divergence of (18.57a), respec-

=’ · (ζ + f ) V [vorticity eq.] (18.59)
Dt = k · — (ζ + f ) V ’ V·V [divergence eq.] (18.60)

This vorticity/divergence form is useful in part because of the exact relations
2 2
ζ= ψ & D= χ

Since the spherical harmonics are the eigenfunctions of the Laplacian, (18.61) is equivalent
to the statement that the spherical harmonic coef¬cients of ζ and D are directly propor-
tional to those of ψ and χ, respectively. Thus, the vorticity and divergence equations plus
the “height” equation (18.59) form a closed set of three equations in three unknowns: ψ, χ,
and φ . Furthermore, all the unknowns are scalars and may be expanded directly in spher-
ical harmonics.
The second virtue of the vorticity-streamfunction procedure is that a semi-implicit algo-
rithm for the shallow water set requires implicit treatment of only (i) the 2 φ term in the
divergence equation and (ii) the 2 χ(= D) term in the height equation. With a ¬nite dif-
ference model, one would have to solve two Poisson equations at every time step because
φ and χ appear only in the form of the Laplacian operator acting upon these unknowns.
With a spectral model using spherical harmonics (but only spherical harmonics!), we do
not have to invert the Laplacian because the spherical harmonics are its eigenfunctions:
the Laplacian applied to a spherical harmonic is equivalent to multiplying that harmonic
by ’n(n + 1) where n is the degree of the harmonic.
This in turn leads to the third trick: applying Galerkin™s method instead of the pseu-
dospectral algorithm. Normally, the latter is simpler and easier to program, but a pseu-
dospectral application of spherical harmonics would require inverting large matrices at
every time step.
The fourth trick is to introduce the modi¬ed velocity components U (≡ u» sin[θ]) and
V (≡ uθ sin[θ]); although not unknowns, they are useful auxiliary quantities for evaluating
the nonlinear terms. If we expand

ψ= ψmn Yn
m=’N n=|m|

and similarly for the other ¬elds, identities and recurrence relations for the spherical har-
monics and (18.58) show that the coef¬cients of U and V are

(n ’ 1) Dmn ψm,n’1 ’ (n + 2) Dm,n+1 ψm,n+1 + im χmn (18.63)
Umn =
= ’(n ’ 1) Dmn χm,n’1 + (n + 2) Dm,n+1 χm,n+1 + im ψmn (18.64)


n2 ’ m2
Dmn ≡ (18.65)
4 n2 ’ 1

The ¬fth trick is the pseudospectral strategy of evaluating the nonlinear terms in grid-
point space rather than by forming convolution sums or inner products of the spherical
harmonic coef¬cients. First, the series for 2 ψ, φ , U , and V are evaluated on a latitude
and longitude grid. [Note that the coef¬cients for the 2 ψ are those for ψ multiplied by
’n(n+1).] Then the ¬ve needed nonlinear products are formed: U 2 ψ, V 2 ψ, U φ , V φ ,
and [(U 2 +V 2 )/2] by multiplying the gridpoint values of the factors. Finally, these products
are reexpanded in spherical harmonics.
In a fully pseudospectral treatment, we omit the step of expanding the nonlinear prod-
ucts in basis functions and instead advance the calculation in time on the grid point values.
However, here it is better to advance the spectral coef¬cients in time. The reason is that the
linear terms in the shallow water wave equations ” including the critical terms that we
wish to treat implicitly ” involve the Laplacian operating on an unknown. If we use the
spherical harmonics rather than the equivalent cardinal functions as the unknowns, then
the Laplacians are reduced to factors of ’n(n + 1) multiplying the coef¬cients.

The end result is to convert the shallow water wave equations to the set of ordinary
differential equations in time:
dψmn 1 dPn
’n (n + 1) imAm , Pn ’ Bm ,
dt 2 dµ
+ [n (n ’ 1) Dmn χm,n’1 + (n + 1)(n + 2) Dm,n+1 χm,n+1 ’ Vmn ]

dχmn 1 dPn
’n (n + 1) m
= imBm , Pn + Am ,
dt 2 dµ
’ [n (n ’ 1) Dmn ψm,n’1 + (n + 1)(n + 2) Dm,n+1 ψm,n+1 + Umn ]
+ n (n + 1) Em , Pn + φmn

dφmn 1 dPn
= ’ imCm , Pn ’ Dm ,
+ ¦ n (n + 1) χmn
dt 2 dµ

where µ = cos(θ), and Am (µ), Bm (µ), Cm (µ), Dm (µ), and Em (µ) are the coef¬cients
of the longitudinal Fourier series of the nonlinear terms: U 2 ψ, V 2 ψ, U φ , V φ , and
(U 2 + V 2 )/2, respectively, and the inner product is
a(µ) b(µ)
a(µ), b(µ) ≡ (18.67)

1 ’ µ2

The inner products are evaluated by Gaussian quadrature using the model grid points,
which are the roots of the Legendre polynomials, Pn (µ).
There are many variants to the semi-implicit time-stepping, but the important point is
that all terms are treated explicitly in (18.66) except for φmn (t) in the divergence equation,
(18.66b), and χmn (t) in the height equation, (18.66c). More details are given in Bourke
(1972, 1974) and Bourke et al. (1977). A well-documented and intensively tested shallow
water code is available from the National Center for Atmospheric Research. The same
algorithm is used in the three-dimensional climate model, the CCM3, which is publicly
available at http://neit.cgd.ucar.edu/cms/ccm3/.
The shallow water wave equations are inviscid, but it is trivial to add viscous damping.
The viscosity operator is the Laplacian and the spherical harmonics are its eigenfunctions.

18.20 Fronts and Topography: Smoothing/Filters
18.20.1 Fronts and Topography
Fluid ¬‚ows spontaneously develop narrow, curved regions of very large gradients. These
regions are known as “fronts” in geophysics and “shocks” in aerospace and plasma physics.
The frontal zones narrow to jump discontinuities in the limit of zero viscosity. For phys-
ically reasonable viscosities, which are very small, the frontal zones are much narrower
than the separation between grid points. Therefore, it is necessary to apply some kind of
¬ltering to prevent computational catastrophe.
Another curse is peculiar to meteorology. The topography of the earth™s surface varies


. 85
( 136 .)