. 88
( 136 .)


normally be vectors themselves. If the modes differ drastically ” for example, the grav-
ity wave oscillations of the shallow water wave equations have very high frequencies in
comparison to the Rossby modes ” than it may be possible to greatly improve numerical
ef¬ciency by choosing vector basis functions that mimic the important modes. When we
expand all the unknowns in the shallow water or primitive equations in spherical harmon-
ics, it is not possible to discriminate in any way between the gravity and Rossby modes.
One can discriminate using vector basis functions.
Several alternatives to ordinary (scalar) spherical harmonics have been discussed. Swarz-
trauber (1981, 1993) and Swarztrauber and Kasahara (1985) have described two different
varieties of “vector spherical harmonics”. The basic strategy is to transform the spherical
components to the corresponding Cartesian components (which are always continuous)
and then expand in (scalar) spherical harmonics. It is convenient to label Swarztrauber™s
new basis set the “vector spherical harmonics” even though he takes pain to note that
this term has been previously applied to a different set: the eigenfunctions of the vector
Helmholtz equation. Swarztrauber has written library software (“SPHEREPACK”) to com-
pute his new basis functions and has applied them to Laplace™s tidal equations, which are
the linearized shallow water wave equations.
Shen and Wang(1999) have recently proposed a new vector spherical harmonic/Legendre
algorithm for the primitive equations. The theoretical analysis and numerical experiments
are promising, but their algorithm has not yet been used for meteorological applications.
The second approach, suggested by T. W. Flattery for data analysis but developed for
time-integration by A. Kasahara, is to use the vector eigenfunctions of Laplace™s tidal equa-
tions as the basis set. In three dimensions, separation-of-variables shows that after the ver-
tical dependence has been eliminated, the two-dimensional equations are simply Laplace™s
tidal equations with a different equivalent depth parameter for each vertical mode. Thus,
although the basis set is two-dimensional, these functions of latitude and longitude be-
come the three-dimensional normal modes when multiplied by the proper vertical struc-
ture function.
The motive for this choice of basis set is that these so-called Hough functions are inti-
mately connected with the physics. In particular, to ¬lter gravity waves, which are “noise”
insofar as the forecast is concerned, one may simply omit the gravity wave Hough modes
from the basis set.
There are alternatives for eliminating the gravity waves; the most famous is the quasi-
geostrophic approximation, which drops terms from the shallow water wave equations
such that the simpli¬ed partial differential equations have only Rossby waves as their free
modes of oscillation. Unfortunately, the errors in quasi-geostrophy (and its generaliza-
tions) are large at low latitudes. Furthermore, deleting terms from the differential equation
is an all-or-nothing approach: the only choice is between models that allow all gravity
modes or none at all.
In contrast, Hough models allow selectivity. There are a few gravity waves with rel-
atively low frequencies which can be retained in a trucation of the Hough basis, but are
always excluded by the quasi-geostrophic approximation.
One obvious dif¬culty is that the Hough functions are the eigenmodes of the linearized
equations, but the dynamics of the global atmosphere is strongly nonlinear. One may well

ask: Doesn™t the nonlinearity produce such strong mode-mode coupling, such strong dis-
tortions of the normal modes themselves, as to invalidate the physical interpretation of the
Hough modes as the gravitational and Rossby oscillations of the atmosphere?
The answer to this question is: Yes, but not so as to diminish the usefulness of Hough
functions as a basis set. The reason is that because of the great differences between Rossby
and gravity waves, the couplings produced by nonlinearity are primarily between modes
in the same class ” Rossby with Rossby, gravity with gravity. What is important for se-
lectively ¬ltering gravity waves is not whether a particular Rossby Hough function closely
resembles a free oscillation as observed in the real atmosphere. Rather, what matters is that
even with nonlinearity, dissipation, and other complications, the Rossby Hough functions
are, to a good approximation, a complete set for expanding the low-frequency, weakly di-
vergent motions that dominate weather and climate and this set is orthogonal to the gravity
Hough functions. The Rossby-to-gravity coupling is weak.
Kasahara (1977, 1978) has clearly demonstrated the feasibility of integrating the hydro-
dynamic equations using Hough functions. It is important to note that one can resolve
a slow ¬‚ow with roughly one-third the number of degrees of freedom of a conventional
spherical harmonics model since there are two gravity modes for every Rossby Hough
function. However, like Swarztrauber™s equally feasible vector basis set, Kasahara™s Hough
functions have not yet caught on, except for the middle atmosphere model of Callaghan et
There is a technical reason for this: when the number of vertical levels is large, the eigen-
values for the higher vertical modes (what meteorologists call the “equivalent depths” of
the modes) become very small. This disrupts the clear separation between Rossby and
gravity waves which exists for the lowest few vertical modes and is the motive for choos-
ing a Hough basis in the ¬rst place.

