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

18.23. RADIAL/VERTICAL COORDINATE: SPECTRAL OR NON-SPECTRAL? 429

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

al.(1999).

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

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

430

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

(1984)

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

CRAY-1.

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.

18.25. NON-TENSOR GRIDS: ICOSAHEDRAL, ETC. 431

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-

laria

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

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

432

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

18.26. ROBERT BASIS FOR THE SPHERE 433

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