. 60
( 136 .)


properties, but the “off-grid” interpolation is where the “conservation” and “positive de¬-
niteness” wars must be fought.

14.10.2 Spectral Off-grid
Spectrally-accurate off-grid interpolation algorithms include the following:
1. Fast Multipole Method (FMM) (Boyd, 1992c, Dutt & Rohklin, 1993, 1994)
2. Partial Chebyshev expansion method (Suli & Ware, 1991, 1992, Ware, 1991,1994)
3. Fine-grid Lagrange interpolation algorithm (Boyd, 1992)
As explained at greater length in Chapter 10, all three methods are asymptotically faster
than direct summation. (See the review by Ware, 1998). However, the proportionality
factors, which vary with the desired error tolerance, are always much greater than for the
Fast Fourier Transform.
The Suli-Ware algorithm is the only one that has been explicitly applied in semi-Lagrangian
calculations. It performed well in both one and two space dimensions.
Boyd™s ¬ne-grid polynomial interpolation algorithm is probably both the fastest and
the simplest to explain. The key idea is that low order polynomial interpolation can be
very accurate if performed on a suf¬ciently ¬ne grid. So, his ¬rst step is to apply the FFT
to transform the data to a ¬ner grid than used for the rest of the calculations. One can then
use ordinary Lagrange interpolation or splines on the ¬ne grid with little loss of accuracy
compared to spectral off-grid interpolation on the original coarser grid.
His algorithm has not been explicitly used in semi-Lagrangian models to date. However,
the ECMWF forecasting model and the NCAR Community Climate Model both employ
his idea implicitly. For reasons of cost, and because fast off-grid interpolations were not yet
known at the time, the programmers used various forms of low-order interpolation and
hoped for the best. (Careful testing showed that this inconsistent use of low order off-grid
interpolation in a high order model does not signi¬cantly degrade the overall forecast.)
However, these models employ dealiasing so that inner products are evaluated on a grid
which is much ¬ner than the actual number of spectral degrees of freedom. Thus, these
models implicitly use Boyd™s strategy.
The inconsistency between spectral calculation of derivatives and non-spectral off-grid
interpolation is in any event not as illogical as it seems. The reason is that to minimize
dispersion, it is important to have accurate computation of derivatives in the linear terms
which are not treated by semi-Lagrangian advection. The off-grid interpolation does not
have to approximate derivatives, but only the velocities themselves, which is easier.
There are two open questions about spectral off-grid interpolation. First, is it cost-
effective in a real climate or forecasting model, or is low order polynomial interpolation
suf¬cient? The second is: can spectral semi-Lagrangian schemes be improved, when nar-
row fronts form, by skillful use of ¬lters or sum-acceleration methods?
We discuss ¬lters brie¬‚y in Chapter 18, Sec. 20, but it is useful to at least hint at the
possibilities. In Fig. 14.9, there is a large spike of error in the immediate neighborhood
of the front in both panels; it is very dif¬cult to eliminate this because even a slight er-
ror (O(h)) in the location of the front will inevitably generate O(1) errors, although the
shape of the computed ¬‚ow may look very similar to the true ¬‚ow. Away from this small
neighborhood, however, Fig. 14.9 shows that ¬ltering can dramatically reduce errors in a
Fourier semi-Lagrangian method. There is only a little experience with ¬ltered spectral
semi-Lagrangian schemes, and ¬lters themselves are a hot research topic. The graph sug-
gests that further work on ¬ltered semi-Lagrangian algorithms would be rewarding.

14.10.3 Low-order Polynomial Interpolation
Finite difference and low order ¬nite element models commonly use polynomial interpo-
lation to approximate quantities off-grid. Cubic splines (Purnell, 1976) have been tried, but

. Filtered
1 .1 Unfiltered

.01 .01

.001 .001

π -π π
-π 0
Figure 14.9: Errors in the solution to the one-dimensional advection equation by Fourier
three-level semi-Lagrangian computation with N = 64 grid points and a time step „ = 1/2.
The initial condition is u(x, 0) = 1 ’ sin(x). Left: Vandeven ¬lter of order p = 6. Right:
un¬ltered numerical solution.

the most common choice for interpolating particle properties (i. e., u’ in the SL schemes
above) is cubic Lagrange interpolation. (Leslie and Purser, 1991, used quintic interpola-
tion.) There is less consistency in calculating displacements because as noted in the pre-
vious subsection, the forecast error is much less sensitive to errors in the velocity in the
¬xed-point iteration for ± than to errors in u’ , F 0 , etc. Some therefore use linear interpo-
lation of u(x ’ ±, tn ) in calculating displacements. McDonald (1986) suggested a hybrid
strategy: linear interpolation on the ¬rst ¬xed-point iteration followed by cubic interpola-
tion on the second (and usually last) iteration.