18.23 Radial/Vertical Coordinate: Spectral or Non-Spectral?
18.23.1 Basis for Axial Coordinate in Cylindrical Coordinates
For cylindrical coordinates, the z-axis is not cyclic, so Chebyshev polynomials (if the cylin-
der is bounded in z) or rational Chebyshev functions (if the cylinder is unbounded) are
best. Thus, a cylindrical basis is usually double Chebyshev (in r and z) and Fourier in the
angle θ.

18.23.2 Axial Basis in Toroidal Coordinates
Toroidal coordinates are similar to cylindrical except that the vertical axis is closed upon
itself to form a circle instead of an in¬nite line. Because this closure-into-a-circle forces the
“toroidal angle” » to be cyclic, one should always use an ordinary Fourier series in this
coordinate. Thus, a toroidal basis is Fourier in θ, Fourier in », and Chebyshev in the third
coordinate, which is radial distance from the center of the torus.

18.23.3 Vertical/Radial Basis in Spherical Coordinates
Although spherical harmonics have been extremely popular in atmospheric science, such
models almost invariably use ¬nite difference approximations in the vertical. It is not that
the idea of a vertical series has been tried and found wanting, but rather it has not been
tried. Francis(1972) and Hoskins(1973) explored Laguerre functions (for an atmosphere of

semi-in¬nite depth), but found that explicit time step limits would be very severe. Machen-
hauer and Daley(1972, 1974) and Eliassen and Machenauer (1974) experimented with Leg-
endre polynomials in depth, but had to use ad hoc modi¬cations; this work was never
published in a journal.
Kasahara and Puri(1981) and Mizzi, Tribbia and Curry(1995) tried the eigenfuctions
of the so-called “vertical structure equation” as the basis functions in height. Mizzi et al.
found that the rate of convergence was lamentably slow compared to the expansion of
test solutions in Chebyshev or Legendre polynomials. This is hardly surprising: using
normal modes for a problem with numerical boundary conditions is akin to using a Fourier
basis for a non-periodic problem and gives the same slow, algebraic rate of convergence.
(Indeed, for idealized temperature pro¬les, the normal modes are sine functions.)
And there the matter has rested except for Shen and Wang(1999). In contrast, Mar-
cus and Tuckerman (1987a, b), Mamun and Tuckerman (1995) and Glatzmaier(1984, 1988)
have computed the ¬‚ows in a rotating spherical annulus using Chebyshev polynomials in
radius. Glatzmaier™s algorithm was so successful that it has been used for studies of stellar
dynamos, mantle convection, and the atmosphere of a gas giant planet totalling at least
thirty-¬ve articles and book chapters. Because Glatzmaier™s algorithm has been used so
widely, it is worth brie¬‚y describing here.

