January 1971 15
UDC 651. M)s. 313:551.611.2
NUMERICAL INTEGRATION OF THE PRIMITIVE EQUATIONS WITH A FLOATING SET
OF COMPUTATION POINTS: EXPERIMENTS WITH A BAROTROPIC GLOBAL MODEL
FEDOR MESINGER 2
Department of Mefeorology, University of California, Los Angeles, Calif.
ABSTRACT
A Lagrangean-type numerical forecasting method is developed in which the computational (grid) points are ad-
vected by the wind and the necessary space derivatives (in the pressure gradient terms, for example) are computed
using the values of the variables a t all the computation points that at the particular moment are within a prescribed
distance of the point for which the computation is done. In this way, the forecasting problem reduces t o solving the
ordinary differential equations of motion and thermodynamics for each computation point, instead of solving the
partial differential equations in the Eulerian or classical Lagrangean way. The method has some advantages over the
conventional Eulerian scheme: simplicity (there are no advection terms), lack of computational dispersion in the,
advection terms and therefore better simulation of atmospheric advection and deformation effects, very little incon-
venience due to the spherical shape of the earth, and the possibility for a variable space resolution if desired. On the
other hand, some artificial smoothing may be necessary, and it may be difficult (or impossible) to conserve the global
integrals of certain quantities.
A more detailed discussion of the differencing scheme used for the time integration is included in a separate sec-
tion. This is the scheme obtained by linear extrapolation of computed time derivatives to a time value of to + aAt where
to is the value of time at the beginning of the considered time step At and where a is a parameter that can be used to
control the properties of the scheme. When choosing a value of a between % and 1, a scheme is obtained that damps the
high-frequency motions, in a similar way as the Matsuno scheme, but needs somewhat less computer time and, with
the same damping intensity, has a higher accuracy for low-frequency meteorologically significant motions.
Using the described method, a 4-day experimental forecast has been made, starting with a stationary Haurwita-
Neamtan solution, for a primitive equation, global, and homogeneous model. The final geopotential height map
showed no visible phase errors and only a modest accumulation of truncation errors and effects of numerical smoothing
mechanisms. Two shorter experiments have also been made to analyze the effects of space resolution and damping in
the process of time differencing. It is felt that the experimental results strongly encourage further testing and inves-
tigation of the proposed method.
1.
2.
3.
4.
5.
6.
7.
CONTENTS
15
16
18
21
21
22
22
22
24
27
28
28
1. INTRODUCTION AND FORMULATION
OF THE METHOD
Ever since the beginning of the development of modern
meteorology, Lagrangean methods have had a strong
appeal to theoretical meteorologists. The laws of motion
and thermodynamics relate the forces and heating to the
individual derivatives of velocity components and po-
1 UCLA Department of Meteorology Contribution No. 180
2 Present affiliation: Department of Meteorology, University of Belgmde, Yugoslsvh
tential temperature; and in this way, the Lagrangean
approach is the most fundamental one. Nevertheless, due
to practical reasons, during the last two decades of
extensive research in numerical prediction and atmos-
pheric general circulation experimentation, the Eulerian
method of formulating and solving the relevant atmos-
pheric equations has been used almost exclusively. There
has been only a modest number of attempfs to use a
Lagrangean or a quasi-Lagrangean technique. Graphical
forecasting methods are essentially Lagrangean (Fjgrtoft
1952, 1955); however, their use has diminished with the
easier availability of fast digital computers. Anf example
in the paper by Welander (1955) illustrates the major
difficulty to be expected in trying to perform a numerical
forecast with the use of the classical Lagrangean appradch :
a chessboard set of fluid elements soon deforms into a
highly irregular pattern, unsuitable for space differentia-
tion along the originally straight material lines. Thus,
frequent redefinition of the material coordinate lines
would be necessary. Besides, the computational complexity
of this approach appears to be considerable; see, for
example, the analysis by Wiin-Nielsen (1959).
To avoid these problems or for other reasons, a number
of investigators have used the Lagrangean formulation of
16 MONTHLY WEATHER REVIEW Vol. 99, No. 1
the advection terms only (qkland 1962, Krishnamurti
1962, Leith 1965). In this way, an advective change at
a space-fixed grid point and in a particular time step
is essentially obtained by finding the location of the
parcel that mill in this time step be carried by the flow
to the grid point and by computing the corresponding
individual value through space interpolation. Since by
this procedure a space interpolated value is advected a t
each grid point and time step, i t is not clear whether the
method has any Lagrangean-type advantages, It does
represent a simple possibility to avoid the “nonlinear”
computational instability; however, it is questionable
whether the gain in simplicity, as compared t o the more
involved application of the Arakawa (1966) concept of
maintaining the integral constraints of physical im-
portance, will often be of more value than the advantages
of the latter alternative. In fact, it has likely been to a
rather large estent due to the wide application of the
Arakama conservation schemes that the remarka,ble
achievements in atmospheric simulation experiments have
been accomplished in recent years. Besides, other less
refined possibilities to overcome t’he nonlinear instability
problem also exist.
Despite this outstanding performance of the present
day Eulerian methods, they do contain a number of
intrinsic deficiencies. It appears convenient to mention
them briefly:
No attempt can be made to conserve, within advection terms,
individual properties following the motion of air particles. That is,
n false dispersion of the wave components of the flow is produced
t )y the finite-difference formulation of the advection terms (Wurtelc
1961, Matsuno 1966b). The error in phase speed of individual waves
increases cata5trophicallg when the wavelength approaches the
smallest possible wavelength resolvable by the grid. This precludes
the formation of sharp gradients of an advected quantity, such as
the gradients observed in frontal and intertropical convergence
zones, and smooths out these zones in case that some initially exist.
Another consequencc of the Eulerian formulation of advection terms
is that the nicntioned process of elongation, tangling, and subse-
quent mixing of fluid elements is not realistically sirnulated (DjuriE
1966).
It seems futile in Eulerian methods to try alleviating this wave
dispersion difficulty by having a higher grid resolution in zones of
sharp gradients. These zones are usually moving and after some
time would so leave the regions of high resolution. Besides, a space
discontinuity in grid resolution results in false reflection of waves
at the boundary separating the two regions.
The task of defining the locations of grid points on spherical
earth so as to accomplish a quasi-homogeneous space resolution,
though feasible in a number of straight-forward ways, is an uncom-
fortable one, leading to laborious computational schemes (Kurihara
and Holloway 1967, Sadourny- et al. 1968).
A Eulerian model accepts information a t predetermined grid-
point locations only, while the observations are and will be made at
relatively randomly spaced points. This restriction may be of some
disadvantage, especially when, as expected, a crucial part of the
atmospheric observations becomes continuous in time.
In trying to remedy the leading of the mentioned
shortcomings, i t was proposed by Charney (1966) to let the
grid points move, as if imbedded in some hypothetical
medium, in such a way that the grid spacing constantly
adjusts to the gradient of the quantity one wishes to
.
resolve. In this way, one would have two simultaneous
flows: the physical flow of the fluid and the flow of this
hypo thetical medium. Charney rejected the possibility of
allowing the two fluids to coincide since, as described
already, the grid lines moving with the physical flow would
soon become highly distorted and useless for the com-
putation of space derivatives. It is suggested here to try the
other alternative: let the computation points move with
the physical fluid and, instead, reject the requirement
for the continuity of their relative position. The necessary
space derivatives have then to be computed by using the
values of the variables at all the computation points that
a t the particular moment happen to be in some way
surrounding the point for which the computation is done.
Even if they initially are not, after some time, say a
number of days, these surrounding points will be dis-
tributed essentially a t random; it is conceivable that
this may result in a significant increase in the number of
computation points necessary to achieve some required
accuracy in space differencing. But, in return, the
prognostic equations reduce to a set of ordinary, instead of
partial, differential equations, with a considerable gain in
simplicity. Besides, no errors are produced by the finite-
difference treatment of the advection terins since there
are no advection terms.
The proposed method has been tested using a global
primitive equations model for a homogeneous and incom-
pressible atmosphere. Following sections describe the
technical details of the computational procedure as well
as the results of a number of performed experiments.
Finally, a comparison of the properties of the method with
those of the conventional Euleriaii schemes will be given.
2. SPACE DIFFERENCING
Throughout this study, we shall deal with a homo-
geneous, incompressible, and frictionless atmospheric
model. Although simple, such a niodel contains a sig-
nificant part of the essential features of the large-scale
atmospheric motions ; see, for example, the discussion in
the paper by Arakawa (1970). For these reasons, i t has
been used many times to test the performance of the
numerical finite-diff erence schemes. The most recent
reference describing such a use is probably the one by
Grammeltvedt (1969) ; it contains a number of the previous
ones. As is customary, the pressure change with height is
assumed in the present study to be given by the hydro-
static equation. The pressure force is then not a function
of height; and if the initial wind is independent of height,
i t will have to remain so. We shall assume this to be the
case. We shall further consider earth to be a perfect
sphere and choose a A, 0, r spherical coordinate system
with, for later convenience, an arbitrary orientation of the
O= f 7r/2 axis. Having a shallow atmosphere, we replace
the radius T with a constant value a, and accordingly
(Phillips 1966) neglect the vertical velocity terms in the
equations of horizontal motion. The governing equations
Fedor Mesinger 17 January 1971
then become
and
Here, U and V are the horizontal velocity components in
the directions of constant 0 and constant A, respectively;
+ is the geopotential of the free surface; t is time; time
derivatives are individual rates of change, following the
motion of the fluid; fi is the magnitude of the earth’s
angular velocity; and p is latitude. We shall further use
h for longitude and denote by u and v the eastward and
northward components of the horizontal velocity,
respectively .
The problem of space differencing now consists of
computing the four space derivatives in eq (1). We want to
do this by making use of the known values of u, v, and +
a t the point for which the computation should be per-
formed, which we shall call the “reference point,” and of
the values of these variables a t a number of surrounding
computation points, which we shall call “neighbors.”
The space distribution of these points is considered to be
arbitrary, except for assuming that they are not or-
ganized in some special inconvenient way (along a line,
for instance). This distribution is given by the known
values of their geographical coordinates h , cp. For a de-
fined orientation of the A, 8 system, the components u, v
and coordinates h , cp can be transformed readily, of course,
into components U, V and coordinates A, e, if this definition
is such as to make this necessary.
Now, we may expect that the property of being a
neighbor has to be a reciprocal one. If it were not, we
would have a situation in which the values of dependent
variables a t a computation point affect the time change
of those a t another one, without the opposite being true.
It would appear that this should lead to an instability of
the computation; indeed, some test runs support this con-
clusion. The simplest way of making the property of
being a neighbor reciprocal is to take as neighbors all the
computation points which a t that particular moment
happen to be within some prescribed distance p of the
considered reference point ; this possibility has been used
in the present study. Now, a method for the computation
of space derivatives will always require a certain minimal
number of neighbors; to make sure that we always have a t
least that many neighbors in a computation with a
randomlike distribution of points, we then have to make
p such that the average number of neighbors is consider-
ably greater than this minimal one. I n this situation, the
method of least squares offers a straightforward way to
compute the space derivatives, and it has been used here.
The accuracy of the least-squares fitting, if performed
in terms of the geographical coordinates A, cp, would ob-
viously break down when approaching the two Poles.
For avoiding this, the least-squares fitting was done in
the coordinate system A, 0, which had its origin always in
the considered reference point, and e= &~/2 axis in its
meridional plane, with 0 increasing northward. This
choice of origin has an additional benefit of simplifying
eq (l ), in that the three tan 0 terms disappear, if computed
forward in space.
The necessary coordinate transformation equations can
be written in the form
A = arctan (Y/X)
e=arctan(Z/JX2+ y2)
X= cos cpo cos cp cos (h - h,) + sin cpo sin cp,
and
where
Y=coscp sin(h--Xo),
Z= -sin cpo cos cp cos(X--Xo) + cos cpo sin cp,
0
and ho, cpo are the geographical coordinates of the refer-
ence point. Except along the meridian through the refer-
ence point, the coordinate lines of the A, 8 system, in
general, make an angle with respect to those of the h, cp
system. Denoting this angle by P and having it increase
in the positive direction, we obtain from eq (2)
sin p=-sin cpo sin(h-Xo)/JX2+Y2
and (3 1
cos p=[cos cpo cos cp + sin cpo sin cp COS(~-A~)~/JX~+Y~.
The velocity components U, V can then be computed
making use of the relations
U=u cos p+v sin P
and . (4 1
V=-u sin B+v cos 0.
After having in this way performed the transformation
of coordinates A, cp and components u, v into coordinates
A, e and components U, V, the least-squares fitting was
done in the experiments reported here by defining
Here, x stands for any of the variables, U, V, and 4, and
18 MONTHLY WEATHER REVIEW Vol. 99, No. 1
xo is the value of x at the reference point; the five coeffi-
cients are computed so as to make 2 differ as little as possi-
ble, in the least-squares sense, from the values of x a t the
individual neighbors. Having, as described, more than five
neighbors, the two will have to be, in general, different.
If we wish, we can prescribe some weights w to the squared
differences a t the individual neighbors and compute the
coefficients so as to minimize the sum
s=2 Wt(X,-;;,)2. i =l
The subscripts here refer to the values of the variables a t
particular neighbors, the total number of neighbors being
n. In the experiments that will be described here, the
same weight has been prescribed always to all the neigh-
bors. It might be possible, however, to optimize the
weights of individual neighbors in some way, say as a
function of their space distribution or the reliability of the
x values a t particular neighbors.
NOW, obviously, we assume
(7) a) a) au av an + ae- --&A , --B+, aK=A,, and - ae --=€ I,;
the computation of space derivatives in eq (1) is in this
way defined.
A few words could be added about the choice of the
neighbor-defining distance p and the search for neighbors.
Preliminary experiments with the present method have
shown that a tendency for the appearance of computa-
tional noise exists, this tendency becoming catastrophic
with the decrease in p. This noise may be partly a result
of the constant change of neighbors; this change must be
associated with some disruption of the balance in the
fields. Another source of the noise should be expected in
the errors in space differencing that are produced by in-
convenient space distribution of the neighbors. For
keeping this noise tendency a t an acceptably low level,
it was found necessary to have a fairly large average
number of neighbors, Z. The experiments reported here
have all been done with the value of
p=a, arccos 1-2 - ( 3
such as to make E equal to 30; here, N i s the total number
of computation points. With such a large value for E ,
no need was felt to check whether the particular values
of n are, a t each reference point and time step, always
greater than or equal to 5.
For enabling a fast search for the neighbors, the globe
was divided into spherical grid boxes; and a t the beginning
of the computation and after each displacement of the
computation points, lists of addresses have been made of
all the points found in particular boxes. For a specific
reference point, then only the points found in boxes that
are within the distance p of the reference point have been
checked for the possibility of being neighbors. This
procedure appeared to be fairly fast, in the sense that it
did not consume a major portion of the total computation
time.
3. TIME DIFFERENCING: THREE-LEVEL
EXTRAPOLATION SCHEME
It has earlier been pointed out by Matsuno (1966a)
that a discontinuity in the grid size of a Eulerian grid can
result in a false generation of high-frequency waves. As
pointed out, this apparently takes place also in our case of
an irregular distribution of computation points. For such
situations, Matsuno recommended use of a computational
scheme that filters out high-frequency oscillations during
the process of time integration. His first-order accuracy
scheme (Matsuno scheme) has since been used extensively
in some of the general circulation studies (Mintz 1965,
Arakawa et al. 1969). It was pointed out by Lilly (1965)
that this scheme may possibly cause a considerable loss in
the energy of the system; on the other hand, for the
purpose of initialization for the primitive equations model,
some authors felt the need for a scheme with even more
powerful damping of the high-frequency motions (Nitta
1969). Use of such high-frequency damping schemes is
likely to become much more extensive with the expected
advent of the time-continuous observations: we may then
wish to use these schemes to speed the adjustment of the
imbalances introduced by continuous input of the data
into a numerical €orecast. In view of this general interest
and possible advantage in having a variable-damping
scheme, a scheme that was analyzed for use in the present
model and can have a variable damping will be described
here in more detail. This is the scheme obtained by per-
forming a linear time extrapolation of computed time
derivatives to a time value of to + aAt where to is the value
of time at the beginning of the considered time step At and
a is a parameter that we can use to control the properties
of the scheme.
We shall write a particular one out of eq (1) in the
general form
ax -=j( a t x, t ), x =x( t ) .
In writing this, we have disregarded the coupling within
the system (1) in exchange for the possibility of considering
x to be complex and thus having eq (9) stand for a pair of
eq (1). If we now use an integer subscript to denote the
number of time steps elapsed, we can write the considered
finite-difference version of eq (9) as
For brevity, we shall sometimes call approximation (10)
January 1971 Fedor Mesinger 19
TLE (three-level extrapolation) scheme. It reduces to the
Euler (forward) scheme when a=O and to the simplified
Adams-Bashforth scheme when a=jh. When a= 1, it
represents a simulation of the backward difference scheme.
If we define as the truncation error e of the considered
scheme, the remainder on the right side of eq (10) when
xr+, are substituted by their true values obtained
through Taylor series expansion, we have
~=-(~-~)f:(At)"-(g+a)f:'(At)~+ . . . .
Thus, approximation (10) is of first-order accuracy, unless
a=%, when its accuracy is of the second order.
To analyze the behavior of the error over a period of
time, however, we have to prescribe first the function
f(x, t ). It is of particular interest to consider the case of
the oscillation equation
ax . z =z w x
since it describes the linearized gravity-inertial motions
present within eq (1). We shall do so throughout the
remainder of this section.
We tentatively write our finite-difference solution of
eq (11) in the form
X,=&XT. (12)
Combining eq (10-12), we then see that eq (12) is the
required finite-difference solution of eq (1 1) , provided that
(-A+B)'/'I X=+[1 f +(A +B y 2 ] +i*[ (1 +a )p i
~
211-Ul
1-a
(13)
where
p=wAt,
A=3[1- (l+~)'p'],
and B=$[1+2( 1 -6a+az)p2+( 1 +a) 4 p 4 ] l/Z*
In eq (13), the upper signs correspond to the "physical,"
and the lower correspond to the spurious "computational
mode" that appears through the use of the three time
levels differencing scheme. Now, we can also write
h=Iileio. (14)
Since the true solution of eq (11) can be written as
we see it appropriate to refer, as customary, to 1x1 as the
amplification factor and to O/wAt as the relative phase
change (per time step) of the finite-difference solution.
For stability, it is now necessary to have both values of
1x1 5 1; for accuracy, it is desirable to have 1x1 and O/wAt
of the physical mode close to unity and to have 1x1 of the
computational mode as small as possible.
We shall first investigate the stability and damping
properties of the considered TLE scheme. From eq (13),
we obtain
'
I X I = &[ 1 + ( 1 + a) 'p2 + B f (A + B) '1'
1--a2
5 --~p(-A+B)"']'~'. (15) 11-4
The right side of eq (15) is equal to unity when
w A t =-(- 1 1-2a ---) "' .
a l+%
.. The TLE scheme is, therefore, stable for values of wAt
equal to or smaller than the right side of eq (16). This
right side has its maximum value
(wAt) ,,,az=42(5&- 11) = 0.6006
when
a= (1 +dg) =0.8090. (17)
Thus, as far as the extent of the stability range is con-
cerned, eq (17) gives the optimum value of the parameter a.
Amplification factors given by eq (15) for this optimum
stability and three other particular values of a are plotted
against wAt in figure 1. It demonstrates a favorable
damping of the computational mode for values of a<l.
However, the computational and physical mode curves
meet for the value of a=l and flip over each other for
larger values of a, the computational mode thereby be-
coming unstable and limiting the stability range of @At.
Still, the stability range decreases rather slowly with the
increase in a, and large values of a may thus be useful if
an especially powerful damping of the high-frequency
motions is desired.
I t is of interest to find the location of the minimum
value of the greater of the amplification factors since the
scheme has the strongest damping at this value of wAt, and
one may wish to choose At so as not to exceed this maxi-
mum damping wAt value. First, for a <l , one finds from
'-i
wAt
FIGURE 1.-Amplification factors of scheme (lo), plotted against
d t , for values of a= 0.50 (simplified Adam-Bashforth scheme),
a =(%)(l +a ) (value giving the largest stable region of oAl),
a= 1.00 (a simulated backward scheme), and a= 1.25.
20 MONTHLY WEATHER REVIEW Vol. 99, No. 1
eq (15) that the two 1x1 = constant lines tangent to the
1x1 curves are
Thus, these are the extreme values of 1x1, the upper sign
again referring to the physical mode. They occur at
- 17 +4pa- 29a2+2a3 F 2'13(3 - 7a+2a2)( 1 l I 2 .
(19)
Hence, for a<l, eq (19) with upper signs gives the maxi-
mum damping wAt value as a function of a. It has a
maximum of about 0.5132 a t a=%.$.
For a>l , maximum damping occurs at the intersection
of the two 1x1 curves. This happens at
wAt=l/(l+a). (20)
3 - 2~ & 23/2( 1 - a ) '/'
1
Both expressions (19) and (20) are of course valid in the
limiting case when a= 1.
A summary of the dcpendence of considered wAt values
on a is given in figure 2. Thc figurc shox-s the extent of
the stability range, eq (16), and wAt values at which the
considered TLE scheme has a maximum damping, eq (19)
with upper signs for a<l and eq (20) for ~2 1 .
To discuss now the accuracy of the finite-difference
solution, me also have to consider the phase change, per
time step, of its physical mode. Using the subscript 1 to
denote the values that pertain to the physical mode, we
can write this phase change as
-+ (2+ (A+ B ) ' 91. (21)
We want t o compare 8' with the phase change, per time
step, of the true solution that is equal to wAt. For that pur-
pose, ratios of the two phase changes as functions of wAt
are plotted in figure 3. This is done for a number of values
of the parameter a and for comparison also for the first-
order Matsuno scheme.
Normally, however, me will only be concerned with the
accuracy of the finite-diff erence solution for small values
of wAt. Thus, we will require IX1( and t?,/wAt to be as close
as possible to unity only when wAt is close t o zero. Figures
1 and 3, while giving a general description of the behavior
of and O,/wAt, are not very helpful in estimating how
well this requirement is satisfied. To this end, considera-
tion of their series expansions in terms of powers of
p=wAt is more convenient; for these, me obtain
lh1~=1+3(1-2a)p2+. . . (22)
STABILITY
w A t
0
FIGURE Z-Curves showing the extent of the stability range of
scheme (10) and oAL values a t which a maximum damping occurs
as functions of the parameter a.
1.6
I .4
z 3
-- 1.2
\
1.0
025 050 0.75 IO0
wnt
FIGURE 3.-Ratios of the phase changes per time step of Anite-
difference solutions and the true solution plotted against o A t .
The solid curves refer to the physical mode of scheme (10) and
values of the parameter a written by t'ne curves; the dashed
curve refers to the first-order Matsuno scheme.
and
(23) e,& 1 - 4 (1 - 6a+3a2)p2+ . . . .
Thus, the amplitude error has a minimum in case of the
Adams-Bashforth scheme when a= )i. The phase error,
on the other hand, has a minimum when
and a maximum in case of the simulated backward
scheme when a= 1.
An appropriate estimate of the accuracy of the TLE
scheme appears to be its comparison with the first-order
Matsuno scheme. The amplitude and relative phase
change expansions for this Matsuno scheme are
Ih,l=l-$pZf . . .
and
eMlp=i+{p2+ . . . .
As could have been expected, these leading terms are
equal t o the values that thc corresponding terms in eq
(22) and (23) take in casc of the simulated backward
scheme. We see that, p being the same and small, the
amplitude error of the TLE scheme is less than that of
the Matsuno scheme for all values of O<a<l; its phase
January 1971 Fedor Mesinger 21
error, moreover, is less than the Matsuiio scheme error
for all reasonable values of a except for a=l. Of course,
p=wAt would not normally be the same; the optimum
choice of at appcars to be such as to obtain maximuni
damping a t the highcst expected frequency. Since the
minimum amplification factor of the Matsuno scheme,
J3/2 ~0 .8 6 6 0 , occurs at wAt= l@, this choice mould
imply having At of the TLE scheme, AtTLE, shorter than
the At of thc Matsuno scheme, AtM, so that
AtTLE At,/&
if a is chosen between about 0.7 and 1. Still, since thc
TLE scheme needs only about half as much computation
timc per time step as does the Matsuno scheme, compu-
tation with thc TLE scheme and this choice of time stcps
would be about 30 percent faster. Having chosen the
valuc of AtTtE, the intensity of damping of thc highest
cxpcctcd frequencies can still bc adjusted within fairly
wide limits by a suitable choice of a ; the valuc of a=l
would, for instance, result in about two to three times
stronger damping of the highest frequencies than with
the Rlatsuno scheme, while the valuc giren by eq (17)
xi-ould already result in a damping slightly weaker Lhan
Matsuno’s. At the same time, for low-frequency motions,
both amplitude and phase errors of the TLE scheme
would be significantly less than the corresponding
h’latsuiio scheme errors. This now includes also thc a= 1
case, since the shortening of the time step leads to an
additional decrcase in total truncation errors, due to the
curvcd shape of thc JX,1 and O,/wAt curves.
A possible disadvantage of the present TLE scheme
may be its rcquiremeiit for storage of time derivatives in
the computer code of the model. However, in the present
floating computation points method, this is offset largely
by saving the need for storage of an extra set of point
coordinates, what would have been required when using a
stable or neutral tmo-levels scheme.
Use of a three time levels schemc, of course, is associated
with a need for a special starting procedure. However,
since for reasonable values of a the computational mode of
the TLE scheme is strongly damped, no dficulties should
be expected if the simplest forn-ard scheme is used for
the first time step. This was, therefore, done in the
experiments reported here.
Finally, it may be worth mentioning that, obviously, a
difference exists between the application of the stability
analysis to the present Lagrangean method and its usual
application to Eulerian methods. The frequency w now
refers to the individual changes and not to the local
changes as in the Eulerian models.
4. REMAINING DETAILS OF THE METHOD
SUPPRESSION OF T H E SPACE N O I S E
Besides the use of the described high-frequency damping
time-diff erencing scheme, two other noise-suppressing
mechanisms were built into the model. While the purpose
of the damping version of the TLE scheme was to suppress
the small-scale features in time changes of the dependent
variables, the purpose of the remaining two mechanisms
was to deal with the apparently more fundamental
problem of the growth of the small-scale features in space
changes of these variables. This spacc noise in the fields
of dependent variables, as discussed at the cnd of section
2, appears to be for the most part a result of the accumu-
lation of random truncation errors in space differencing.
In Eulerian models, accumulation of such errors is nor-
mally prevented by the interaction of the neighboring
grid points; each irregularity then behaves as a small-
scale gravity-inertial wavc and is thus dispersed and
possibly also damped in the process of time differencing.
This may not bc possible in the present method since the
space differencing by the least-squares fitting to a large
number of neighbors makes the interaction of the neigh-
boring cornputation points extremely meak-at least when
the weight of the neighbors is not madc to depend on their
distance from the reference point. With thc usual mech-
anism for thc suppression of space noise thus mostly
lacking, an alternative procedure had to be constructed.
First, a constraint \i-as imposed 011 the magnitudes of the
components of gradients of dependent variables, as com-
puted using the values of the variables a t reference points
and their nearest neighbors. These magnitudes were not
allou ed to surpass some arbitrarily prescribed very large
d u e s , believed to be impossible in an error-free develop-
ment of the flo\\-. A check on their actual d u e s was
performed at each computation point and time step! before
the computation of space deriviatives; for the case that
this constraint is rioltLted, instructions \\-ere provided to
change the dependent variables a t the t x o points so as to
reduce the gradient component to the allowed magnitude,
\\ ithout changing the average V ~I L W of the variables. This,
in fact, represented an enforcemelit of the requirement that
the dependent variables a t two computation points ap-
proach each other in case that the distance between these
points happens to approach zero. However, the experi-
ments reported here have all been started with initial
point distributioiis such that 110 two points \\-ere close to
each other; as a result, only rarely did any tn-o points
approach each other very closely, and the described nearest
neighbor smoothing most likely, in these experiments, has
never been performed.
The second, less exclusive, smoothing mechanism con-
sisted of adding a diffusion simulating term to the right
sides of each of the three eq (1). This term had the form
.;.(\ - k d (X - X)
where k , is a constant and the circumflex represenh an
average over all the neighbors that are within some pre-
scribed distance p, of the reference point, this distance
being less than the distance p defined by eq (8). I n present
22 MONTHLY WEATHER REVIEW Vol. 99, No. 1
experiments, the distance Pa was defined so as to cover an
area equal to one-third the area covered by the neighbor-
defining distance p , and the constant kd was defined so as
to give unity if multiplied by the time step of 4 hr. The
effect of the terms in (25) is expected to be mostly smooth-
ing of the space irregularities in the x fields that are formed
by the individual values of x. Computation of the x changes
due to the terms in (25) was not done by using the TLE
time-differencing scheme, but rather with simple forward
steps.
TRAJECTORY COMPUTATIONS
The most straightforward approach to the computation
of trajectories of computation points is to consider that it
consists of solving the system of equations
ax
dt a, cos Q -=u(x, Q, t ),
(26)
It has previously been shown (Djurib 1961, Mesinger
1965) that the truncation errors of trajectory computations
are rather small, even with time steps of the order of 1 hr
or more. The present method might have high demands
with respect to the accuracy of the trajectories; however,
since maintenance of linear stability requires much
shorter time steps, it appears that even the demands for
an extremely high accuracy of the trajectories should be
.satisfied with any reasonable scheme. One. problem,
though, is to avoid the polar singularities in eq (26). In
computations on a limited domain, this has customarily
been accomplished by transforming eq (26) into its
counterpart on the image surface of a conformal projec-
tion. This seems to be impractical in a global model
since no single conformal projection can be used for the
entire sphere. However, a simple alternative exists and
has been used here: trajectories can be computed using
spherical geometry, assuming that each trajectory seg-
ment is an arc of a great circle. This procedure is especi-
ally comfortable in the present prognostic calculation
where it involves no space translation of the velocity
vector.
To display the necessary formulas, let us first observe
that the displacement of the trajectory points is centered
in time if performed using the wind vector
Here, dv/dt is the time-extrapolated value of the accelera-
tion as described in the previous section. The angular
displacement of a trajectory point in the time step At is
then
a=Iv,+tIAt/a,.
We now consider the spherical triangle with two of its
vertices a t the positions of this trajectory point at time
levels T and T+ 1 and with the third vertex at the South
Pole. We denote its interior angle a t the initial position of
this point by y ; then
Using the law of cosines for a spherical triangle, we can
now write
sin pr+l=cos a sin p7-sin a cos pr cos y
and
coslAhI=(cos a-sin sin vr)/(cos CP++~ cos pT).
With the addition of the obvious relation
the displayed formulas define, in general, the computation
of the r.+l time level coordinates of the considered
point. When ur+*=O, the considered triangle collapses
into a line, but the first eq (27) is still valid; the second
eq (27) is not needed in that case. Finally, another special
case of a trajectory point coming exactly to one of the
Poles, or to its extreme vicinity, might be harmful; this
was avoided by moving the point slightly away from the
Pole should such a case tend to occur.
TRANSFORMATION OF THE VELOCITY
AND ACCELERATION COMPONENTS
When a computation point is moved to its new ~+1
time level position, its velocity components, the extra-
polated acceleration contributions already added, have
to be transformed to the new orientation of the geographi-
cal A, Q coordinate lines. Since a one-step displacement
is relatively smdl, we assume no change in the orientation
of the A, 0 coordinate lines. (This is equivalent to forward
computation of the small tan e terms in eq (1). One could,
of course, compute their space-centered values if desired.)
We can then compute the new u, v velocity components
by using eq (2), (3), and
u=u cos 8-vsin 8,
and (28)
v=u sin p+v cos 8.
Parallel to this, transformation of the nonextrapolated
acceleration components also has to be performed so that
they can be used for time extrapolation of the acceleration
components a t the T + 1 time level position of the consid-
ered point. This transformation is, of course, done in the
same way as the velocity components transformation,
except for the acceleration components replacing the
corresponding velocity components in eq (28).
5. INITIAL CONDITIONS
A number of test computations were done with the
January 1971 Fedor C
initial conditions similar to those used by Phillips (1959)
and Krishnamurti (1962). They are given by the stream
function
(29) $=-a: sin p(k,-k2 cosR cp cos Rh)
where k, and k, are constants and R (integer) is the wave
number. As shown by Haurmitz (1940) for the linearized
case and by Neamtan (1946) for the nonlinear one, the
flow pattern given by eq (29) will in a nondivergent
barotropic atmosphere move from west to east without
change of shape and with the angular velocity
R( 3 + R) k, - 2 Q '= (l+R)(2+R)
'
Having a primitive equations model, we need as initial
conditions the initial velocity components and height
of the free surface. The velocity components equivalent
to eq (29) are obtained by
The height field associated with eq (29) can be determined
by solving the balance equation; as given by Phillips
(1959), one obtains
where
2(Q+k1)kz c "[(R ~+~R +~)-((R +~)~c ~], B(cp)'(R+ 1) (R+2)
C(cp)= ~~:c "[(R +l )~~- (R+2)],
and
c=cos (p.
It has been pointed out by Phillips that eq (29-32)
will also be a solution of the corresponding divergent
primitive equations system only if the constant kl happens
to be such that the wave is stationary. This will occur
when
ki=2Q/[R(3+R) 3. (33)
If this is not so, the solution of the primitive equations
system is unknown. Thus, it appears that the best test of
the performance of a barotropic numerical primitive
equations model can be made by choosing eq (33) for the
definition of kl as originally done by Charney (1955).
Then we know that any departures of the numerical
Jlesinger 23
I
FIGURE 4.-Distribution of the height of the free surface at the
beginning (upper map) and a t the end (lower map) of the &day
forecast, shown on a cylindrical meridionally equidistant projec-
tion. The maps are printed with a shading interval of 160 gpm.
solution from the initial fields of eq (29), (31), and (32)
represent numerical errors. Hence, kl was here defined by
eq (33); the constants R and do were given the same
values as those by Phillips and Krishnamurti, that is,
R=4 and &,=9.8X8000 m2 s -~,
while the constant k, was defined at about half of their
value, that is,
to obtain velocities more comparable to those observed
in the atmosphere. The initial distribution of 4, obtained
with these values of the constants, is shown in the upper
part of figure 4. The initial distributions of u and v are
similar to those shown in the paper by Krishnamurti,
with u values approximately equal to half of those dis-
played in his figure 7 and values, since they do not
depend on kl, equal to almost exactly half of his values
in figure 8.
24 MONTHLY WEATHER REVIEW VOl. 99, No. 1
6. RESULTS OF THE EXPERIMENTS
FIGURE 5.-Space distribution of the computation points at the
beginning (upper map) and at the end (lower map) of the 4-day
forecast, shown on a cylindrical equal-area projection. The pro-
jection is conformal at latitudes of about 45'.
The present method, as defined in sections 2 to 4,
involves no special choice of the initial space distribution
of computation points, as long as each computation point
has a sufficient initial number of neighbors. Thus, an
initial space distribution of computation points also has
to be prescribed. This was done by choosing a distribution
between the random and the hexagonal one. Such quasi-
hexagonal initial distribution of computation points was
accomplished by placing the points a t constant increments
in longitude and at random on every second of 2L equal
intervals of the sine of latitude, going repeatedly K times
from the South Pole to the North Pole in a chessboard
fashion (placing the points at successive sweeps on odd
and on even of the ZL intervals). Computation points
then form an irregular K X L point grid on an equal-area
cylindrical projection of the globe with consecutive some-
what slanted columns of points shifted relative to one
another so as to form a quasi-hexagonal pattern. A par-
ticular space distribution of 160x25 points obtained in
this way is shown as the upper part of figure 5; i t was
used for the 4000-point experiments presented here.
A 4-day forecast was computed using the described
method and initial conditions. It was made with 4000
computation points, time steps of 12 niin, and time
extrapolation parameter a equal to the value of 0.25
(1 + &), which accomplishes the maximum stability
range of At. Values of the remaining numerical constants
have aIready been given in the preceding sections.
Besides this 4-day experiment, two shorter 1-day runs
have been performed: one differing from the 4-day run
only in having a coarser space resolution of 3000 points
arranged initially in a 120x25 point grid and the other
only in having the parameter a equal to X, that is, by
being computed using the Adams-Bashforth time-differ-
encing scheme.
Results of the 4-day experiment will be illustrated
here by a number of maps and diagrams. As first of these,
figure 4 shows in its lower part the geopotential field at
the end of the 4-day period; it can be compared with the
initial field shown in the upper map, both on a cylindrical
meridionally equidistant projection. These maps were
obtained by performing a space interpolation of the
geopotential values to points of the spherical 72x46
point grid (5' by 4" in longitude and latitude) and then
using a usual mapping program with a 160 gpm contour
interval. The space interpolation was done by finding
always three nearest computation points to a particular
spherical grid point and then fitting a linear polynomial
to the geopotential values at these three points-with an
additional constraint that the interpolated value must
remain within the range between the lowest and the
highest of the considered geopotential values. This con-
straint had the purpose of protecting against an unfavor-
able space distribution of the three nearest points when
they are organized nearly along a line; otherwise, a
meaningless result could be obtained in such a case. A
slight amount of noise can be introduced occasionally
into the maps by such interpolation as can be seen by
close inspection of the upper map showing the initial
field that contained no noise in the computation points
data. An alternative interpolation procedure would be to
perform a least-squares fitting to the values at the four
nearest computation points; in that way, a smoothing
effect would be obtained.
An outstanding feature of the 4-day map is the fact that
no phase errors can be noticed: locations of the trough and
ridge lines appear to be exactly the same as on the initial
map. The amount of noise present, despite all precautions,
may seem somewhat discouraging. However, most of it
appears in polar regions where the geopotential field was
initially rather flat, with the lowest 8000 gpm isolines
running along the polar lines of the projection. Thus, only
a small occasional decrease in the geopotential heights was
sufficient to produce noise in the map shading in these
regions; moreover, this noise is shown on the map with
exaggerated space dimensions due to the large area magni-
January 1971 Fedor Mesinger 25
05‘ I 2 3 4
TIME (DAYS)
FIGURE 6.--Mean value of the distances from each computation
point to its nearest neighbor as a fuiiction of time for the 4-dny
experiment.
fication in the vicinity of a Pole. In other n-ords, the un-
favorable visual irnl>rcssioi~ as to the amount of noise
present in the final gcopotential field is to some extcut
deceiving. A fair amount of noise can be seen along the
Equator \\-here thc field is also flat; but the remaining
contour lines in regions of stronger gradients are often
remarkably smooth.
In figure 5, the initial (npper map) and final (lower map)
space distribution of computation points is shown on an
cyual-area cyliiiclrical projection so as to reflect their
number per unit area of the globc. One can notice that, in
the 4-day map, the points exhibit a tendcncy to organize
dong irregular lines; this must be a consequence of the
fact that, except for somc iioisc, 110 turbulence \vas pre-
sent in the computed flo\v fields. A more irregular final
distribution of points should of course be expected when
the points arc advccted by a less orgunizcd flow.
Starting with an initial distribution of points that is
better than random, onc must expect i t to deteriorate
with time and eventually become random or possibly 11-orse
than random--“bctter” and “v-orsc” meaning that the
points are more and less c\-enly dispersccl than in a random
distribution , rcspcctivcly. In fact, somc earlier cxperimcnts
(Mesinger 1965) have sho\\-ii that the iiondivcrgcnt part
of an observed t\\-o-dimciisioiial atmospheric flo\\- had thc
property of producing as well as maintaining a raiidom
distribution of initially organized and floating particles
while divergent flows maintained distributions worse than
random in the above sciisc. The time rate of this process
is of‘ interest and d l be discussed here along with its
effect on the accuracy of the forecast.
For giving a statistical description of the distribution
of our computation points, it is convenient to consider
the mean distance from each point to its nearest neighbor.
It has been shown by Hertz (1909) that, when a very
large number of points is distributed at random on a plane
and the unit of length chosen so as to make their average
concentration equal to 1, this mean distance will be equal
to 0.5. If, however, the points are distributed at random
only relative to one another but with a concentration
that is a function of space, this mean distance has to be
less than the random distribution value of 0.5 (Mesinger
fn zs n
I 2 3 4
TIME (DAYS1
FIGURE 7.-Root-mean-squcwe individual change per time step in
the geopotential height as a function of time. The values refer
to the 4-dag- experiment.
1965). Finally, if the points are organized relative to one
another, this mean distance can of course be greater;
it is, for example, equal to &/4.&=1.0746 when the points
are arranged into a hexagonal grid. Apparently when
the mean nearest neighbor distance is greater, me
might expect a better performance of a particular distri-
bution of computation points. This mean distance was
computed at 1-hr intervals in the 4-day experiment: its
initial value, corresponding to the distribution shoxn in
the upper map of figure 5, was equal to about 0.7555;
i t decreased rather slonrly later on, reaching a value of
about 0.6319 at the elid of the experiment. The obtained
mean distances as a function of time are shown in figure 6.
It can be compared with figure 7 in the aforementioned
paper by the author; in that previous experiment, the
iiondivergent part of an observed flow needed much less
time, only about half a day, to accomplish the same de-
crease in the mean nearest neighbor distance. The dif-
ference, again, is due to the lack of turbulence in the
present experiment.
For obtaining a more quantitative idea of the amount of
various kinds of noise present in the 4-day experiment,
a record of the root-mean-square (rms) individual change
per time step in the geopotential heights is shown in
figure ’i-also based on computations at 1-hr intervals.
Computed per unit time and divided by the rms value of
the geopotential height, this quantity would give the
mass-weighted rms value of div (divergence) v. Thus,
i t is equal to zero in the analytic solution and does not
depend on time.
I n the experiment, the initial value of the rms individual
geopotential change was equal to about 3.38 m2 s-’/12
min; this corresponds to the value of about 0.5X10-’ s-l
for the mentioned value of div v. Compared with the
Pole-to-Equator change in the free-surface height of
somewhat more than 1600 gpm, this initial error of about
% gpm/12 min appears to be moderate; in fact, some
caution should be exercised here since these numbers
may already be influenced by the round-off error of the
performed IBM single precision computations. Since
there is no noise in the initial fields, this error-assuming
no round-off influence-reflects only the difference be-
26 MONTHLY WEATHER REVIEW Vol. 99, No. 1
zi
I 2 3 4
TIME IDAYS)
F~GURE 8.-Average height of the free surface as a function of
time for the &day experiment.
tween the analytic space derivatives of velocity com-
ponents and those computed by the least-squares fitting of
a second-degree polynomial to the values a t available
neighbors. Thus, it shows a combined effect of the incon-
venient space distribution and number of neighbors and
of the truncation of higher order terms in the polynomial.
Some feeling for the relative significance of these two
factors can of course be obtained through suitable experi-
ments. For instance, the 3000-point run, with points
arranged into a 120x25 pattern, had this initial error
equal to 3.97 m2 r2/12 min; this must have been mostly
a consequence of the increased importance of the neglected
higher order terms. On the other hand, a different organi-
zation of the 4000 points into a 100x40 instead of 160 X25
pattern with a 0.7997 mean nearest neighbor distance
was found to result in a reduction of this error to the value
of 2.75 of the same units.
During the first day or so, the considered quantity is
seen to undergo vigorous gravity-type oscillations with
a period of very nearly 3 hr. These oscillations must
be set off by small systematic differences between the
analytic and computed space derivatives. By the end
of about 1 day, the oscillations are damped out through the
process of time differencing. The parallel run, differing
from the 4-day run only by being computed with the
Adams-Bashforth time-diff erencing scheme, did not ex-
hibit such damping; and the oscillations persisted without
noticeable change in amplitude throughout the l-day
period of the experiment.
At times after about 1 day, the rms individual geo-
potential change seen in figure 7 settles at the value of
about twice the initial one, showing the contribution of
the random space noise in the velocity fields. After about
2 days, a slow rise in this rms change can be noticed,
interrupted by three irregular sudden increases. The
gradual rise in this quantity should probably be attributed
to the slow deterioration of the space dispersion of compu-
tation points, and the three peaks to occasional occur-
rences of very inconvenient configurations of neighbors.
Performance of the model can further be analyzed by
examination of the behavior of global integrals of a number
of relevant physical quantities. Conservation of none of
these integrals is formally guaranteed in the present
method, and a change in conservative quantities thus
reflects the effect of numerical errors or d the artificial
smoothing mechanisms. A record of one of such quantities,
the average height of the free surface, is shown in figure 8,
again based on values computed a t l-hr intervals. This
average was obtained by using the geopotential values
interpolated to points of the spherical 72x46 point grid,
by interpolation described earlier in this section. The
interpolation protects the computed average values
against the effect of possible systematic variations in the
space density of computation points; this protection
seems necessary since such variations would be correlated
with those in the height of the free surface. The time
changes in the average free-surface height shown in
figure 8 appear to be random and of a modest amplitude.
The observed amplitude of less than 2.1 gpm represents
about 0.02 percent of the total average height. The
considered average is of course proportional to the total
mass of the fluid.
It may be noted that, using eq (32), one can obtain
(34)
where
is the global average value of 4. With the values of
constants given in section 5, eq (34) gives 9057.33 gpm
for the average free-surface height. The initial value of
9057.02 gpm entered in figure 8 happens to be slightly
less than this "exact" value, the difference showing the
amount of the sampling error. This exact value could
also have been entered t o start the curve in figure 8,
what would have further reduced its amplitude.
Finally, a record of the average values of total, kinetic,
and available potential energy computed also at 1-hr
intervals is shown in figure 9. Denoting the total kinetic
energy averaged per unit mass by K, and the available
potential energy averaged in the same way by A,, we have
and
(35)
(36)
where S represents the total area of the globe. I n computing
the values of K, and A,, shown in figure 9, these intergals
and global average of 9 in eq (36) were substituted by
corresponding sums and by the average of based on
values of 4 and v2 interpolated again to points of the
spherical 72x46 point grid. Thus, the values obtained
for Kl and A, were also protected against systematic
January 1971 Fedor Mesinger 21
TOTAL
650
,n
r AVAILAELE POTENTIAL
1 2 3 4
TIME (DAYS1
FIGURE 9.-Average values per unit mass of the total, kinetic, and
available potential energy as functions of time. The solid curves
refer to the Cday experiment; the dashed curves refer to the
parallel 3000-point experiment, differing from the +day one only
in the number of computation points; and the dotted curves refer
to the parallel a= 54 (Adams-Bashforth scheme) experiment, dif-
fering from the &day one only in the value of this parameter
of the timedifferencing scheme.
variations in the space density of computation points.
Such a protection seems desirable since it appears that
there should exist a definite correlation between these
variations and those in the contributions to total energy
values. However, it is this space interpolation that gives
rise to the time noise seen in the kinetic (and thus also
total) energy curves in figure 9; no time noise appears
when the integrals in eq (35) are computed by taking
simple sums of corresponding quantities over all the corn-
putation points. Of the energy curves, the solid ones
obviously refer to the 4-day experiment; the dashed curves
refer to the “coarse resolution’’ 3000-point run, and the
dotted ones to the run performed by using the nondamp-
ing Ad ams-Bashf orth time-diff erencing scheme.
An outstanding characteristic of the energy curves is
a periodic exchange between the kinetic and available
potential energy. These oscillations are apparently brought
forth by the truncation of the higher order terms in the
least-squares polynomial. Namely, it would seem that in
this way one computes pressure gradients that are pre-
dominantly smaller, in their absolute values, than the
analytic ones. Since the initial Coriolis forces are not
affected by the space-differencing method, an imbalance is
obtained such as to set off inertia-gravity oscillations in
which the fluid particles initially move predominantly
toward higher free-surface heights; this is associated
with an increase in available potential and a decrease in
kinetic energy. The amplitude of these oscillations should
then increase with an increase in the neighbor-defining
distance due to the increased importance of the neglected
higher order terms; this is confirmed by the dashed
kinetic energy curve. Namely, since the average number
of neighbors E was the same in all of the three experi-
ments, the neighbor-defining distance of the 3000-point
experiment was greater than that of the two 4000-point
ones; accordingly, this 3000-point experiment exhibits
a greater amplitude of the discussed oxcillations. I n the
4000-point 4-day experiment, this amplitude initially
amounts to about 1.1 percent of the total kinetic energy
and decreases later on. Damping in the time differencing
and physical dispersion of the waves are responsible for
this decrease in amplitude. The period of the oscillations
is less than the pure inertial one, in accordance with the
usual linear theory; see, for example, the paper by Arakawa
(1970).
The energy curves in figure 9 exhibit a general de-
creasing trend, showing the loss of energy due to the
artificial space smoothing and damping in the time
differencing. During the 4-day experiment, the total
kinetic energy decreases by about 3.71 percent of its
initial value. One part of this energy is transformed into
the available potential energy through the described
inertia-gravity oscillations and subsequent geostrophic
adjustment; it somewhat more than compensates for the
loss in the latter to the existing space-time smoothing
mechanisms, leaving the potential energy a t a value
slightly higher than the initial one. The remaining part
is dissipated through the numerical smoothing processes.
It is tempting to try estimating the relative contribu-
tions of space smoothing and time differencing to the
observed kinetic energy loss. One may attempt to do this
with the help of the results of the parallel run computed
using the Adams-Bashforth time-differencing scheme.
This run was terminated at the l-day time; a t that
moment, its loss in kinetic energy amounts only to some-
what less than one-third of that demonstrated by the
solid kinetic energy curve in figure 9. However, some
caution should be exercised in interpreting this difference
since there is a possibility that it appears partly as a
consequence of the weak instability of the Adams-Bash-
forth scheme. No visual sign of such instability, though,
could be noticed by careful comparison of the two 1-day
geopotential height maps; thus, a value for the parameter
a smaller than the one used for the 4-day experiment
appears preferable.
7. CONCLUDING REMARKS
The proposed Lagrangean method offers a number of
distinct advantages over the conventional Eulerian tech-
nique. We shall summarize them here:
Formulation of the problem involves no advection terms; this
makes the governing equations simple, free of nonlinear instability,
and guarantees, except for truncation errors, the conservation of
MONTHLY WEATHER REVIEW Vol. 99, No. 1
individual properties following the motion of fluid particles. In
other words, there is no computational dispersion in the advection
terms. This should improve the simulation of mesoscale phenomena,
such as fronts, jets, and intertropical convergence zones, associated
with sharp gradients of an advected quantity. Further, this enables
a realistic simulation of the deformation-dependent process of
lateral diffusion.
The simulated flow is not subjected t o systematic space changes
in grid resolution and properties. With computation points moving
with the fluid, a quasi-homogeneous resolution is accomplished
without any geometrical considerations.
The two mentioned properties significantly reduce the amount
of work needed to code a computer program for the model. Absence
of considerations related to the geometry of a particular grid and
to the treatment of advection terms relieves the programmer of a
major effort that these considerations require in a Eulerian global
model.
It would appear easy to accomplish a variable resolution, such
as t o have a higher density of computation points in regions of
rapid spatial change in appropriate quantities. Such regions exhibit
a strong tendency to move with the fluid and would thus maintain
this increased resolution for an extended time. A density dependent
neighbor-defining distance could then be used, while still observing
the principle of property of being a neighbor having to be a reciprocal
one.
The model can accept information at arbitrary points in space,
and they can be prescribed different weights, depending on their
reliability. This seems ideally suited to the expected future global
observation system (ICSU/IUGG 1967).
On the other hand, a number of objections.could be
raised :
The present experiments showed a tendency for the appearance
of computational noise, and artificial smoothing was needed t o
suppress this tendency. However, the situation should be better
when a simulation of the real atmospheric diffusion is incorporated
into the model.
It may be difficult or impossible to construct schemes that would
conserve global integrals of relevant quantities.
It seems that more computation time and more points are needed
than in a comparable Eulerian model. The 4000-point experiments
described here required about 50 min of the IBM 360/91 computer
time per 1 day of simulated time.
Finally, some uncertainty remains as to horn severe the
performed test was; it would appear that a task of main-
taining a stationary state, with computation points
floating through it, is fairly demanding. Still, one feels a
need for a more nonlinear two-dimensional test and of
course a three-dimensional one. The amount of effort
invested into the present method is exceedingly small
compared to that which brought the Eulerian methods to
the present state; this, coupled with the results of the
present study, mould seem to highly encourage further
investigation of the proposed method.
ACKNOWLEDGMENTS
The author is grateful t o Profs. Yale Mintz and Akio Arakawa
for their encouragement of this research and for some helpful
suggestions. M&ssrs. Wayne Schubert and Robert Haney kindly
helped t o prepare the manuscript by proposing numerous im-
provements in its text. Computing assistance was obtained from
the Campus Computing Network, University of California, Los
Angeles.
REFERENCES
Arakawa, Akio, “Computational Design for Long-Term Numerical
Integration of the Equations of Fluid Motion: Two-Dimensional
Incompressible Flow, Part I,” Journal of Computational Physics,
Vol. 1, No. 1, J ~l y 1966, pp. 119-143.
Arakawa, Akio, “Nunierical Simulation of Large-Scale Atmospheric
Motions,” Numerical Solution of Field Problems in Continuum
Physics, American Mathematical Society, Providence, R.I., 19710,
Arakawa, Akio, Katayama, Akira, and Mintz, Yale, “Numerical
Simulation of the General Circulation of the Atmospherc,’’
Proceedings of the WMOlIUGG Symposium on Numerical Weather
Prediction, Tokyo, November 26-December 4, 1968, Japan Mcte-
orological Agency, Tokyo, Mar. 1969, pp. IV-7-IV-12.
Charney, Jule G., “The Use of the Primitive Equations of Motion
in Numerical Prediction,” Tellus, Vol. 7, No. 1, Feb. 1955, pp.
Charney, Jule G., “Some Remaining Problems in Numerical Weather
Prediction,” Advances in Numerical Weather Prediction, Travelers
Research Center, Inc., Hartford, Conn., 1966, pp. 61-70.
DjuriC, DuBan, “On the Accuracy of Air Trajectory Computations,”
Journal of Meteorology, Vol. 18, No. .5, Oct. 1961, pp. 597-605.
DjuriC, DuSan, “The Role of Deformation in the Large-Scale Dis-
persion of Clusters,” Quarterly Journal of the Royal Meteorological
Society, Vol. 92, No. 392, Apr. 1966, pp. 231-238.
Fj$rtoft, Ragnar, “On a Numerical Method of Integrating the Baro-
tropic Vorticity Equation,” Tellus, Vol. 4, No. 3, Aug. 1952,
FjGrtoft, Ragnar, “On the Use of Space-Smoothing in Physical
Weather Forecasting,” Tellus, Vol. 7, No. 4, Nov. 1955, pp. 462-
480.
Grammeltvedt, Arne, “A Survey of Finite-Diff erence Schemes for
the Primitive Equations for a Barotropic Fluid,” Monthly Weather
Review, Vol. 97, No. 5, May 1969, pp. 384-404.
Haurwitz, Bernhard, “The Motion of Atmospheric Disturbances on
the Spherical Earth,” Journal of Marine Research, Vol. 3, 1940,
Hertz, Paul, “Uber den gegenseitigen durchschnittlichen Abstand
von Punkten, die mit bekannter mittlerer Dichte im Raume
angeordnet sind” (On the Average Mutual Distance of Points
Distributed in Space With a Known Average Density), Mathe-
matische Annalen, Vol. 67, Springer-Verlag, Berlin, 1909, pp. 387-
398.
ICSU/IUGG-Committee on Atmospheric Sciences, Report of the
Study Conference Held at Stockholm, 28 June-11 July 1967, on the
Global Atmospheric Research Programme (G A R P ), World Mete-
orological Organization, Geneva, 1967, 144 pp.
Krishnamurti, Tiruvalam N., “Numerical Integration of Primitive
Equations by a Quasi-Lagrangian Advective Scheme,” Journal
of Applied Meteorology, Vol. 1, No. 4, Dec. 1962, pp. 508-521.
Kurihara, Yoshio, and Holloway, J. Leith, Jr., “Numerical Integra-
tion of a Nine-Level Global Primitive Equations Model Formu-
lated by the Box Method,” Monthly Weather Review, vel. 95,
Leith, Cecil E., “Lagrangian Advection in an Atmospheric Model,”
World Meteorological Organization Technical Note No. 66, Geneva,
Lilly, Douglas K., “On the Computational Stability of Numerical
Solutions of Time-Dependent Non-Linear Geophysical Fluid
Dynamics Problems,” Monthly Weather Review, Vol. 93, NO. 1,
Jan. 1965, pp. 11-26.
Matsuno, Taroh, “Numerical Integrations of the Primitive Equa-
tions by a Simulated Backward Difference Method,” Journal of
the Meteorological Society of Japan, Ser. 11, Vol. 44, NO. 1, Feb.
1966a, pp. 76-84.
pp. 24-40.
22-26.
pp. 179-194.
pp. 254-267.
NO. 8, Aug. 1967, pp. 509-530.
1965, pp. 168-176.
January 1971 Fedor Mesinger 29
Matsnno, Taroh, “False Reflection of Waves at the Boundary Duc
to the Use of Finite Differcnces,” Journal of the Meteorological
Society of Japan, Ser. 11, Vol. 44, No. 2, Apr. 1966b, pp. 143-157.
Mesinger, Fedor, “Behavior of a Very Large Number of Constant-
Volume Trajectories,” Journal of the Atmospheric Sciences, Vol. 22,
No. 5 , Sept. 1965, pp. 479-492.
Minte, Yale, “Very Long-Term Global Integration of the Primitive
Equations of Atmospheric Motion,” World Meteorological Organi-
zation Technical Note No. 66, Geneva, 1965, pp. 141-167.
Neamtan, S. M., “The Motion of Harmonic Waves in the Atmos-
phere,” Journal of Meteorology, Vol. 3, No. 2, June 1946, pp. 53-56.
Nitta, Takashi, “Initialization and Analysis for the Primitive
Equation Model,” Proceedings of the WMOII UGG Symposium on
Numerical Weather Prediction, Tokyo, November 26-December 4,
1968, Japan Meteorological Agency, Tokyo, Mar. 1969, pp.
VI-1 1-VI-20.
qkland, Hans, “An Experiment in Numerical Integration of the
Barotropic Equation by a Quasi-Lagrangian Method,” Geo-
jysiske Publikasjoner, Vol. 22, No. 5, Oslo, Jan. 1962, 10 pp.
Phillips, Norman A., “Numerical Integration of the Primitive
Equations on the Hemisphere,” Monthly Weather Review, Vol. 87,
Phillips, Norman A., “The Equations of Motion for a Shallow
Rotating Atmosphere and the Traditional Approximation,”
Journal of the Atmospheric Sciences, Vol. 23, No. 5, Scpt. 1966,
Sadourny, Robert, Arakawa, Akio, and Minte, Yale, “Integration
of the Nondivergent Barotropic Vorticity Equation With an
Icosahedral-Hexagonal Grid for the Sphere,” Monthly Weather
Review, Vol. 96, No. 6, June 1968, 351-356.
Welander, Pierre, “Studies on the General Development of Motion
in a Two-Dimensional, Ideal Fluid,” Tellus, Vol. 7, No. 2,
May 1955, pp. 141-156.
Wiin-Nielsen, Aksel, “On the Application of Trajectory Methods in
Numerical Forecasting,” Tellus, Vol. 11, No. 2, May 1959,
Wurtele, Morton G., “On the Problem of Truncation Error,’’ Tellus,
NO. 9, Scpt. 1959, pp. 333-345.
pp. 626-628.
pp. 180-196.
Vol. 13, NO. 3, Aug. 1961, pp. 379-391.
[Received March 5, 1970; revised May 6, 19701
I
CORRECTION NOTICE
Vol. 98, No. 7, July 1970: p. 531, photo of fig. 3 should be interchanged with
photo of fig. 4.
I