√

m ’ ∞, n ¬xed, φ ∼ O

m

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

18.17. ASYMPTOTIC APPROXIMATIONS, II 413

9

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

9

[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

m

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.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

414

9

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.

18.18. SOFTWARE: SPHERICAL HARMONICS 415

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.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

416

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

1

= ’(ζ + f ) k — V ’ V·V (18.57a)

Vt φ+

2

=’ · (φ V ) ’ ¦ D (18.57b)

φt

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-

tively:

=’ · (ζ + f ) V [vorticity eq.] (18.59)

ζt

1

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

φ+

2

This vorticity/divergence form is useful in part because of the exact relations

2 2

(18.61)

ζ= ψ & D= χ

18.19. SEMI-IMPLICIT: SHALLOW WATER 417

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

N N

m

(18.62)

ψ= ψ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)

Vmn

where

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.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

418

The end result is to convert the shallow water wave equations to the set of ordinary

differential equations in time:

m

dψmn 1 dPn

’n (n + 1) imAm , Pn ’ Bm ,

m

(18.66a)

=

dt 2 dµ

+ [n (n ’ 1) Dmn χm,n’1 + (n + 1)(n + 2) Dm,n+1 χm,n+1 ’ Vmn ]

m

dχmn 1 dPn

’n (n + 1) m

(18.66b)

= imBm , Pn + Am ,

dt 2 dµ

’ [n (n ’ 1) Dmn ψm,n’1 + (n + 1)(n + 2) Dm,n+1 ψm,n+1 + Umn ]

1

+ n (n + 1) Em , Pn + φmn

2

m

dφmn 1 dPn

= ’ imCm , Pn ’ Dm ,

m

(18.66c)

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

1

a(µ) b(µ)

a(µ), b(µ) ≡ (18.67)

dµ

1 ’ µ2

’1

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