. 59
( 136 .)


everywhere, and many points have errors of 1/10 or 1/20, even far from the shock.

The crucial data is in the left panel of each graph. When the time step is large, there
is a gap in the departure points in the sense that all trajectories originate at a considerable
distance from the front. This is important because the pseudospectral interpolation, which
for these ¬gures was performed without smoothing or special tricks, becomes increasingly
inaccurate as the front is approached because of the oscillations of Gibbs™ phenomenon.
When the time step is small, some particles originate very close to the front, so the inter-
polation gives poor approximations to the velocities that advect the particles towards the
front from either side. The result is large errors.
SL schemes that use various low order polynomial interpolations suffer the same curse;
close proximity to a jump discontinuity is an ill wind for all approximation schemes. Nu-
merical experiments such as Kuo and Williams (1990) have con¬rmed that SL methods
have their smallest errors as a function of time step for rather large time steps ” a couple
of hours in meteorological applications.
Fig. 14.7 illustrates another odd property of SL schemes. If one plots only the grid point
values of a pseudospectral SL scheme, one ¬nds, for moderate „ , a good approximation to
the exact solution. However, the Fourier series can be summed to evaluate u(x) at arbitrary
points not on the canonical grid. If we do this for many intermediate points, we obtain the
very oscillatory graph shown in the right panel. The pseudospectral solution is much better
on the grid than off.
Finite element methods often exhibit a similar “superconvergence” at the collocation

2 2



1.2 0.6 0.8 1.0
0.6 1.0
Figure 14.7: Same case as the two previous ¬gures, but a “zoom” which shows only part of
the x-interval around the shock. In both panels, the dashed line is the exact solution. On
the left, the solid line-with-circles shows the semi-Lagrangian Fourier solution at the points
of the collocation grid. The solid-plus-crosses on the right is the same solution except that
the Fourier interpolant has been summed at many points between those of the collocation
grid. The semi-Lagrangian solution exhibits “overconvergence”, that is, the solution is
much more accurate AT the collocation points than BETWEEN them.

points. The existence of this phenomenon has been rigorously proved for some ¬nite ele-
ment algorithms, but not as yet for spectral.

14.8 Two-Level SL/SI Algorithms
The three-level SL/SI is simple and robust, but it has drawbacks. The ¬rst is that it re-
quires storing the solution at three different levels. The second is that the implicit terms
are treated with an effective time step of 2„ whereas a standard Crank-Nicholson scheme
would average the unknowns at time levels tn+1 and tn with a time step of „ and an error
four times smaller.
Fortunately, it is possible to replace Robert™s scheme with a two-level method that has
a much smaller proportionality constant. The key idea is to apply the three-level formula
with u’ interpreted as the solution at time level tn rather than tn’1 . The intermediate
velocity at t = (n + 1/2)„ is obtained by extrapolation from previously known values of
the velocity.
The displacements are calculated from

±j = („ /2) uext (xj ’ ±j , tn+1/2 ) (14.45)

where ± is de¬ned to be the displacement over a time interval of „ /2, subscript “ext” de-
notes a quantity which must be extrapolated and j is the index of (spatial) grid points.
The actual time-marching step for Du/Dt + G(u, t) = F (u, t) is identical in form to the
three-level scheme:

u+ + („ /2)G+ = u’ ’ („ /2) G’ + („ /2) Fext

As before, u+ = u(x, tn + „ ), but

≡ F (x ’ ±, tn+1/2 )
u’ ≡ u(x ’ 2 ±, tn ) (14.47)

The most widely used extrapolation formula is

15 10 3
u(x, tn ) ’ u(x, tn’1 ) + u(x, tn’2 )
uext (x, tn+1/2 ) = (14.48)
8 8 8
which Temperton and Staniforth (1987) found preferable to the two-level extrapolation

uext (x, tn+1/2 ) = (3/2)u(x, tn ) ’ (1/2)u(x, tn’1 ) (14.49)

