. 89
( 136 .)


degree of non-uniformity is so small as to be almost invisible when the number of grid
points is large. The rub is that the grid nevertheless is non-uniform and this leaves two
problems: (i) programming is much more dif¬cult than for a uniform grid and (ii) the ac-
curacy of ¬nite differences deteriorates because the grid is not entirely regular. In addition,
recent models have been bedeviled by stability problems and spurious excitation of lon-
gitudinal wavenumber ¬ve disturbances. Heikes and Randall(1995a,b) used a “twisted”
icosahedral grid which respects the parity symmetry with respect to the equator, but at
the expense of more complication. Algorithms that replace the icosahedron with another
Platonic solid, such as a cube or octagon, have been tried, too (Table 18.8).
All of these variants work, but all are messy and most have required some tweaking.
These faults of almost-uniform grid methods are some consolation for the intricacies of
spherical harmonics.

18.26 Spectral Alternatives to Spherical Harmonics, I: Robert
The Fast Fourier Transform (FFT) converts grid point values of a function to its Fourier
coef¬cients in only O(N log2 N ) operations. Since the spherical harmonics are an ordinary
Fourier series in longitude, we can apply the FFT in », but not to the Gegenbauer polyno-
mials in θ. This inspired serious efforts in the late 60™s and early 70™s to ¬nd substitutes for
the latitudinal basis functions.

Andr´ Robert suggested the basis functions

sin|m| (θ) cos(nθ) [“modi¬ed Robert functions”] (18.89)

These basis functions reproduce the polar behavior of the spherical harmonics. The high
order zeros for large wavenumber m eliminate the “pole problem”; one can integrate using
(18.89) with a fairly large time step and no special tricks. Merilees (1973a) pointed out
the fatal ¬‚aw of the “Robert” functions: they are horribly ill-conditioned. The problem is
that for large m, the sin|m| (θ) factor is zero except in a narrow boundary layer around the
equator. Within this layer

n2 θ2
cos(nθ) ≈ 1 ’ ≈1 if |θ| (18.90)
so that the lowest few Robert functions are graphically indistinguishable for large m. This
translates into large round-off errors when a function is expanded in terms of these (al-
most) linearly dependent functions. The result is that the Robert functions have never
been employed for large-scale numerical models.

18.27 Spectral Alternatives to Spherical Harmonics, II: Parity-
Modi¬ed Fourier Series
Merilees (1973b) and Orszag (1974) have discussed parity-modi¬ed Fourier series in which
the latitudinal dependence of the zonal wavenumber m component of a scalar function is
expanded in the form

fm (θ) = sin (θ) an cos(nθ)


m even
s≡ (18.92)
n odd

The pseudospectral grid points are

π (2i ’ 1)
θi ≡ (18.93)
i = 1, . . . , N,
that is, an evenly spaced grid with all points on the interior of θ ∈ [0, π]. This implies
that we may freely divide out sins (θ) and then apply the FFT (a Fast Cosine Transform,
actually). The coordinate singularities of the differential operators at θ = 0 and π cause no
major problems because we do not need to evaluate any quantities at the poles.
One obvious concern is that the spherical harmonics (and all well-behaved functions on
the sphere) have m-th order zeros at the poles whereas the modi¬ed Fourier basis functions
in (18.91) have at most a ¬rst order root. However, to quote Orszag (1974): “The detailed
behavior expressed by the sin|m| (θ) factor is usually of little direct interest in a numerical
simulation . . . The above argument suggests that surface harmonic-series contain much
information on the behavior of functions near the poles that is not of primary interest.”
Boyd (1978b) offers an amusing illustration of Orszag™s point by applying the pseudospec-
tral method with 60 basis functions to calculate the eigenvalues of the two-dimensional

Table 18.9: Spectral Alternatives to Spherical Harmonics: Bibliography
Note: “Double Fourier” denotes a Fourier basis in both latitude and longitude.
References Comments
Robert (1966) “Robert” function basis, but these are very ill-conditioned
Ill-conditioning of “Robert” basis sinm (θ) cos(nθ)
Merilees (1973a)
Merilees (1974) Double Fourier series with high latitude ¬ltering
Orszag (1974) Comparison of spherical harmonics with double Fourier series
Boyd (1978b) Comparison of Chebyshev, Fourier & Associated Legendre
for boundary value and eigenvalue problems
Yee(1981) Poisson equation on sphere; double Fourier basis
Spotz&Taylor high-latitude-¬ltered double Fourier series
& Swarztrauber(1998)
Shen(1999) double Fourier basis
Cheong(2000, 2001) double Fourier series