Moisture has a very patchy distribution in the atmosphere; the air above the Mediter-
ranean Sea is very humid because the warm water evaporates like crazy, but the Mediter-
ranean Sea is bordered by the Sahara desert, which is not humid at all. Global spectral
models, with a resolution of only O(100 km), exhibit Gibbs™ phenomenon in the vicinity
of such moisture fronts. The resulting regions of negative water vapor are decidely un-
physical and have a bad effect on model climate. SL schemes do much better. (Rasch and
Williamson,1990b,1991 and Williamson, 1990).

The moisture transport problem has motivated studies which replace ordinary inter-
polation by a procedure in which the interpolant is adjusted by side constraints, such as
smoothness and positive-de¬niteness. Such “quasi-monotone” or “shape-preserving” SL
schemes have been investigated by Williamson and Rasch(1989), Rasch and Williamson(1990a)
and Bermejo and Staniforth (1992).

Adjusted or constrained interpolation can also be used to enforce conservation of total
mass. Priestley (1993) and Gravel and Staniforth (1994) have generalized SL to be quasi-

14.10.4 McGregor™s Taylor Series Scheme: an Alternative for Computing
Departure Points
The only disadvantage of polynomial interpolation is that if done accurately, it is rather
expensive, especially in multiple space dimensions. McGregor (1993) suggested an alter-
native procedure, which is to calculate the trajectory by a Taylor series in time about the
departure point. The catch is that because the trajectory is the goal, the time derivatives
must be total time derivatives. Thus, in one space dimension,

„2 „3 ‚ ‚
± = ’„ u + (ut + uux ) ’ (14.53)
+u (ut + uux ) + . . .
2 6 ‚t ‚x

Calculating the local time derivatives above is very hard. However, McGregor argues
that these can be neglected compared to the advective terms so that

„2 „3 ‚
± ≈ ’„ u + uux ’ u (14.54)
(uux ) + . . .
2 6 ‚x
where the velocities are to be evaluated at the spatial location of the arrival point at the
midpoint in time (tn for the three-level scheme, tn+1/2 for a two-level method).
His method is an order of magnitude cheaper than bicubic interpolation and accurate to
second order. Adding additional terms does not increase the formal time accuracy beyond
second order, but McGregor shows that empirically, including partial third order terms
does improve accuracy considerably for ¬nite „ .

14.11 Higher Order SL Methods
Purser and Leslie (1994) have developed an ef¬cient third-order SL scheme. Earlier, Maday,
Patera and Rønquist (1990) and Suli and Ware (1991, 1992) and Ware (1991, 1994) have
developed a family of algorithms which can be extended to arbitrary order. Building on
Holly and Preismann (1977) and Yang and Hsu (1991), Yang, Belleudy and Temperville
(1991) developed a ¬fth order accurate SL scheme.
The bad news is that these high order schemes are all rather complicated. The assump-
tion that a particle is advected at constant speed by the velocity at the midpoint of the
trajectory is suf¬cient to give second order accuracy in time. At higher order, the displace-
ment on t ∈ [tn’1 , tn ] is no longer the same as t ∈ [tn , tn+1 ], and a simple ¬xed point
iteration no longer gives direct access to the trajectory.
The Maday-Patera-Rønquist and Suli-Ware algorithms require solving auxiliary time-
marching problems, k problems for a (k + 1)-th order method. Maday et al. show that
one can “sub-cycle”, integrating these auxiliary advection problems with a short time step
while applying the implicit time-marching formulas (for non-advective terms) only over a
long time interval.
Even so, these algorithms are uncompetitive with the second order SL schemes for me-
teorological applications. We shall therefore not describe them in detail. These are likely
to become more important as resolutions improve so that time-truncation errors begin to
catch up with space truncation errors. The latter are currently much worse for forecasting
and climate models than the temporal differencing errors.

14.12 History and Relationships to Other Methods

Table 14.1: A Selected Bibliography of Semi-Lagrangian Methods
SL is an abbreviation for “semi-Lagrangian”; SI for “semi-implicit”.