18.24 Stellar Convection in a Spherical Annulus: Glatzmaier
This model simulates stellar convection by solving the anelastic equations of motion for a
spherical shell heated from below. The inner wall is at 0.53 of the solar radius and the outer
wall at 0.93 of the solar radius so that there are seven scale heights within the shell.
The latitude-longitude basis was spherical harmonics in a triangular T31 truncation.
The radial basis was Chebyshev polynomials up to and including T16 (r).
The time-integration scheme was semi-implicit. The explicit Adams-Bashforth scheme
was applied to all terms except diffusion, which was treated implicitly with the Crank-
Nicholson scheme. Without magnetic forces, the cost was about 1.0 sec per timestep on a
The implicit Crank-Nicholson scheme does not require splitting because the spheri-
cal harmonics are eigenfunctions of the horizontal part of the diffusion operator. Because of
this, it is only necessary to solve a one-dimensional boundary value problem in radius only
to implement the Crank-Nicholson method. Ordinarily, this would require storing the
LU decomposition of O(L2 ) Chebyshev collocation matrices, one for each horizontal basis
function. However, the eigenvalue of the diffusion equation is - n(n + 1), independent of
the zonal wavenumber m. Consequently, Glatzmaier is only obliged to perform (N +1) LU
decompositions in a preprocessing stage and then store these thirty-two 17 x 17 matrices
throughout the run. The cost of the backsolves for the Crank-Nicholson treatment of diffu-
sion is O(L2 N 2 ) whereas the number of gridpoints is only O(L2 N ), so we have lost a factor
of N relative to a ¬nite difference method. This will be only a small part of the overall work
“ this is a very complicated set of equations because of the magnetic and Coriolis forces “
unless N is a good deal larger than 17.
The smallness of N is a tribute to the accuracy of pseudospectral methods. The density
varies by a factor of exp(’7) over the shell,4 and there are viscous boundary layers at both
the top and the bottom, but seventeen grid points in radius suf¬ce. The model does not use
a variable grid because the built-in quadratic stretching of the Chebyshev grid is suf¬cient.
4 This is the physical meaning of the statement that the shell is “seven scale heights” thick.

Glatzmaier reports that a test calculation with N = 65 differed from the results for N = 17
by only 2 % after 300 time steps.
Although meteorological spherical harmonics codes are almost always dealiased, the
stellar convection model is aliased; no 2/3™s rule here! Glatzmaier states the position of that
large group of modellers who ignore aliasing: “An unaliased method, which would require
signi¬cantly more computer time and memory, would ensure that the error due to the
truncation of the harmonics is normal [i. e., orthogonal ] to the retained set of harmonics.
However, a truncation error would still exist... If the two methods [aliased and unaliased]
produce solutions that are not negligibly different, both solutions are bad. Only when
enough harmonics are retained so that both methods produce negligibly different solutions
are the solutions good”.
The incompressibility constraint is enforced by decomposing the three dimensional ve-
locity into

— — (W r) + — (Z r) (18.88)
v= ˆ ˆ

where r is a unit vector in the radial direction and ( —) denotes the curl. Through elemen-
tary vector identities, one may verify that the divergence of (18.88) is zero. The functions
W and Z are scalar, and generate the “poloidal” and “toroidal” components of the ¬‚ow.
The incompressibility has reduced the number of scalar functions needed to specify the
¬‚ow from three to two, and similarly for the magnetic ¬eld, which is also nondivergent
(unless the Sun has magnetic monopoles!) This poloidal/toroidal decomposition is the
three-dimensional, spherical analogue of the usual vorticity streamfunction formulation,
but the pressure is not eliminated and remains one of the six unknowns: W , Z, the two
similar magnetic functions, p, and the entropy s.
The continued presence of the pressure as an unknown creates a problem. At the rigid
boundaries, the constraints of no normal ¬‚ow and no stress imply four conditions on W ,
which satis¬es only a second order equation in r. However, the pressure p also satis¬es a
second order equation, and there are no intrinsic boundary conditions on p. Glatzmaier
therefore treats part of the radial diffusion for W explicitly by including the terms that
depend on p. He then demands that p have whatever values at the boundary are necessary
to make the explicit terms for W equal to zero so that W remains consistent with the four
boundary conditions upon it.
This sounds awkward when expressed in words; it is not a lot of laughs for the pro-
grammer either. The pressure is a perpetual nuisance in models which lack a prognostic
equation for p. Pressure gradients appear in the governing equation, just as gradients of
the velocity do, but we have physical boundary conditions only on v. The pressure is a
passive quantity which varies so that the incompressibility constraint is always satis¬ed,
and we have to live with this awkwardness.