Laplacian, i. e. the spherical harmonics themselves, for m = 49. Although the 60 term ap-
proximation may have at most 61 zeros (counting multiplicity and the two roots of sin(θ))
and the spherical harmonics have 98 + (n ’ m) roots, the lowest twenty eigenvalues are
calculated to within an accuracy of 1 part in 70 million.
For so-called “jury” problems, that is to say, for boundary value and eigenvalue prob-
lems, Boyd (1978b) argues that the parity-modi¬ed Fourier series is the best basis set. It is
much easier to program than spherical harmonics and the accuracy does not seem poorer.
For time-dependent or “marching” problems, however, the “pole” problem returns in
full force. Orszag (1974) has shown that the time-step limits for realistic calculations will
be quite severe. He suggested several remedies. First, one may apply a polar ¬ltering ev-
ery ten or twenty time steps. Second, one may split the convective terms into (i) a surface
harmonic of degree 2 that gives the ¬‚ow over the poles and (ii) the rest of the convective
terms and then treat (i) implicitly; the bulk of the terms, i. e. (ii), may be treated explic-
itly without instability because (ii) does not involve ¬‚ow over the pole. Third, treating
viscous terms implicitly ” almost always necessary with spectral and pseudospectral al-
gorithms ” tends to stabilize the convective terms, too, at least for a Backwards-Euler or
other implicitly-dissipative scheme.
Even so, longitude/latitude double Fourier series have been little used for time-dependent
calculations in the quarter-century after Orszag™s article. Part of the reason is inertia; as ex-
plained earlier, the relatively slow cost of associated Legendre transforms has not been
a major liability for spectral models even at the highest resolution employed in weather
forecasting at the dawn of the third millenium, T 639.
In recent years, however, double Fourier series have made a comeback in the articles
of Shen(1999), Cheong(2000, 2001) and Spotz, Taylor and Swarztrauber(1998). Spotz et al.,
for example, revived an idea of Merilees(1974) with an improved ¬lter, and showed that
latitudinal Fourier series can be both fast and accurate when the polar ¬ltering is strong

18.28 Grid-to-Spherical-Harmonics-to-Grid Projective Filter
for Time-Dependent Latitudinal Fourier Series
The ultimate polar ¬lter is to transform grid point values on a tensor product latitude/longitude
grid to a triangular truncation of spherical harmonics and then sum the spectral series to
obtain a different set of grid point values. The output of the fore-and-back transform will

inherit the equiareal resolution property of spherical harmonics. The pole problem is com-
pletely suppressed.
This combination of transforms has been described as both a “spherical ¬lter” and a
“projection” in the literature; we shall call it a “projective ¬lter”. The reason for the termi-
nological confusion is that the projective ¬lter is different from other ¬lters in that its pur-
pose is not to provide arti¬cial viscosity or suppress small-scale noise and aliasing instabil-
ity. Rather, the grid-to-spectral half of the projective ¬lter is a non-dissipative transform. It
is, in the usual parlance of functional analysis, a “projection” of the function f (», θ) onto the
subspace spanned by the spherical harmonics. The spectral-to-grid half of the transform
is merely a summation of the truncated series. The grid-to-grid transform is more than a
projection (because it is a summation, too), but in some sense less than a ¬lter (because the
only ¬ltering is a truncation).
The weakness of a latitude/longitude Fourier method is that projective ¬ltering has a
cost which is the same order of magnitude as an Associated Legendre transform. However,
because only the unknowns themselves must be projectively ¬ltered ” no derivatives need
be transformed ” double Fourier algorithms have the potential of extending the lifetime
of global spectral methods on the sphere for at least another decade.
The method which Spotz and Swarztrauber (2000) label both “naive” and “standard”
is to compute both halves of the projective ¬lter through separate Parity Matrix Multipli-
cation Transforms (PMMT) in latitude. (All transformations in longitude are performed
by the Fast Fourier Transform (FFT), which convert grid point values into functions fm (θ)
which describe how the coef¬cient of cos(m») or sin(m») varies with colatitude.) The ad-
jective “Parity” means that grid point values of f (θ) are ¬rst split into components that
are symmetric and antisymmetric with respect to the equator θ = π/2 through f S ≡
{ f (θ) + f (π ’ θ) } and f A ≡ { f (θ) ’ f (π ’ θ) }. Each component is then projectively
¬ltered separately; the even and odd parity ¬ltrates are recombined by adding their grid
point values in the northern hemisphere and subtracting the ¬ltered values of f A from the
¬ltrate of f S in the southern hemisphere.
Yarvin and Rokhlin(1998) noted that one could reduce the cost by 25 % by combining
the forward and backward transform matrices into a single N — N matrix for small m
while performing the forward and backward transforms separately for large m. (Note
that the forward and backward transform matrices are rectangular with smaller dimension
(N ’ m); for the extreme of m = N , the transform matrices are only vectors; two vector-
vector multiplies are obviously cheaper than an N — N matrix-vector multiply.)
Swarztrauber and Spotz (2000) developed the “Weighted Orthogonal Complement”
projective ¬lter which has only half the cost of the “naive”/“standard” method and only
two-thirds the expense of Yarvin and Rokhlin™s “direct” method. Better yet, they prove
that the WOC matrices can be written as subsets of a couple of “master” matrices, reduc-
ing storage for transform/projection matrices from O(N 3 ) to only O(N 2 ).
For the shallow water wave equations, nine Legendre transforms are needed at each
time step with a spherical harmonic basis. With a double Fourier expansion and the WOC
projective ¬lter, the transform cost is reduced to the equivalent of only three Legendre
Jakob-Chien and Alpert(1997) developed a projective ¬lter based on Fast Multipole
Method (FMM) ideas; this was improved by Yarvin and Rokhlin(1998). As with other
multipole transforms, the FFM projective ¬lter is faster than the WOC or any MMT-type
method in the limit N ’ ∞, but the cross-over point where the FMM scheme is faster is
likely to be at very large N .
Spotz and Swarztrauber(2000) carefully compare all these projective ¬lters for various
N , both with and without the application of the Two-Thirds Rule. (Dealiasing truncation
reduces the advantages of the newer MMT schemes relative to the “naive” method.) They