Reference Comments
Fjørtoft(1952, 1955) First Lagrangian algorithm (meteorology)
Wiin-Nielsen (1959) Earliest SL, at least in meteorology
Sawyer (1963) coined “semi-Lagrangian”
Holly&Preissman (1977) SL variant using Hermite interpolation (hydrology)
Haltiner&Williams(1980) Good review of (singly-upstream) SL
Robert (1981, 1982) Earliest pairing of SL with SI
Douglas&Russell(1982) Invention of an SL scheme
Pironneau (1982) Invention of SL as “transport-diffusion” algorithm
Bates&McDonald (1982) First multiply-upstream SL
Bates (1984) More ef¬cient multiply-upsteam SL
Morton (1985) Invention of SL as “Lagrange-Galerkin” scheme
Robert&Yee&Ritchie(1985) SL with SI
Pudykiewicz et al.(1985) Convergence of ¬xed point iteration
Temperton & Staniforth(1987) Independent invention of two-level SL
McDonald & Bates (1987) 2d order accuracy through extrapolated velocities
Ritchie (1986) Invention of “non-interpolating” SL
Ritchie (1987,1988,1991) Non-interpolating SL applied to a sequence of
increasingly sophisticated spectral models
McCalpin (1988) Analysis of dissipation inherent
in SL interpolation
McDonald&Bates(1989) Showed SL with time-splitting was inaccurate
Bates&Semazzi& Finite difference SL shallow water model
Higgins&Barros(1990) Accuracy improved by avoiding time-splitting
Cot´ &Gravel&Staniforth(1990) Finite-element SL with “pseudostaggering”
Williamson&Rasch(1989) SL with shape-preserving (positive de¬nite)
Rasch&Williamson(1990a) advection
Rasch&Williamson(1990b,1991) moisture transport in climate & weather models
Purser&Leslie(1988,1991,1994) Not spectral but improvements in SL
Leslie&Purser(1991) latest article is 3rd order in time
Yang&Belleudy Fifth order SL
Yang&Hsu(1991) Improved Holly-Preismmann & Yanget al. schemes
Smolarkiewicz & Rasch(1990) Generalization of non-interpolating SL; monotone
Maday&Patera&Rønquist(1990) SL of arbitrarily high time-order
Suli&Ware(1991,1992) Fourier and Chebyshev SL/SI, proofs and examples
Ware (1991,1994) independent invention of arbitrary time-order SL
Baker&Suli&Ware(1992) & fast off-grid interpolation scheme
Staniforth& Cot´ (1991) Highly readable review
Gravel&Staniforth&Cot´ (1993) Stability analysis

Table 14.1: Bibliography of Semi-Lagrangian Methods (continued)

Reference Comments
Bermejo&Staniforth(1992) Quasi-monotone SL
Priestley(1993) Quasi-conservative SL
Gravel&Staniforth(1994) Quasi-conservative SL for shallow-water system
Inexpensive O(„ 2 ) formula for departure pts.
Pellerin&Laprise Experiments with simple advection/
&Zawadzki(1995) condensation system
Cot´ &Gravel Extra time level scheme to defeat
&Staniforth(1995) orographic resonance
Ritchie&Tanguay(1996) Comparison of spatially averaged Eulerian
& Lagrangian treatments of mountains
Bottcher (1996) Modi¬ed exponential splines for interpolation
Makar&Karpik(1996) Periodic B-spline interpolation on sphere for SL;
comparisons with spectral methods

The Lagrangian (particle-tracking) formalism is a couple of centuries old. Fjørtoft (1952,
1955) showed that particle-tracking could be used for weather forecasting. However, We-
lander (1955) showed through a vivid and widely- reproduced graph that blobs of ¬‚uid
were rapidly stretched and deformed into a very convoluted shapes. Today, we would
say that ¬‚uids can be chaotically mixed, in the sense of “chaos” theory, and a blob of ¬‚uid
initially bounded by a small square is bounded, for later times, by a “fractal” curve of
ever-increasing length.
Wiin-Nielsen (1959) therefore introduced a “semi-Lagrangian” scheme: just as in mod-
ern schemes, each set of particles is tracked only over a single small time interval of length
„ . A new set of particles is chosen for later time steps: those which will arrive at the
points of an evenly spaced grid. Later meteorological papers by Krishnamurthi (1962) and
Sawyer (1963) and others established that SL algorithms are useful, but did not identify
overwhelming advantages.
The method become much more popular when Robert (1981, 1982) showed how to
combine semi-Lagrangian advection with semi-implicit treatment of the other terms. This
launched an explosion of activity in meteorology which has continued to the present.
Bates and McDonald (1982) made another important advance by introducing “multiply-
upstream” CFL methods. These allowed the use of a time step which far exceeded the
advective limits of conventional explicit schemes.
Temperton and Staniforth (1987) and McDonald and Bates (1987) independently de-
vised two-level schemes of second order time accuracy. (Robert™s method required three
time levels while the earlier two-level algorithms of Bates and McDonald were only ¬rst or-
der.) These allow a doubling of the time step with no loss of accuracy. However, two-level
schemes are somewhat less stable than Robert™s method, so several later articles have de-
scribed slight modi¬cations, such as a “theta” treatment of the implicit schemes, a damping
of divergence, and so on to render such algorithms more stable and robust.
Ritchie (1986) developed the “non-interpolating” version of SL, and then applied it to
a hierarchy of increasingly complicated models in later papers. Smolarkiewicz and Rasch


. 60
( 136 .)