It is possible to extrapolate along trajectories (that is, shifting x positions by the displace-
ment ± as one moves back in time), but Temperton and Staniforth (1987) found it was more
accurate (and cheaper) to extrapolate along grid points.
Extrapolation has a bad effect on stability, as noted above. However, the two-level
method is stable if the extrapolation is done carefully. With complicated physics ” moun-
tains, complicated temperature dependence and so on ” it is sometimes necessary to add
extra damping, such as replacing Crank-Nicholson by the “theta” scheme (which is a lin-
ear combination of Crank-Nicholson with Backwards-Euler that weakly damps high fre-
quency). It is also a good idea to transfer as many terms as possible from the extrapolated,
explicit quantity Fext to the implicitly-treated term G. These empirical ¬xes are problem-
speci¬c and model-speci¬c, so we must defer the details to the papers catalogued in the
bibliographic table below.

The two-level scheme is in theory twice as cheap as the three-level scheme in the sense
that one can achieve the same level of time accuracy with twice as long a time step in the
two-level model as in the older algorithm. However, because the three-level scheme is
more stable and does not need time-extrapolation, it remains competitive with two-level

14.9 Noninterpolating SL Methods and Computational Dif-
The error in interpolating from the evenly spaced grid to the “feet” of trajectories is a source
of computational diffusion. McCalpin (1988) has given a thorough analysis for polynomial
interpolation schemes of various orders, but it would take us too far a¬eld to reproduce his
analysis here. Instead, we shall offer only a heuristic argument.
Consider the simplest case: the pure advection equation in one dimension. The formula
for computing the new unknown in the three-level SL scheme is then

u(x, tn+1 ) = u(x ’ 2±, tn’1 ) (14.50)

Interpolation errors will add a small random component to the displacement ±. Sometimes
the ¬‚uid blob will be advected a little too far, sometimes not far enough, because interpo-
lation errors generate inexact particle velocities and displacements. To the SL scheme, the
foot of a trajectory is not known exactly, but only to within a line segment (or circle or
sphere) which is the size of the average interpolation error. The advected ¬elds are blurred
by this randomness, and blurring is equivalent to a computational damping.
McCalpin™s analysis shows linear interpolation generates ordinary, second derivative
diffusion; cubic interpolation gives a biharmonic (fourth derivative) damping and so on.
For ocean models where the scale of the energy-containing eddies is an order of magni-
tude smaller than in the atmosphere, and the affordable number of grid points-per-eddy-
diameter correspondingly smaller, too, McCalpin argues that, arithmurgically speaking,
this computational diffusion would turn the seas into great viscous pools of honey.
Fig. 14.8 illustrates his point. The semi-Lagrangian solution is damped much more than
that of the standard pseudospectral method even though the explicit, physical viscosity is
the same.
Ritchie (1986) was thus motivated to invent a “non-interpolating” SL algorithm. The
basic idea is very simple. At large Courant number, split the advecting velocity into two

uj = uj,SL + uj,E

where uj,SL is chosen so that the “foot” of the trajectory is exactly on a point of the regular
grid when the particle is advected by uj,SL only. Of course, the actual foot of the trajectory
almost never coincides with a point of the regular grid, but this is not a serious limitation
since the SL algorithm is really just a transformation to a moving coordinate. There is
no law of nature or arithmurgy that says that the velocity of the coordinate change must
exactly match that of the ¬‚uid.
However, Ritchie™s modi¬ed transformation eliminates only part of the advection, so
the remainder must be treated by the standard pseudospectral method, just as if we used
an Eulerian (¬xed grid/trajectory-ignoring) coordinate system throughout. It is important
to choose uj,SL to advect to the grid point nearest the actual foot of the trajectory so that
uj,E is as small as possible. The residual Courant number is then always less than one so

1 1




-1 0.5
-2 0 2 2 2.5 3
x x

Figure 14.8: Comparison of the semi-Lagrangian/semi-implicit Fourier method (thick
lower curve) with the standard Fourier pseudospectra method (upper curve with circles)
for Burgers™ equation, ut + uuxx = ν uxx where ν = h/2 and h = π/32 is the grid spacing.
The SL time step was „ = 1/2; that for the standard method was determined by an adap-
tive 4th/5th order Runge-Kutta subroutine. The solution from u(x, 0) = sin(x) “breaks”
at t = 1, that is, would develop an in¬nite slope at that time were it not for the viscosity.
The solutions are shown at t = 2. In spite of the long time step, which tends to slow the
diffusion of the νuxx term (which is treated implicitly), the SL/SI algorithm is much more
strongly damped than the standard pseudospectral method, which treats the viscosity hon-
estly by using a short time step. The SL/SI damping comes primarily from interpolation