18.25 Non-Tensor Grids: Icosahedral Grids and the Radio-
Spherical harmonics, in which every part of the globe directly in¬‚uences every other part,
have always alarmed proponents of parallel computing because of the need to share lots
of grid point values between different processors of a massively parallel machine. This has
not become a dif¬culty as early as expected because even at T639 resolution, preliminary
tests have shown that Legendre transforms cost only 19% of the overall CPU time using 48
vector processors. (The model at T213 resolution runs quite ef¬ciently on 500 processors of
a Cray T3E.) Nevertheless, there has been a search for ¬nite difference and ¬nite element

alternatives to spherical harmonics. Unfortunately, the globe imposes special dif¬culties
even on low-order methods: there is no free lunch on the sphere.
One obvious alternative to pseudospectral methods is to simply lay down an evenly
spaced grid or triangulation of the sphere and use ¬nite differences or ¬nite elements.
Unfortunately, an evenly spaced grid and a triangulation composed of identical triangles
are both impossible except at very low resolution.
The choice of an evenly spaced latitude-longitude grid creates a severe pole problem
for all non-spherical harmonic schemes, both low order differences and high order spectral
methods, unless (i) the time-stepping is fully implicit, which is prohibitively expensive or
(ii) strong ¬ltering is applied near the poles.
The other grid option is to abandon a tensor grid and instead simply spread grid points
as uniformly as possible over the sphere. As reviewed by Saff & Kuijlaars(1997), this turns
out to be very dif¬cult.
The problem was identi¬ed by Leonard Euler in the middle of the eighteenth century,
but it ¬rst became important in science when marine biologists on H. M. S. Challenger
dredged up samples of microscopic ocean life a century ago. Fig. 18.18 shows the “radio-
laria”, Aulonia hexagona. At ¬rst glance, the siliceous skeleton appears composed of an even
net of hexagons. One might imagine that one could circumscribe the skeleton with a globe
and subdivide each hexagon into six equilateral triangles to obtain a uniform triangulation
of the sphere, or compute ¬nite differences directly on the hexagonal grid. Unfortunately,
the appearance of a uniform net of hexagons is an illusion. W. D™Arcy Thompson (1917)
shows that every such hexagonal array must have 12 pentagons inserted in the lattice.
The vertices of the ¬ve Platonic solids, in contrast, are uniformly spaced. However,
the icosahedron, with twelve vertices and twenty equilateral triangles, is the Platonic solid
with the largest number of sides. One can hardly hope to model the weather well with just
twelve grid points!
The best one can hope for is to create an almost uniform triangulation. Baumgardner and
Frederickson (1985) show how one can use the icosahedron as a starting point for covering
the globe with little triangles. Each of the twenty sections is bounded by a triangle of
identical size, but the triangles within each icosahedral face are of different sizes.

Table 18.8: Gridpoint Methods on the Sphere Bibliography

References Comments
Taylor (1995) proves spherical harmonics interpolation problem
is insoluble;
expansion must be computed
by integration with extra points
Sloan(1995) interpolation and “hyperinterpolation” on sphere
Baumgardner&Frederickson(1985) Icosahedral/triangular grid
Heikes&Randall(1995a,b) Twisted icosahedral grid
Ronchi&Iacono New ¬nite differences (“cubed sphere”)
& Paolucci(1997) compared with spherical harmonic code
Ranˇ i´ &Purser
cc Expanded cube grid: comparisons of
Mesinger(1996) gnomonic and conformal coordinates
Purser&Rancic(1997) Finite difference grid based on conformal octagon
Swarztrauber&Williamson icosahedral grid; 3D Cartesian method
&Drake(1997) projected onto the surface of a sphere
Thuburn(1997) icosahedral-hexahedral ¬nite difference grid
Saff & Kuijlaars(1997) Popular review of the mathematical dif¬culties
in computing nearly-uniform grids on the sphere
Stuhne&Peltier(1999) icosahedral grid

Figure 18.18: The skeleton of aulonia hexagona. After Growth and Form by W. D™Arcy Thomp-
son (1917). Although the lattice appears to be composed entirely of hexagons, Thompson
showed that there must be at least a dozen pentagons, too. It is impossible to embed the
sphere in a lattice of identical polygons when the number of polygons is greater than 20.

These nearly-uniform grids work well ” as evident in the picture of the radiolaria, the


. 88
( 136 .)