Table 18.10: Projective Filters Bibliography

References Comments
Jakob-Chien&Alpert(1997) Grid-to-spherical-harmonics-to-grid
Yarvin&Rokhlin(1998) projective ¬lter for longitude/latitude
time-dependent double Fourier series
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
Cheong(2001) Approximate projective ¬lter (hyperviscosity)
which combines damping with projection

¬nd that their WOC algorithm is faster than the “direct” and “standard” PMMT and the
FMM at least up to N = 200. The WOC method is moderately faster than formal operation
counts would suggest; the very low storage of the WOC leads to more ef¬cient utilization
of the on-chip cache memory of current machines. The WOC method also parallelizes very
Cheong(2001) employs an alternative: Polar ¬ltering by a combination of fourth order
and sixth order hyperviscosity. Recall that Yn is an eigenfunction of the two-dimensional
Laplace operator on the surface of the sphere with the eigenvalue »n = ’n(n+1). It follows
that a dissipation which is proportional to a high power of the Laplacian can be tuned to
strongly damp all spherical harmonics with n > N while having only a slight effect on
harmonics retained within a triangular truncation. The Fourier-Galerkin discretization of
the Laplace operator is a tridiagonal matrix, so the ¬lter is very fast.
In contrast to the Spotz-Swarztrauber and FMM methods, Cheong™s hyperviscosity is
not an exact projection onto spherical harmonics. Since most time-dependent models re-
quire a little scale-dependent dissipation for noise-suppression anyway, an approximate
projective ¬lter is satisfactory. However, there has not been a detailed comparison of the
relative speed and accuracy of exact projective ¬lters versus approximate ¬lters that damp
and project simultaneously.
Swarztrauber and Spotz (2000) note, “The goal of an O(N 2 log N ) harmonic spectral
method remains elusive; however, the “projection” method provides a new avenue of re-
search. Perhaps the development of a fast projection [projective ¬lter] will prove to be
easier than the development of a fast harmonic transform.” Suf¬cient progress has al-
ready been made to encourage the development of double Fourier methods on the sphere.
However, as of 2000, latitude/longitude Fourier series have not yet been applied to an
operational forecasting or climate model.

18.29 Spectral Elements on the Sphere

Another alternative to global spectral methods is domain decomposition. A moderately
high order Legendre series is used in each subdomain, and then the pieces are stitched
together as in ¬nite elements to obtain a global approximation. The great virtue of spec-
tral elements is that each subdomain can be assigned to a different processor of a multi-
processor machine. The articles listed in Table 18.11 show that the algorithm parallelizes
very well, and seems promising for both atmospheric and oceanic modelling.

Table 18.11: Spectral Elements on the Sphere Bibliography

References Comments
Ma(1992,1993a,b) Ocean basin models
Iskandarani&Haidvogel global ocean model
Levin&Iskandarani ocean model
Taylor&Tribbia & atmospheric model;
Iskandarani(1997) shallow water equations
Curchitser&Iskandarani Ocean modelling on a variety of
& Haidvogel(1998), massively parallel computers
&Iskandarani &Hughes


. 89
( 136 .)