that the integration is stable with very long time steps where the residual Courant number

CoE ≡ max(|uE |) (14.52)
where hmin is the smallest distance between any two points on the regular pseudospectral
Ritchie™s later papers (1987, 1988, 1991) show that the “non-interpolating” SL method
works just ¬ne with both spectral and ¬nite difference spatial discretizations, even in three-
dimensional models. However, his idea has only partially replaced older SL algorithms.
There are two reasons. First, treating some of the advection by a non-SL scheme brings
back some of the disadvantages and poor front resolution of conventional methods. Sec-
ond, for atmospheric forecasting models, the resolution is suf¬ciently high and the length
of the integration is suf¬ciently short that damping is not a problem. However, climate
models, which are run at lower spatial resolution for long time integrations, may bene¬t
from non-interpolating SL.
For example, his comparisons with non-SL models showed that at a resolution of “T106”,
then the standard for forecasting, both interpolating and non-interpolating SL models lost
only about 1% of the initial energy. This was actually better than the standard non-SL

spherical harmonics model of the European Centre for Medium-Range Weather Forecast-
ing, which lost about 3% of its energy due to an (unphysical) biharmonic damping which
is included to suppress small scale noise and give better front resolution. The SL models,
in spite of their inherent implicit computational dissipation, were actually less damping
than the non-SL algorithm because an explicit arti¬cial viscosity was unnecessary for the
SL schemes.
At lower resolution, as in a climate model or an ocean model, Robert™s SL algorithm
is more disspative because the interpolation errors worsen as the grid coarsens. At T63,
Ritchie found that his non-interpolating model lost only 1% energy, but the interpolating
SL lost 5%.
Smolarkiewicz and Rasch (1990) have generalized non-interpolating SL schemes. They
show that their generalization can be applied to many existing, non-SL algorithms to relax
time step limitations while retaining good properties of the methods such as that of being
“positive-de¬nite monotone”.
The conclusion is that non-interpolating SL may be very useful for ocean and climate
modelling, or whenever the resolution is suf¬ciently poor that interpolation errors create
too much damping. However, the errors are a function of the order of the interpolation as
well as the grid spacing. If one uses high order interpolation, such as fast spectral off-grid
interpolation, or even quintic polynomial interpolation, instead of the cubic interpolation
which is now the standard in atmospheric models, the computational diffusion is greatly
reduced. Thus, it is only in the context of speci¬c problems with speci¬c grid size and a
particular interpolation order that one can intelligently answer the question: Are interpo-
lating or non-interpolating schemes better?

14.10 Off-Grid Interpolation: Shape-Preserving, Conserva-
tion and All That
14.10.1 Off-Grid Interpolation: Generalities
Off-grid interpolation is rather expensive, if done in the obvious way. For spectral methods,
off-grid interpolation cannot be done by the Fast Fourier Transform because the departure
points are spaced irregularly between the points of the canonical grid. Summing an N -term
spectral series at each of N points costs O(N 2 ) operations.
There are two subproblems which require off-grid interpolation: (i) calculating the dis-
placements and (ii) calculating properties at departure points. The interpolation of mid-
trajectory velocities can be done by a lower order method than interpolating velocities,
temperatures, ozone concentration and so on at the departure points. The reason is that
the displacement is the velocity multiplied by the time step, and so errors in the displace-
ment are smaller than those in u itself.
It is much more important to interpolate accurately when calculating the properties of
a ¬‚uid blob at its departure point. The reason is that errors in this step are not multiplied
by the time step, but rather are errors in the answer itself. An error of a tenth of a degree
in interpolating the temperature at the departure point will be carried through the entire
remainder of the computation.
The third comment is that SL schemes do not automatically conserve either mass or
energy, nor preserve the positive de¬niteness of densities like mass and water vapor which
are always non-negative. It is possible to modify SL procedures to recover many of these


. 59
( 136 .)