1 50 MONTHLY WEATHER REVIEW Vol. 97, No. 2
UDC 551.509.313:551.511
A TRUNCATIONERRORREDUCINGSCHEMEFORBALANCEDFORECASTMODELS
ARNE SUNDSTROM
International Meteorological Institute, Stockholm, Sweden
ABSTRACT
Various ways to improve the numerical accuracy of solutions t o balanced forecast equations are discussed and
compared. Among these, the most efficient method seem to be a correction operator technique with an associated
space-truncation error of fourth order in As. Results from a number of real data short-range forecasts with this
method are also presented and discussed.
1. INTRODUCTION
In balanced models for routine numerical weather
prediction, commonly used values for the horizontal grid
dist'ance and the time step are As =300400 km, At"0.5-1
hr. With this choice, the numerical errors are dominated
by the space-truncation error, defined as the At-in-
dependent part of the total truncation error. (The re-
maining part is called the time-truncation error.) Its
influence on the solution is of ten noticeable in 24-hr
forecasts and may be considerable for extended forecast
ranges.
As a simple example, the Rossby wave solution
+=-Uy+1C/'eik(z-c2) to the finite-difference form of the
divergent barotropic vorticity equation
a - (W")="(+, W+f> at (1)
with the 5-point approximations
1
m&+ij= a s 2 (1C/$+lj+$i-1 ~+J.~i+l++ii--1"4+ij) (2)
Jw(% P)<j=J++(a, P )i j =i (%+r-~t-lj)(P*,+l 4As l (
-~i j -~)-(~~j +,-~~j -,)(P i +l ~- t -l j
)>
(3)
to the Laplace and Jacobi operators, moves with a phase
speed
C=6---" k2YU-P
k2r+ P
where
Y= COS klAs)+2(1-~0~ k2As)
(kAs)
1
6= k,As sin k,As +k2As sin kzAs
(kA5)
1
kl=k cos cp, kz=k sin cp, and cp is the angle between the
x-axis and the i-direction. This speed may be much less
than its correct value, given by
k'l7-P c=-.
k2+ P
For a grid oriented along the x- and y-axes, assuming
A~=300 km, u=18 m sec-l, p=0.8.10-12 m-2, P=1.62
.lo-" m-l sec-l, the error is 10 percent for a vavelength
X=2n/k of 2600 km, 15 percent for X=2000 km, and 20
percent for X=1800 km (see fig. 1-3). The time truncation
error is here neglected.
For other orientations of the grid the phase-speed errors
are smaller, reaching a minimum (about half of the maxi-
mum value) when the grid is at a 45" angle to the s-y axes.
This directional dependence of the error causes a
spurious distortion of the solution pattern, in addition t o
the reduction of the phase speeds. To avoid this type of
error, slightly more complicated finite-difference approxi-
mations to the Laplace and Jacobi operators have to be
t' tsec
-l r -20
FIGURE 1.-Phase speed for a Rossby-wave solution t o the diver-
gent barotropic vorticity equation. Dotted line is for the exact
equation, solid line for the second-order approximation, and
dashed line for the new fourth-order method. As=300 km,
U= 18m sec-1, q=O.8. 10-12m--2, B= 1.62.10-11 m-1sec-1.
February 1969 Arne Sundstrom 1 51
-30 %
FIGURE 2.-Relative phase-speed error for the second-order method
(solid line) and for the new fourth-order method (dashed line).
used. Of all possible expressions of the values at the point
(i, j ) and the eight adjoining points, the Laplacian with
minimum orientat,ional truncation error is (Miyakoda, 1960)
11100. phose *,.Or)
dmys
1 i
10-
5 -
05
c kbs
0
1 5 r h I .B I aom cmo 3000 z y a 2000 m o i m 3mkm
FIGURE 3.-Time for obtaining a phase error of 180'. Solid line
shows values for the usual second-order method, dashed line
values for the new method.
The corresponding Jacobian is not uniquely defined, and any member of the one-parameter set
may be used. Although there is no restriction on the value
of a, the cases a=O, 1, or 2 are of main interest. For a=l
we have the symmetric conservative Jacobian proposed
by Arakawa (1966) to guarantee computational stability.
(This property has been confirmed by Sundstrom, 1969.)
The approximate independence of the grid orientation
follows from the Taylor series expansions
Unfortunately, the maximum phase-speed error still
remains the same. On the average, the phase-speed error
may thus actually increase, and it is mainly in combination
with one of the following error-reducing methods that the
orientation-independent operators are of importance.
9 . POSSIBLEWAYSTOIMPROVE
THENUMERICALACCURACY OF
THEFINITE-DIFFERENCEAPPROXIMATIONS
1) The simplest way to diminish the truncation errors
is of course simply refining the grid, which should obviously
be done in such a way that the total computation time is
kept as small as possible. Since we cannot a priori tell
whether As and At should both be reduced, we must first
seek for the choice giving the smallest possible total
truncation error (E) within a given total computation time
(T), or alternatively, the smallest T for a given error E.
If the leap-frog method is used for the time differencing, E
is of second order in As and At while T is approximately
proportional to the number of time steps (T/At) , the num-
ber of grid points in the horizontal (S/As2), and the number
of iterations used for inverting the Helmholz operator a t
each time step. For the extrapolated Liebmann method,
the optimal choice of As and At is then the one making the
space truncation error three times as large as the time
truncation error (if the choice is not in conflict with the
stability conditions). The proof is found in the Appendix
(Case 1). It is also shown there that the minimum compu-
152 MONTHLY WEATHER REVIEW Vol. 97, No. 2
tation time is proportional to E-2-to reduce the error by
a factor of two, the computation time increases by a factor
of four. For the Alternating-Direction Implicit (ADI)
method the relation is slightly more favorable and T is
nearly proportional to Anyhow, refining the grid is
only possible if the memory capacity of the computer is
abundant, since the utilized capacity increases as E-1.
2) If the computation is made with two different grid
distances, one may extrapolate the results to zero grid
distance (the Richardson-extrapolation technique, see
Bulirsch and Stoer, 1966). To facilitate the programming,
the larger grid distance (As") is usually chosen as a multiple
of the smaller one (As'). The most favorable case is
As"=2Asr,At"=2At'. The minimum total computation
time is then proportional to E-l, if the extrapolated Lieb-
mann method is used, and nearly to for the AD1
method (see the Appendix, Case 2). In spite of these
favorable results,' the method is presently of limited
interest. The asymptotic error formula is not valid for
kAs larger than about 1.5, giving an upper bound for As
around 400 k m . With Asr=200 km, As"=400 km, the
computation time is five times as large as for the conven-
tional scheme with As=300 km, and the gain in accuracy
hardly justifies such an increase.
3) The two remaining methods are based on approxi-
mations with a local. space truncation error of fourth
order in As (but where the usual leap-frog scheme is used
for the time differencing).
a) The easiest way to obtain such a scheme is inserting
fourth-order approximations to the Laplace and Jacobi
operators in the forecast equation. These finite-difference
operators must use values a t points outside the central
9-point set and may be constructed from the previously
given 5- or 9-point formulae as
4 1 V,"+i,= -V(AS)+~,- 3V2(2As)+ij 3 (8)
J*(% P )i j = 3 J(As)(a, 01ij- ,J(2As)(0r1 P>,j. (9)
and
4 1
Schemes of this type have been thoroughly tested by
Miyakoda (1960, 1962), showing that they accomplish the
desired error reduction as long as the shortest wavelength
exceeds 4As. As shown in the Appendix (Case 3), the total
computation time is proportional to with the opti-
mum choice ofAs and At, provided that the number of
iterations used for inverting the Helmholz operator with
the extrapolated Liebmann method is proportional to As-l.
This is, however, not strictly true, and the computation
time may therefore increase more rapidly. Futhermore,
solution values must be prescribed both a t the boundary
and a t grid points just outside it, which may cause stability
problems. To avoid these difficulties, Miyakoda investi-
gated a scheme with a fourth-order Jacobian but with a
conventional Laplace operator and found that it gave a
noticeable improvement of the phase-speed errors; such a
simplified scheme was introduced in the barotropic fore-
casting routine of the Japan Meteorological Agency in
1960. This change increased the total computation time
by only about 20 percent for short-range forecasts and
even less for longer integration times.
b) By a 'correction-operator technique, showing Some
similarity to Thompson's (1955) inverse-averaging pro-
cedure, it is however possible to find a fourth-order
scheme where the Helmholz operator is only a [)-point
expression, which may be inverted with the usual extra-
polated Liebmann method without any boundary value
problems. The derivation of this method is described in
section 3, and the associated phase-speed error for a
Rossby-wave solution is discussed in section 4 of this
paper. It has been successfully tested for a number of
cases both with a simple channel-flow model and the 24-hr
and 48-hr barotropic forecast routine of the Swedish
Meterological and Hydrological Institute. Some of these
integrations were extended to 72-hr and 96-hr forecast
time. The average reduction of the root-mean-square
error in 24-hr forecasts was about 10 percent. In the 48-hr
forecast, the improvement was less than this value and it
even failed t o appear in the 72-hr forecasts. The reason
for this was found to be the inability of the barotropic
scheme to predict the motion of planetary-scale systems,
causing an erroneous retrogression of the waves. This
retrogression is partly reduced by the influence of the
semiempiric divergence parameter p, but also by the
phase-speed errors, which act in the same direction. With
the new scheme, a larger optimum value of p should
therefore be expected. This value depends upon the size
of the forecast region, the length of the forecast, and the
desired accuracy for different scales of motion. After such
a modification, the correction-operator technique is
probably the most efficient method of reducing the
numerical errors of balanced forecast schemes.
3. DERIVATION OF
THECORRECTION-OPERATORSCHEME
From the Taylor series expansions (6, 7), it follows that
a.nd
(where u is the parameter in the Jacobian (5)) if the func-
tions +, a, and p have derivatives of sufficiently high order.
With the aid of the relations
February 1969 Arne Sundstrom 153
with where
2(1--cos k As) sin k As 7 6=- Y= (k As) k As
and
are easily found to be true for a grid oriented in the x- and
(15) y-directions, and the phase speed is consequently
it may be shown that the equation
gives an approximation to equation (I) with a total truncation
error 0(At2)+O(As4), if the leap-frog method is used For small k AS, this gives c =k z +o ((k ~~)4 ), 8,s
for t h e differencing. A comparison with the k +a
normal scheme based on 9-point operators shows that the expected*
only new features are the two correction operators The Rossby-wave phase speed was computed by this
1--(1-u)- v(i) and 1-(1+u)- v&, the factor 1-p- formula for the case a=O with As=300 km, u=18 m As2 A s 2
12
" l2 sec-',q=0.8.10-12m-2, and P=1.62.10-11m-' sec" (corre- in front Of the 'Perator, and the use Of
Of *, f. For that scheme, One usuallY *,f instead sponding to 450 ]at.); see figure 1. The relative phase-speed
* and '(i)lC. error and the time required for obtaining a phase error of t o be prescribed a t the boundary points. Here,
for other choices of a.
additional 1800 are also shown in figures and 3. For kAs<0.75, the
of kAs, the error is always negative and approximately 25
given percent of the error of the conventional scheme. In this
conditions must be imposed, One if Or and two is alu,aysless than 1 percent, and for larger values
If a # -l J One may assume that '&* is
at the boundary points, when computing the corrected wavelength range, u is much larger than p/k2, and it may
1 - (1 + u )g ~ 12 g)) (V (g) ~+j ). then be shown that a=O gives the smallest phase-speed
error. Although the Jacobian with a=l has the merit of
scheme, no corresponding stability theorem has been
found for the fourth-order method. The value a=O was
therefore used in the experiments, and the results showed
2, If a # the Of the right-hand side for points giving computational stability with the usual second-order next to the boundary involves boundary-point values of
Since they Occur in a term by As2, Only an error
Of Order is committed if they are taken to be no sign of instability provided the initial field was
ciently smooth and At was less than about (1.5u)-lA8.
(which is known at the boundary) as 5. APPENDIX long as these values actually produce a smooth solution. With a sufficiently fine grid, the total computation time
4. PHASE-SPEED ERRORS FOR THE NEW SCHEME r for a numerical integration of the barotropic vorticity
equation is approximately proportional to the number of A Of the of finite-difference time steps (T/At), the number of iterations used for equation (l ) is the Phase- inverting the Laplace or Helmholz operator at each time speed error for a Rossby-wavesolution $=-uY+# e '
step, and the number of grid points in the horizontal
scheme With 5- 0' 9-Point Laplace and Jacobi operators asymptotic rate of convergence is of order ()(AS). The
to the
1 i k (z -c t )
The magnitude of the errors for the usual second-order (S/As2). For the extrapolated Liebmann method, the
was discussed in the Introduction. By direct insertion number of iterations is consequently o(As-l), so that
in (4) and (5), the relations rr. =KlA~-3At". v(i1 ,iXCz-cU ,-k2 Y e t k (z -c l ) (17) For the AD1 method, the rate of convergence is of order
and (log S/As2)" a n d t h u s
J(s", (y, e i L (z -c l ))=-i l c 6 e i k (Z -C 1 )
(18) TADI K ~A s -~ log (&'/AS2) At-'.
1 54 MONTHLY WEATHER REVIEW Vol. 91, No. 2
1) For the conventional scheme, the total truncation If this is in conflict with the stability condition, choose
error E is given by E=AAt2+B As2+higher order terms; the largest possible -=x. At Then
A/B>O. B y t h e method of Lagrange multipliers, the As
minimum value of TL for a given E is found t o be
‘Lmtn = Klx-1(n2+n-2) (&x4+ B2+ C2x2) E-l (28)
(TL)mtn = (g) 16A2K1E-2
3 I2
with F o r t h e AD1 method,
for sufficiently small E, provided that the value for
At/As=(B/3A)’I2 is not in conflict with the stability condi- (T ~~~)~~~=K ~K ~~~[(~~~~+~~~’ log S/Ad2
tions. Otherwise, choose the largest possible At/As=x giving -n3/2 log n2]($~(2B2$Q~~o~1))314~-314 (30)
giving
TLmin =: (A+ B /x ~)~x ~K ~E -~
with
E =: (Ax2+ B) AS^.
Similarly,
with
1 + $4 log S/As2
E= l+log &‘/As2
(23) with K~~~ approximately given by
corresponding to (20, 21). Formulae analogous to (22, 23)
are obtained in the same way.
2 ) To obtain the same formulae for the Richardson
extrapolation method, the fourth-order terms i n t h e
expansion for E must be considered: E=A1At2-A2AP+
BlAs2-B2As4-Q2At2As2+ higher order terms, with A2, B2,
and C2 of the same sign. Assume now
At’At” ””
As!-Asf laK’
Then TL=K~K-’(AS”~+AS’”~)
and
neglecting the higher order terms. If As” is a multiple of
As’, As”=n As’,
and
TL L- Klr-1(n2+n-2) (A2rc4+B2+ C2~2)E-1
E=n2(A2~4+B2+C2~2)A~‘4.
The optimal choice of K is
3 ) For a scheme with a fourth-order space truncation
error but only a second-order time-truncation error
E=AAt2+BAs4+higher order terms, A/B>O,
the optimal time-step is At= - As2 if the extrap-
olated Liebmann method is used. For sufficiently small
As, this is not in conflict with any stability condition of
(3 1 2
the type - <x. Then, At
As -
and
(33)
REFERENCES
Arakawa, A., “Computational Design for Long-Term Numerical
Integration of the Equations of Fluid Motion: I. Two-Dimen-
sional, Incompressible Flow,” JOUTnal of Computational Physics,
Vol. I, No. 1, Academic Press, New York, Aug. 1966, pp. 119-
143.
Bulirsch, R., and Stoer, H. J., “Numerical Treatment of Ordinary
Differential Equations by Extrapolation Methods,” Num-
erische Mathematik, Vol. 8, No. 2, Springer-Verlag, Berlin,
Miyakoda, K., “Numerical Calculations of Laplacian and Jacobian
Using 9 and 25 Gridpoint Systems,” Journal of the Meteoro-
logical Society of Japan, Ser. 2, Vol. 38, No. 2, Apr. 1960,
Miyakoda, K., “A Trial of 500 Hour Barotropic Forecast,” Pro-
ceedings, International Symposium on Numerical Weather
Prediction, Tokyo, Nov. 1960, Tokyo, Mar. 1962, pp. 221-240.
Sundstrom, A., “Stability Theorems for the Barotropic Vorticity
Equation,” Monthly Weather Review, Vol. 97, 1969, (to be
published).
Thompson, P. D., “Reduction of Truncation Errors in the Com-
putation of Geostrophic Advection and Other Jacobians,”
Technical Memorandum No. 1, US. Joint Numerical Weather
Prediction Unit, Washington, D.C., Jan. 1955, 16 pp.
1966, pp. 1-13.
pp. 94-106.
