NASA Technical Memorandum 104606, Vol. 2

Technical Report Series on

Global Modeling and Data Assimilation

Max J. Suarez, Editor

Goddard Space Flight Center

Greenbelt, Maryland

Volume 2

Direct Solution of the Implicit Formulation of

Fourth Order Horizontal Diffusion for

Gridpoint Models on the Sphere

Yong Li

General Sciences Corporation

Laurel, Maryland

S. Moorthi

National Meteorological Center

Camp Springs, Maryland

J. Ray Bates

Goddard Space Flight Center

Greenbelt, Maryland

National Aeronautics and

Space Administration

Goddard Space Flight Center

Greenbelt, Maryland

1994

Abstract

High order horizontal diffusion of the form KV 2m is widely used in spectral models

as a means of preventing energy accumulation at the shortest resolved scales. In the

spectral context, an implicit formulation of such diffusion is trivial to implement. The

present note describes an efficient method of implementing implicit high order diffusion

in global finite difference models. The method expresses the high order diffusion equation as a sequence of equations involving V _. The solution is obtained by combining

fast Fourier transforms in longitude with a finite difference solver for the second order

ordinary differential equation in latitude.

The implicit diffusion routine is suitable for use in any finite difference global model

that uses a regular latitude/longitude grid. The absence of a restriction on the timestep

makes it particularly suitable for use in semi-Lagrangian models. The scale selectivity

of the high order diffusion gives it an advantage over the uncentering method that has

been used to control computational noise in two-time-level semi-Lagrangian models.

°°°

111

PAGE IILANK NOT FILM_"I_

1'

Contents

1 Introduction

2 Description of The Solver

2.1 Equations .....................................

2.2 Longitudinal and Latitudinal Discretization ..................

2.3 Boundary Conditions at The Poles .......................

2.4 Numerical Solution ................................

2.5 Coding Strategy .................................

2.6 Extension to Higher Orders ...........................

3 Justification for Using Fourth Order Diffusion

4 Application to A Semi-Lagrangian Shallow Water Model

4.1 Initial Condition of McDonald and Bates (1989) ...............

4.2 Real Data Case ..................................

5 Application to a 3-D Semi-Lagrangian GCM

6 Concluding Remarks

APPENDIX A: Fortran Program of The Solver

References

FIGURES

2

2

3

4

5

6

7

7

9

9

10

10

11

13

28

3O

v
List of Figures

4

5

The damping rates as a function of wavenumber for the second and the fourth

order diffusion in both implicit and explicit formulations ...........

(a) Spherical harmonic wave number one; (b) Same as (a), but with 2Ax

waves embedded; (c) After applying the semi-implicit formulation of fourth

order diffusion to (b); (d) After applying the implicit formulation to (b); (e)

After applying both the semi-implicit and the implicit formulations to (b)..

(a) Geopotential at 50 days for the integration with the initial condition of

McDonald and Bates (1989) and a value of 0.03 for the uncentering parameter(contours every 60.0 meters); (b) Same as (a) but with the fourth order

diffusion applied to geopotential(contours every 60.0 meters); (c) Same as

(a) but with the fourth order diffusion applied to divergence (contours every

60.0 meters); (d) Same as (a) but with the fourth order diffusion applied to

both geopotential and divergence(contours every 60.0 meters) .........

(a) Geopotential of the ECMWF analysis of 0000 UTC, 1 December, 1992

(contours every 60.0 meters); (b) Geopotential after 3 days of integration with

a value of 0.05 for the uncentering parameter(contours every 60.0 meters);

(c) Same as (b) but with a value of 0.015 for the uncentering parameter and

the the fourth order diffusion applied to both geopotential and divergence

(contours every 60.0 meters) ...........................

50 mb temperature field from NASA/GLA 3-D semi-Lagrangian semi-implicit

GCM without the application of the horizontal diffusion. The value for the

off-centering parameter is 0.05. (a) 3-day forecast; (b) 5-day forecast; (c)

7-day forecast; (d) 10-day forecast. (Contour interval is 3 degrees) ......

50 mb temperature field from NASA/GLA 3-D semi-Lagrangian semi-implicit

GCM with the application of the horizontal diffusion. The value for the offcentering parameter is 0.05. (a) 3-day forecast; (b) 5-day forecast; (c) 7-day

forecast; (d) 10-day forecast. (Contour interval is 3 degrees) .........

3O

31

34

36

38

4O

vi
1 Introduction

Diffusion and filtering have been common features of numerical models since the beginning

of numerical weather prediction (Shuman, 1957). The reason for this is two-fold. First,

diffusion and filtering can serve as crude parameterizations of subgrid scale processes that

are not represented by the model. In the real atmospheric flow, there is a cascade of energy

to smaller and smaller scales due to nonlinear interactions, until the motion is dissipated

by friction. But in a model, even with the largest computer available, only scales far larger

than the inertial subrange of turbulence can be resolved. In other words, what needs to

be parameterized is much more than molecular diffusion. Secondly, models are usually

subject to undesirable high wavenumber noise due, for example, to numerical discretization

procedures. Diffusion and filtering can serve to reduce or eliminate this kind of noise.

Explicit formulations of diffusion equations are only conditionally stable (Mesinger and

Arakawa, 1976). As a consequence, the coefficient of diffusion cannot be increased beyond a

critical value that depends on gridlength and timestep. This imposes a stringent restriction

in application to atmospheric models, especially those such as semi-Lagrangian models that

allow the use of large timesteps: the allowed diffusion coefficient is too small to eliminate

the high wavenumber noise. Such a restriction is particularly marked in global models using

a uniform latitude/longitude grid where the convergence of the meridians leads to a grid

size much smaller at high than at low latitudes. A diffusion coefficient that is barely large

enough to ensure smooth solutions at low latitudes can far exceed the critical value for

computational stability at high latitudes. This restriction becomes even more marked if a

diffusion equation of higher order is used.

An alternative to diffusion are spatial filters; see Raymond and Garder (1991) for a review.

The most commonly used in meteorological models is the Shapiro filter (Shapiro, 1971),

which has the advantages of being very simple in formulation and of being easily extendible

to higher orders. The filter is designed in a manner analogous to that used in solving

the explicit formulation of a diffusion equation of equivalent order with the computational

stability criterion satisfied. When large timesteps are used, the implied diffusion coefficient

is so small that high wavenumber noise cannot be effectively reduced. In addition, the

preferred spherical formulation of the filter, which involves simply averaging among the

surrounding points, does not preserve area-weighted mean values. This implies violation of

global conservation of the quantity being filtered. All attempts at using the filter in the

context of spherical geometry have been shown to lead to degraded performance at high

latitudes (Shapiro, 1979).

In some semi-Lagrangian models, an uncentered discretization has been employed to damp

computational noise (Tanguay, et al., 1992; Bates, et al., 1993; Gravel, et al., 1993). This,

however, can reduce the accuracy of those terms discretized semi-implicitly (the second

order accuracy is decreased to first order). Also, being insufficiently scale selective, the

uncentering can have a negative impact on the large scale features in both amplitude and
phase,ashasbeenrecognizedin our experience and also by Semazzi and Dekker (1994).

Implicit high order horizontal diffusion is trivial to implement in spectral models since the

model variables are expressed in spectral harmonics which are eigenfunctions of the spherical

Laplace operator. In this context, implicit diffusion is widely used (e.g., Williamson and Olson, 1994; Ritchie, et al., 1994). The need for an implicit formulation of horizontal diffusion

in gridpoint models has also been recognized, but the solution is much less straightforward.

Jakimow et al. (1992) presented an implicit formulation for second order horizontal diffusion, and showed that its application to a semi-Lagrangian hemispheric model gave positive

results. It was found that the optimum value for the diffusion coefficient exceeds the value

allowed by the stability criterion of an explicit formulation. McDonald (1994) has developed an implicit solver for an approximated form of V 2m diffusion for limited area finite

difference models, and has found the diffusion to be effective in reducing noise in the semiLagrangian context. Our experience with the NASA/GLA Semi-Lagrangian GCM (Bates

et al., 1993; Moorthi et al., 1994) indicates that at moderate spatial resolution the intrinsic

smoothing associated with the interpolations, along with a small uncentering, can prevent

the growth of high wavenumber noise. But in high resolution integrations (e.g., 1°× 1.25 ° in

the horizontal with 46 vertical layers), computational noise is present even with increased

uncentering, especially at high latitudes and near orography (Moorthi, 1994). To alleviate

this problem, the present implicit formulation of fourth order horizontal diffusion for a regular latitude/longitude grid on the sphere has been developed. Some justification for the

use of fourth order rather than second order diffusion is also presented.

A detailed description of the solver is presented in Section 2. The justification for using fourth order diffusion is discussed in Section 3. Results of applying the solver to a

global shallow water semi-Lagrangian model are given in Section 4. Concluding remarks

are presented in Section 5. The Fortran program of the solver is attached to the note as an

Appendix.

2 Description of The Solver

2.1 Equations

The fourth order horizontal diffusion equation takes the form

and can be discretized as

cn+l __ Ca

At _ _KV4 ((1- _)¢n+l +_¢n)

(1)

(2)

2
where_ is a parameterwhosevalueis determinedby the type of temporaldiscretization

used;it is semi-implicitif _ is 1⁄2and implicit if _ is 0. (This terminologyis adoptedfor

consistencywith themeteorologicalliterature;thetermsCrank-Nicholson and Euler implicit

are corresponding descriptions used elsewhere.)

Integrating (2) involves the solution of the following equation, in a manner similar to that

used by Lindzen and Kuo (1969):

i. e.,

V4¢"+1 -]- KAt (1 - _) - KAt (1 - _) ¢" (1 - _) (3)

V4¢-t- re = f (4)

where r is a prescribed parameter, f is the forcing term and the superscript (n + 1) has been

omitted. Numerical solution of (4) on the sphere requires a conservative discretization of

the V 4 operator, with a consistent treatment of the polar boundary conditions. Difficulties

associated with these as well as with the cross derivative terms can be avoided by writing

the equation in the form of two second order equations:

where

v2z + r¢ = f (5)

v2¢- z = 0 (6)

V2 - 1 02 1 oq /cosg0__

-- a 2cos 290A 2 + a 2cos/_08 \ t_0]

2.2 Longitudinal and Latitudinal Discretization

Since the coefficients in the above equations are independent of longitude, we can perform

Fourier transforms using the Fast Fourier Transform (FFT) algorithm in the A-direction:

= (i) (7)

Then the partial differential equations (5) and (6) lead to the ordinary differential equations

(ODE) in vector form:

L)( + Q)( = F (8)

where

L- a2cosSd8 cos0

3
k 2

q -- a 2 cos 2 0"

The discretized forms of L and q at the interior points are:

(Ly)j = ujyj+l + vjyj + wjyj-l

sin2 (_)

qJ= 2

( _ ) a2 cos2 0j

Here, y is any variable operated on by L and

, cos Oj_ 1⁄2 cos 0j+ _ cos 0j+ _ + cos 0j_ _ w_ = (9)

uj = a 2 cos0j (A0) 2 vj = a2c°sSJ (A0) 2 a2c°sSJ (A0) 2

for j = 1,2,3,...,N - 1 (j = 0 and j = N represent the south and the north poles,

respectively). We denote

Aj = wjI Bj = vii + Qj Cj = ujl (10)

where I is the unit diagonal matrix. Thus, the discretized equation for the interior points,

j = 1,2,3,...,N - 1, can be written as:

Aj._j-1 + Bj._j + Cj._j+I = -Pj. (11)

2.3 Boundary Conditions at The Poles

The boundary conditions at the south and the north poles are defined as:

A0)¢0 + B0)(, = F0 (12)

AN.XN-I + BN.XN = fiN. (13)

At each pole, there should be a single unique value for the solution Z and ¢. For all

longitudinal asymmetric components (k > 0), the boundary conditions for the Fourier

components )( at the poles simply reduce to:

The Fourier component (k = 0) is the zonal mean, whose solution is derived from the ODE

for it. The boundary conditions at the poles can be determined simply by discretizing (5)

and (6) at the poles. We arrive at

4
_(_): _ _(0) + 0 (_: _ _(:);

for the south pole, and similarly

- : o (_):

(f (0N)) (15)

for the north pole. The overbar denotes zonal average. These correspond to the boundary

conditions of the form in (12) and (13) with the following coefficient matrices:

Bo

_o

AN

((_°): o)

= _ ; (:7)

o (_o)

((_°)_ o)

= _ ; (19)

o (_)

BN :- 2 ); (20) -(_o)

fiN: (](oN)) • (21)

2.4 Numerical Solution

We again adopt the method of Lindzen and Kuo (1969) to solve the system (ll) coupled

with (12) and (13). Assume that the solution of (11) takes the following form:

£j = _j£j+l + _j (22)
where c_j is a 2 × 2 matrix and/3j is a vector. We know from (12) that

ao = - Ao-l Bo

_o --Ao-1*oo.

If we substitute

(23)

(24)

Xj-1 = OLj-IXj +/_j-1

into (11), we have

x, = -(Ajwj_, + Bj)-' CjXj+I -[- (AjoLj-1 -_ Bj) -1 (Fj - Ajt_j-1). (25)

Comparing (25) and (22) leads to:

oej = - (Ajoej_, -'b Bj) -1 Cj (26)

/_j: (Aj_j-I -I-Bj)-I (*Oj- Aj_j-x) (27)

for j : 1, 2, 3,..., N - 1. With the following two equations:

AN.)(N-1 "q-BN._N-- *ON (28)

.eYN- 1 = OtN-l £ N "_ I_N_I (29)

we can obtain the solution for the north pole as:

)?N= (ANOzN-I + BN) -1 (*ON -- ANflN-a). (30)

Now we recursively use (22) to obtain the complete set of Fourier components for each

j = N - 1,N - 2,...,2, 1,0. Finally, we use inverse FFT to obtain the solution at grid

points.

2.5 Coding Strategy

The coding strategy used here is similar to that of Moorthi and Higgins (1992, 1993). The

FFT has been made very efficient by using the CRAY subroutine RFFTMLT to transform

(both forward and backward) all latitudes simultaneously. In axtdition, the code for solving

the second order O.D.E's for the FFT components and the zonal mean is vectorized over the

wave numbers. After vectorization, a speed-up by a factor of 40 was observed on the CRAY

C-90 computer. The Fortran program with comprehensive documentation is supplied in the

Appendix.
2.6 Extension to Higher Orders

Extendingthesolverpresentedaboveto higherorder is straightforward, both in formulation

and in coding. Suppose we are to solve a mth order equation (m being an even number). In

the place of (2.5) and (2.6), we now have a set of m/2 second order equations. The FFT,

latitudinal discretization and numerical procedures are almost the same. The matrix size

of 2x2 then becomes m/2 x m/2. As long as m is not too large, say, not larger than 8, the

computational cost is not expected to increase dramatically.

3 Justification for Using Fourth Order Diffusion

It is well known that higher order diffusion can effectively damp the high wavenumber noise

but leave the low wavenumbers very little affected. Here we present a discussion of the

damping properties of the second and fourth order diffusion equations, in the context of

both implicit and semi-implicit formulations. Though these properties have been discussed

in earlier works (e.g., Mesinger and Arakawa, 1976), we present these here for the convenience of reference. We use the 1-D diffusion equation for illustration and then discuss the

implication for spherical applications.

The second order diffusion equation in its discretized form can be written as:

¢7+1 - ¢7

where _ = 0 corresponds to the implicit scheme and _ = 1⁄2 the semi-implicit scheme. The

notation 5_x is the finite difference approximation of the second order derivative, which is:

5x_A = Ax+l + Aj-I - 2.0AI

(Ax)2

Assume the solutiontakesthe followingform:

¢7 = ¢0 h"e/k/_'.

With second order centered difference for the spatial discretization, the damping rate, A2,

becomes:

1 - 4_r2 sin2 (_) (32)

A2= l+4(1-()r2sin2(_) "

_At 21r

where At-_¥) is wave r2 = _ and m (= the number measured by grid intervals.

The fourth order diffusion equation in its discretized form is

At (33)

7
where

_r_=_:rA= (_xA)_+l + (_r_A)/-1 - 2.0 (_=_=A)l

(Ax) 2

is the finite difference approximation of the fourth order operator. The damping rate, A4,

is

A4= 1- 16_ r4 sin4 (_) (34)

1 + 16(1 - _) r4sin 4 (_)"

KAt where r4 = ,_-:-_, •

V'x)

From (32) and (34) it is clear that for both second order and fourth order diffusion, the fully

implicit treatment (_ = 0) is unconditionally stable. When 0 < _ < 1, there is a critical

value for K above which there is a reversal of the sign of the damping rate. When _ > 1⁄2,

there is also a critical value of K above which the scheme becomes unstable, this critical

value being the smallest when the scheme is fully explicit (i.e., _ = 1). When c is nonzero,

it is possible to choose a value of K such that 2Ax waves are completely eliminated.

Figure 1 shows the damping rates as a function of wave number for the second and the fourth

order diffusion equations in both implicit and semi-implicit formulations. The parameters

r2 and r4 are so chosen that the 2Ax waves are completely eliminated with the semiimplicit formulations. Three features are observed: 1) In all cases, IAI < 1, which indicates

unconditional stability; 2) The fourth-order equation is much more scale selective than the

second order one; while the second order diffusion damps the 20Ax waves by about 5%,

the fourth order leaves it almost unchanged; 3) The semi-implicit scheme can completely

eliminate the 2Ax waves but the implicit scheme damps it only by a factor.

The above results might mislead one to believe that the semi-implicit scheme is better than

the implicit one since the former has the potential of eliminating the 2Ax waves completely

with an appropriate diffusion coefficient. But for our case on the sphere with a uniform

latitude/longitude grid, it may well be the situation that K is too large near the poles but

too small near the equator. When K is very large, the 2Ax waves just change sign instead

of being smoothed (see (34)).

Figure 2a is the geopotential field of the initial conditon described in McDonald and Bates

(1989), consisting of spherical harmonic wave number one. Figure 2b is the same field

with 2Ax waves embedded. The result of applying the semi-implicit fourth order diffusion

equation, with r4 so chosen that the 2Ax waves at the equator are eliminated, is shown in

Figure 2c. The noise away from the equator remains. Careful examination shows that sign

changes occur. Figure 2d represents the result after applying the implicit diffusion with

the same value for the diffusion coefficient. Now the noise near the poles is damped, but

some noise near the equator remains. Figure 2e is obtained by applying the semi-implicit

formulation followed by the implicit formulation. It is almost identical to the original field

in Figure 2a. Therefore, the use of the semi-implicit formulation of the diffusion equation

alone should be avoided. The implicit formulation is preferred. If some noise near the
equatorneedsto bedamped,thesemi-implicitformulation,in addition to the implicit one,

canbeusedif the additionalcomputationalcost is not a concern.

4 Application to A Semi-Lagrangian Shallow Water Model

In this section, we present some results of applying the fourth order diffusion to the shallow

water model based on the momentum formulation, termed the (u, v) model as in Bates et

al. (1994). This is the 2-D version of the NASA/GLA semi-Lagrangian GCM (Bates et al.,

1993).

4.1 Initial Condition of McDonald and Bates (1989)

For a resolution of 512x257 on a uniform latitude/longitude grid and a timestep of one hour,

a minimum value of 0.03 for the uncentering parameter is needed to suppress noise in the

case of the geostrophically balanced initial condition of McDonald and Bates (1989). The

initial geopotential is defined as:

¢ (A, 0) = ¢ + 2f_av0sin 3 t?cos 0sin A

where f_ and a are, respectively, the earth's rotation and radius. The values ¢ = 5.7879 × 104

m2/s and v0 = 20 m/s are used. The geopotential at 50 days of the integration is presented

in Figure 3a.

Now we replace the uncentering by the fourth order diffusion. In order to remove the residual

noise near the equator, we apply the semi-implicit formulation in addition to the implicit

one as discussed in the previous section. The diffusion coefficient for the semi-implicit

formulation takes a value of 1.32 × 10 is m4/s, which is just large enough to completely

eliminate the 2Ax waves at the equator. The value needed for the implicit formulation is

found to be dependent on the way the diffusion is applied. If the diffusion is applied to the

geopotential only, a minimum value of 9.25 × 1015 m4/s is needed. The geopotential at 50

days of the integration thus obtained is shown in Figure 3b. If the diffusion is applied only

to the divergence, the minimum value needed is 11.825 × 10 Is m4/s. The geopotential at 50

days of the corresponding integration is shown in Figure 3c. Remember that the divergence

is not a basic prognostic model variable, rather it is an intermediate variable. This may be

the reason why a relatively larger minimum value is needed. If we apply the diffusion to

both the divergence and the geopotential, the minimum value for the diffusion coefficient

becomes 5.28 × 10 l_ m4/s, which is smaller than the value used in either of the previous

two cases. The geopotential at 50 days of the integration is now presented in Figure 3d.

The results shown in Figures 3b-d, which were obtained by applying the diffusion in various

ways, demonstrate less damping of the long waves than is the case when using uncentering

9
(Figure 3a). Of the three cases (diffusion applied to geopotential alone, to divergence

alone, and to both geopotential and divergence), it seems that the last two cases have the

least undesirable damping effects. In addition, the results in Figures 3b-d all demonstrate

the same phase speed, but the predicted geopotential of the uncentering case shows faster

westward propagation of the long wave. It has been observed in our numerical results

(not shown) that the westward phase speed monotonically increases as the uncentering

parameter increases. This clearly demonstrates the negative impact of the uncentering on

the long waves in both amplitude and phase. Improvement is obtained by replacing it by

the fourth order diffusion.

Examination of the evolution of the global mean geopotentiai indicates that the fourth order

diffusion, regardless of the magnitude of the diffusion coefficient, has almost no impact on

global mass conservation.

4.2 Real Data Case

The ECMWF analysis of 0000 UTC, 1 December 1992, is chosen. Again, the horizontal

resolution is 512x257 points on a uniform latitude/longitude grid, but the timestep is chosen

to be 30 minutes. The data are initialized using a digital filter (Lynch and Huang, 1992),

with a cutoff period of 6 hours. The initialized geopotential is shown in Figure 4a.

Without any diffusion, a value of 0.05 is needed for the uncentering parameter. The geopotential at day 3 of the integration is shown in Figure 4b. When we apply the fourth order

diffusion to both geopotential and divergence, the model can be integrated stably without

uncentering. The diffusion coefficient for the semi-implicit formulation is 1.32 x 10is ra4/s,

and that for the implicit formulation is 3.96 x 10 is rn4/s. However, there is still some residual noise near the equator, which can be eliminated most effectively by using a small value

of uncentering parameter (0.015 in this case). Figure 4c is the geopotential at day 3 of the

integration thus obtained. The improvement of the representation of the resolved features

in Figure 4c over Figure 4b is recognizable (The high in the upper left part of the model

domain is stronger in Figure 4c). The advantages of using the fourth order diffusion should

become more obvious when orography is present, in which case a much larger value for the

uncentering parameter (undesirable) is otherwise needed to control the high wavenumber

noise in mountainous regions.

5 Application to a 3-D Semi-Lagrangian GCM

We have incorporated the fourth-order horizontal diffusion into the NASA/GLA semiLagrangian semi-implicit (SLSI) general circulation model(GCM) with full physics. The

adiabatic part of this GCM is described in Bates et al. (1993) and the version with full

10
physicsis describedin Moorthi et al. (1994). Here, we present some results to illustrate

the use of the fourth-order diffusion.

The version of the model used here has twenty sigma layers in the vertical and a horizontal

resolution of 2 degree latitude by 2.5 degree longitude. The initial state for the experiments

is taken from the NASA/GLA Data Assimilation Office analyses for November 1, 1992,

performed during the 1992 fall SPADE mission, using the GLA Eulerian GCM. The fourth

order diffusion in its implicit formulation is applied to both temperature and moisture. The

diffusion coefficient is taken as K = rain ( °'25-x1°15, 5.0 × 1015). Here al is the value of the \ at

vertical coordinate sigma at the middle of the lth layer. This choice is made particularly

due to the fact that there exists stronger undesired noise in the model integration above the

100 hPa level than at lower levels. The associated e-folding time for the shortest resolvable

waves at the equator is approximately 2 days near the bottom and 0.1 day near the top.

We have performed four ten-day runs, two without diffusion but with values of 0.05 and

0.07 for the uncentering parameter, and two with diffusionand with a value of 0.05 for the

uncentering parameter. Both runs without diffusion show significant noise at the upper

level of the model. An example is shown in Figure 5, where the temperature field at 50

hPa level is shown for the forecasts at days 3, 5, 7, and 10 respectively for the run with

the uncentering parameter value of 0.05. It can be seen that the temperature field is noisy

and the noise progressively increases near the equator as the integration proceeds. The

corresponding field for the forecast with diffusion is shown in Figure 6. It can be seen that

the inclusion of diffusion has effectively reduced the high wavenumber noise. Of course, one

can obtain even better smoothing by increasing the value of the diffusion coefficient. The

results for the uncentering parameter value 0.07 are similar to those shown above. Because

the noise in the upper level equatorial region (which may be the result of the gravity waves

generated by convective processes), cannot be efficiently eliminated through judicial choice

of the uncentering parameter, inclusion of a higher order diffusion becomes very important.

6 Concluding Remarks

A direct solver for implicit fourth order horizontal diffusion on the sphere has been developed. It has been made efficient by vectorizing over the FFT wave spectrum. Extension

of the formulation to the solution of even higher order diffusion is straightforward, and the

additional coding effort needed is trivial.

Applying the solver to a global shallow water semi-Lagrangian model leads to positive

results. It gives a more accurate representation of larger scales in both amplitude and

phase than does the use of the uncentering method of damping. It is found that horizontal

diffusion has almost no impact on the global mass conservation. The solver is also applied

to the 3-D NASA/GLA semi-Lagrangian semi-implicit GCM and found to be helpful in

11
controllingthehighwavenumbernoise. The solver takes approximately 4% of the total CPU

time on a single processor of the CRAY C98 computer when applied to the temperature

field in the adiabatic version of the 3-D GCM.

Acknowledgements. This research was supported by the NASA Global Atmospheric Modeling and Analysis Program under Grant 578-41-16-20

12
APPENDIX A

Fortran Program of The Solver

Page 1

SUBROUTINE DIRSOL4(IMA, IM,JM,AE, INITI,INIT2,DIFF, FI,FO,XK,DT, IMPL)

**************************************************************************

C Subroutine DIRSOL4

C MITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994

C

C PURPOSE: This routine solves the fourth-order diffusion equation on the

C sphere implicitly or semi-implicitly (Crank-Nicolson scheme), that is,

C to solve the following equation:

C

C (\partial phi over \partial t) = - K \delta^4 (phi)

C

C where K is the diffusion coefficient.

C

C Implicit Soluton: (phi^(n+l)-phi^(n))/DT = - K \del^4 (phi^(n+l))

C

C Semi-Implicit: (phi^(n÷l)-phi^(n))/DT =

C - K (1/2) \del^4 (phi^(n+l)+phi^(n))

C

C Each of the above discrete equation converts to the solution

C of a linear fourth-order partial differential equation on the shpere

C of the following form:

C

C \delta^4 (phi)^(n+l) - r*phi^(n+l) = f

C

C The strategy of solution is similar to that in Lindzen and Kuo (1969,

C Mon. Wea. Rev., p732-734), which is summarized as:

C

C i) . Convert the \del^4-equation into a set of two \del^2 equations.

C Now the model variables become a vector of

C (\del^2(phi)^(n+l), (phi)^(n+l)) ;

C

C 2). Since the coefficients are only latitude-dependent, FFT along the

C latitudes are performed to transform the second-order vector PDE

C to a series of second-order vector ODE's for the FFT coefficients;

C

C 3). Solve the set of second-order ODE's;

C

C 4). Perform inverse FFT to obtain the solution (phi)^(n+l) on gridpoints.

C

C

C CALLED FROM: Any routine which applies this smoother.

c

C SYSTEM ROUTINES USED:

C FFTFAX: Cray routine which initializes the FFT routines;

C RFFTMLT: Cray FFT routine (Cray Manual: SR-2081 8.0)

C

C SUBROUTINES CALLED:

C CAL_R: calculate the coefficient r from K and and DT;

C FORCE: calculate the forcing term f for the \del^4 equation;

C S20DE: solve the set of the second-order ODE's for FFT

C coefficients k > 0;

C SOMEAN:solve for the zonal mean;

C WRAPD: wrap ghost points along the latitudes.

C

C VARIABLES :

INPUT:

AE-radius of earth (meters);

DIFF-logical variable meaning discrete FFT when it is .TRUE.;

DT-time step (seconds)

FI-input field before the smoothing (phi)^(n);

IM-number of grid points along latitudes;

IMA-dimension in X-direction (=IM + # of ghost points);

IMPL-fully implicit scheme (=I) and semi-implicit scheme (=0) ;

INITl-logical variable (need to be TRUE);

INIT2-1ogical variable (need to be TRUE);

XK-diffusion coefficient K (m^4/s).

13
Page 2

C

C

c

C

C

c

c

c

C

C

C

C

C

C

C

C

C

C

C

C

C

C

c

C

c

C

C

C

C

c

C

C

C

C

C

C

C

C

C

C

OUPUT:

FO-output field after smoothing (phi)^(n+l).

LOCAL:

A-array storing the FFT coefficinets;

AM-coefficient array for the 2nd-order ODE (for the zonal mean)

at j-I point;

AMK-coefficient array for the 2nd-order ODE for the FFT coefficients

at j-I point;

B-working array for the FFT;

BM-coefficient array for the 2nd-order ODE (for the zonal mean)

at j point;

BMK-coefficient array for the 2nd-order ODE for the FFT coefficients

at j point;

CM-coefficient array for the 2nd-order ODE (for the zonal mean)

at j+l point;

CMK-coefficient array for the 2nd-order ODE for the FFT coefficients

at j+l point;

CPH-COS(theta) where theta is the latitude;

DNl-real part of the forcing term for the 2nd-order ODE or forcing

for the equation for the zonal mean;

DNIK-real part of the forcing term for the 2nd-order ODE or forcing

to the equation for the FFT coefficients;

DN2K-imaginery part of the forcing term for the 2nd-order ODE or forcing

to the equation for the FFT coefficients;

F-array storing the forcing term for the \del^4 equation;

FM-zonal mean of the input field (phi)^(n);

FNl-the solution to the 2nd-order ODE for the zonal mean;

FNIK-real part of the solution to the 2nd-order ODE for the FFT

coefficients;

FN2K-imaginery part of the solution to the 2nd-order ODE for the FFT

coefficients;

IX-working array for the FFT routine;

TR-working array for the FFT routine.

U-coeffcients for intermediate calculation;

V-coeffcients for intermediate calculation;

W-coeffcients for intermediate calculation;

WVSQ-array storing coefficients for the X-direction derivative term

in the \del^2 equation;

INTEGER FORWARD, BACKWARD

PARAMETER(FORWARD = -i, BACKWARD = i)

DIMENSION FI(IMA,JM), FO(IMA, JM)

DIMENSION IX(IM), TR(3*(IM+2)+I), A(IM+2,JM-2), B(IM*2,JM-2)

DIMENSION AM(2,2,JM), BM(2,2,JM), CM(2,2,JM), F(IMA,JM)

#, DNI(2,JM), DN2(2,JM), FNI(2,JM), FN2(2,JM)

#, FM(JM), CPH(JM,2), U(JM), V(JM), W(JM)

#, WV(IM/2), WVSQ (IM/2)

DIMENSION AMK(IM/2,2,2,JM), BMK(IM/2,2,2,JM), CMK(IM/2,2,2,JM)

#, DNIK(IM/2,2,JM), DN2K(IM/2,2,JM)

#, FNIK(IM/2,2,JM), FN2K(IM/2,2,JM)

LOGICAL INITI, INIT2, DIFF, FIRST

PI = ACOS(-I.0)

TWOPI = 2.0 " PI

PI02 = PI/2.0

IMB2 = IM/2

IJM = (IMB2-1)*JM

IJMM = IJM*22 + JM

IJM2 = IJMM + (IM+2)'(JM-2)

14
Page 3

CALL CAL R(XK, DT, R, IMPL)

CALL FORCE (IMA, IM, JM, FI, F, AE, R, IMPL)

CALL WRAPD (IMA, IM, JM, F)

IF (INITI) THEN

CALL FFTFAX (IM, IX, TR)

JMMI = JM - 1

JMM2 = JM - 2

IMP1 = IM ÷ 1

IMP2 = IM + 2

LEN = IMB2 - 1

FIM = 1.0 / FLOAT(IM)

DLM = TWOPI*FIM

DPH = PI/FLOAT(JMMI)

RDLM = 1.0/DLM

RDPH = 1.0/DPH

DY = AE*DPH

AESQ = AE*AE

DO J = I, J_

TEM = - PIO2 + (J-I)*DPH

CPH(J,I) = COS(TEM+0.5*DPH)

CPH(J,2) = COS(TEM)

ENDDO

CPH(I,2) = 0.0

CPH(JM,2) = 0.0

FIRST = .TRUE.

INIT2 = .TRUE.

ENDIF

IF (INIT2) THEN

DO J = 2, JMMI

U(J) = CPH(J,I)/CPH(J,2)

W(J} = CPH(J-I,I)/CPH(J02)

V(J) = - 1.0 * (U(J)+W(J))

ENDDO

IF (DIFF) THEN

DO I = I, IMB2

WVSQ(I) = (SIN(0.5*I*DLM)*( RDLM*2"0))t*2

ENDDO

ELSE

DO I = I, IMB2

WVSQ(I) = I'I

ENDDO

ENDIF

ENDIF

A(:, :) = 0.0

DOJ = I, JM

FM(J) = 0.0

DOI = i, IM

FM(J) = FM(J) + F(I,J)

ENDDO

ENDDO

DOJ = I, JM

FM(J) = FM(J)*FIM

ENDDO

DO J = 2, JMMI

DOI = i, IM

A(I,J-I) = F(I,J)

ENDDO

FM(J)

15
Pago 4

A(IM÷I,J-I) = A(I,J-I)

A(IM+2,J-I) = A(2,J-I)

ENDDO

CALL RFFTMLT(A,B,TR, IX, I,IMP20IM, JMM2,FORWARD)

A(I, :) = 0.0

A(2, :) = 0.0

A(IMP2, :) = 0.0

DO J = 2, JMMI

DO I = i, IMB2

AMK(I,I,I,J) = W(J)

AMK(I,2,2,J) -- W(J)

AMK(I,I,2,J) = 0.0

AMK(I,2,I,J) = 0.0

CMK(I,I,I,J) = U(J)

CMK(I,2,2,J) = U(J)

CMK(I,I,2,J) = 0.0

CMK(I,2,I,J) = 0.0

BMI4(I,I,I,J) = - (DY**2)'WVSQ(I)/(AE*CPH(J,2))**2

÷ V(J)

BMK(I,I,2,J) = R*(DY**2)

BI_IK(I,2,I,J) = - DY*'2

BMK(I,2,2,J) = - (DY**2)*WVSQ(I)/(AE*CPH(J,2))**2

+ V(J)

DNIK(I,I,J) = (DY**2)'A(I+I+I,J-I)

DNIK(I,2,J) = 0.0

DN2K(I,I,J) = (DY**2)*A(I+I+2,J-I)

DN2K(I,2,J) = 0.0

ENDDO

ENDDO

DNIK(:, :,i) = 0.0

DN2K(:, :,i) = 0.0

AMK(:,I,I,I) = 1.0

AMK(:,2,2,1) = 1.0

AMK(:,I,2,1) = 0.0

AMX(:,2,1,1) = 0.0

BMK(:, , :,i) = 0.0

DNIK(: :,JM) = 0.0

DN2K(: :,JM) = 0.0

AMK(:,:,:,JM) = 0.0

BMK(:,I,I,JM) = 1.0

BMK(:,2,2,JM) = 1.0

BMK(:,!,2,JM) = 0.0

BM_(:,2,I,JM) = 0.0

CALL S2ODE(IMB2,DNIK,DN2K, FNIK,FN2K,AMK,BMK, CMK,JM)

DO J = 2, JMMI

DO I = I, IMB2

A(I+I+I, J-l) = FNIK(I,2,J)

A(I+I+2, J-l) = FN2K(I,2,J)

ENDDO

ENDDO

DO J = 2, JMMI

16
P&go 5

AM(I,I,J) = W(J)

AM(2,2,J) = W(J)

AM(I,2,J) = 0.0

AM(2,I,J) = 0.0

CM(I, l,J) = U(J)

CM(2,2,J) = U(J)

CM(I,2,J) = 0.0

CM(2,I,J) = 0.0

BM(I,I,J) = + V(J)

BM(I,2, J) = R* (DY*_2)

BM(2,I,J) = - DY**2

BM(2,2,J) = BM(I,I,J)

DNI(I,J) = (Dy**2)*FM(J)

DNI(2,J) = 0.0

ENDDO

DNI(I,I) = (Dy**2)*FM(1)

DNI (2, i) = 0.0

AM(I,1,1) = - 4.0

AM(I,2,1) = R-(DY'*2)

AM(2,1,1) = - (DY**2)

AM(2,2,1) = - 4.0

BM(I,I,I) = 4.0

BM(I,2,1) = 0.0

BM(2,1,1) = 0.0

BM(2,2,1) = 4.0

DNI(I,JM) = (DY*-2)*FM(JM)

DNI(2,JM) = 0.0

AM(I,I,JM) = 4.0

AM(I,2,JMJ = 0.0

AM(2,I,JM) = 0.0

AM(2,2,JM) = 4.0

BM(I,I,JM) = - 4.0

BM(I,2,JM) = R*(DY**2)

BM(2,I,JM) =- (DY**2)

BM(2,2,JM) = - 4.0

FNI(:, :) = 0.0

CALL SOMEAN(DNI, FNI, AM, BM, CM, JM)

CALL RFFTMLT (A, B, TR, IX, I, IMP2, IM, JMM2, BACKWARD)

DO J = i, JMM2

DOI=I, IM

FO(I,J+I) = A(I,J) ÷ FNI(2,J+I)

ENDD0

ENDDO

[901=1, IM

FO(I,I) = FNI(2,1)

FO(I,JM) = FNI(2,JM)

KNDDO

CALL WRAP1 (IMA, IM, JM, FO)

RETURN

KND

SUBROUTINE CAL_R(XK, DT, R, IMPL)

AB = (I.0+FLOAT(IMPL))'XK'DT

17
Page 6

R = 2.01AB

RETURN

END

SUBROUTINE DEL2(IMA, IM, JM,A, DA,AE)

C Subroutine DEL2

C WRITTEN BY: Yong Li, GSFC code 910.3, Apt 1994

C

C PURPOSE: calculate the \del^2(A) on the sphere.

C

C CALLED BY:

C DEL4

C

C SUBROUTINES CALLED: None

C

C VARIABLES :

C

C INPUT :

C A-input field;

C AE-radius of earth;

C IM-number of total grid points in X-direction;

C IMA-dimension size in X-direction;

C JIM-number of total grid points in Y-direction;

C

C OUTPUT:

C DA-\del^2 of A

dimension a(ima,jm), da(ima,jm), dx(jm)

dimension al(500), a3(500), bp(500), bm(500)

impl = im + 1

imp2 = im + 2

pi = acos(-l.0)

twopi = 2.0*pi

dlm = _wopi/float(im)

dph = pi/float(jm-l)

dy = ae*dph

a2=l.0d0/dy**2

do j = 2, jm-I

the_a = - pi*0.5 ÷ float(j-l)'dph

bp(j) = cos(theta+0.5*dph)/cos(theta)

bm(j) = cos(theta-0.5*dph)/cos(theta)

dx(j) = ae'dlm*cos(theta)

al(j) = 1.0/dx(j)'*2

a3(j) = 2.0*al(j)+a2*(bp(j)+bm(j))

enddo

30

20

do 20 i=2,impl

do 30 j=2,jm-i

da (i, j) =al (j) * (a (i+l, j) +a(i-l, j) )

# +a2*(bp(j)*a(i,j+l)÷bm(j)*a(i,j-l))

# -a3 (j)*a(i, j)

continue

continue

do j = 2, jm-i

da(l,j) = da(im+l,j)

do i = im+2, ima

da(i,j) = da(i-im, j)

enddo

enddo

Laplace at the south pole

AAC = 0.0

18
Page 7

4O

5O

C

C

C

60

do 40 i = 2, impl

AAC = AAC + a(i,2)

continue

da(l,l) = (2.0/dy)**2*(AAC/float(im)-a(l,l))

do 50 i = 2, ima

da(i,l) = da(l,l)

continue

Laplace at the north pole

AAC = 0.0

do 60 i = 2, imp l

AAC = AAC + a(i,jm-l)

continue

da(ima,jm) = (2.0/dy)**2*(AAC/float(im)-a(l,jm))

do i = i, ima-I

da(i,jm) = da(ima, jm)

enddo

return

end

SUBROUTINE DEL4(IMA, IM, JM,A, DA,AE)

C Subroutine DEL4

C WRITTEN BY: Yong Li, GSFC code 910.3, Apr 1994

C

C PUEPOSE: calculate the \del^4(A) by applying a conservative \del^2 twice.

C

C CALLED FROM:

C FORCE: calculate the forcing term

C

C SUBROUTINES CALLED:

C DEL2: \del^2 operator;

C WRAPD: wrap up the ghost points in X-direction.

C

C VARIABLES:

C

C

C

C

C

C

C

C

C

C

INPUT:

A-input field;

AE-radius of earth;

IM-number of total grid points in X-direction;

IMA-dimension size in X-direction;

JM-number of total grid points in Y-direction.

OUTPUT:

DA- \del^4 (A)

DIMENSION A(IMA, JM), DA(IMA,JM), B(IMA, JM)

CALL DEL2(IMA, IM, JM,A,B,AE)

CALL WRAPD(IMA, IM,JM, B)

CALL DEL2 ( IMA, IM, JM, B, DA,AE)

CALL WRAPD (IMA, IM, JM, DA)

RETURN

END

SUBROUTINE FORCE (IMA, IM, JM, A, F, AE, R, IMPL)

*****************************************************************

C Subroutine FORCE

C WRITTEN BY: Yong Li, GSFC code 910.3, Apt 1994

C

19
Pmge 8

C PURPOSE: calculate the forcing term for the \del^4-equation.

C

C CALLED FROM: DIRSOL4

C

C SUBROUTINES CALLED:

C DEL4: calculate the \del^4 operator of field A;

C WRAPD: wrap up the ghost points in X-direction;

C

C VARIABLES:

C

C INPUT:

C A-array storing the field to be smoothed;

C AE-radius of earth;

C IM-number of total grid points in X-direction;

C IMA-dimension size in X-direction;

C IMPL-for implicit (=i), for semi-implicit (=0);

C JM-number of total grid points in Y-direction;

C R-coefficient for the zero-order term in the \del^4 equation;

C

C OUTPUT:

C DA-array storing \del^4(A);

C F-array storing the forcing term for the \del^4 equation

DIMENSION A(IMA,JM), DA(IMA,JM), F(IMA, JM)

IF(IMPL.EQ.I) THEN

F(:, :) = R'A(:, :)

ELSE

CALL DEL4(IMA, IM,JM,A, DA,AE)

CALL WRAPD(IMA, IM, JM,DA)

F(:, :) = R'A(:, :) - DA(:, :)

ENDIF

RETURN

END

SUBROUTINE INVERSEI(A,B)

C Subroutine INVERSE1

C

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994

C

C PURPOSE: Find the inverse of a 2x2 matrix

C

C VARIABLES:

C INPUT:

C A-input matrix to be inverted;

C OUTPUT:

C B-inverse of matrix A;

DIMENSION A(2,2), B(2,2)

DD = A(I,I)*A(2,2)-A(I,2)'A(2,1)

B(I,I) = A(2,2)/DD

B(I,2) = - A(I,2)/DD

B(2,1) = - A(2,1)/DD

B(2,2) = A(I,I)/DD

RETURN

END

SUBROUTINE INVERSE2(LEN,A,B)

C Subroutine INVERSE2

C

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994

C

2O
Page 9

C PURPOSE: Find the inverse of a a number (LEN) of 2x2 matrix

C

C VARIABLES :

C INPUT :

C A-input of LF--N 2x2 matrices to be inverted simultaneously;

C OUTPUT :

C B-inverse the LEN 2x2 matrices A;

Dimension A(LEN, 2,2), B(LEN, 2,2), DD(LEN)

DO KK = I, LEN

DD(KK) = A(KK, I,I)*A(KK,2,2)-A(KK, I,2)*A(KK,2,1)

ENDDO

DO KK = I, LEN

B(KK, I,I) = A(KK,2,2)/DD(KK)

B(KK01,2) = - A(KK, I, 2) /DD(KK)

B(KK,2,1) = - A(KK,2,1)/DD(KK)

B(KK,2,2) = A(KK, I, I) /DD(KK)

ENDDO

RETURN

END

SUBROUTINE MULTAB (A, B, C, N)

C Subroutine MULTAB

C

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994

C

C PURPOSE: Multiplication of two NxN matrix C = A B

C

C VARIABLES :

C INPUT :

C A-input mmtrix;

C B-input matrix;

C OUTPUT :

C C-output matrix.

DIMENSION A(N,N), B(N,N), C(N,N)

DO I = I, N

DOJ = i, N

C(I,J) = 0.0

DOK = I, N

C(I,J) = C(I,J) + A(I,K)*B(K,J)

ENDDO

ENDDO

ENDDO

RETURN

END

SUBROUTINE MULTABK (LEN, A, B, C, N)

C Subroutine MULTABK

C

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994

C

C PURPOSE: Multiplication of a series of two NxN matrix C = A B

C

C VARIABLES :

C INPUT :

C A-input of LEN 2x2 matrices;

C B-input of LEN 2x2 matrices;

C OUTPUT :

C C-output matrix.

****************************************************************

DIMENSION A(LEN,N,N), B(LEN, N,N), C(LEN,N,N)

#, WRKI (LEN) , WRK2 (LEN)

21
Page 10

DO I = i, N

DOJ = I, N

C( : , I,J)

DOK = I,

= 0.0

N

WRKI(:) = C(:,I,J)

DO KK = i, LEN

WRK2(KK) = A(KK, I,K)*B(KK,K,J)

ENDDO

C(:,I,J) = WRKI(:) + WRK2(:)

ENDDO

ENDDO

ENDDO

RETURN

END

SUBROUTINE MULTAV(A, V, W, N)

C Subroutine MULTAV

C

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994

C

C PURPOSE: Multiplication of an NxN matrix A and a vector V: W = A V

C

C VARIABLES :

C INPUT:

C A-input matrix;

C B-input vector;

C OUTPUT :

C C-output vector.

DIMY, NSION A(N,N), V(N), W(N)

DOI = I, N

W(I) = 0.0

DOJ = i, N

W(I) = W(I) + A(I,J)*V(J)

ENDDO

ENDDO

RETURN

END

SUBROUTINE MULTAVK (LEN, A, V, W, N)

C Subroutine MULTAVK

****************************************************************

C

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994

C

C PURPOSE: Multiplication of an NxN matrix A and a vector V: W = A V

C

C VARIABLES :

C INPUT:

C A-input of LEN 2x2 matrices;

C B-input of LEN vectors;

C OUTPUT :

C C-output vector.

****************************************************************

DIMENSION A(LEN,N,N), V(LEN, N), W(LEN,N), WRK(LEN)

DOI = i, N

W(:,I) = 0.0

DO J=l, N

WRK( :) = W(:,I)

W(:,I) = WRK(:) + A(:,I,J)*V(:,J)

ENDDO

ENDDO

RETURN

END

22
Page 11

SUBROUTI_ S20DE(LEN, DNI, DN2, FNI, FN2, AM, BM, CM, JM)

C Subroutine S2ODE

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apt 1994

C

C PURPOSE: This routine solves the second-order ODE which originally results from

C the \del^4 equation on the sphere. But appplication is not so restricted

C as far as the discretized equation is in the format as (i).

C The unknown to be solved is a 2xl vector (Lindzen and Kuo, 1969, MWR).

C Therefore this routine is used here to solve for the FFT coefficients (k > 0)

C in the \del^4 direct solver if the the B.C's are chosen to be zero for the poles.

C

C The discretized form of the ODE is:

C

C A(j) X(j-l) + B(j) X(j) + C(j) X(j+l) = F(j) ...... (I)

C (for j = 2, 3, 4 ....... JM-I)

C

C The boundary conditions are:

C

C A(1) X(1) + B(1) X(2) = F(1) ........ (2.1)

C

C and A(JM) X(JM-I) + B(JM) X(JM) = F(JM) ........ (2.2)

C

C The coefficient A, B and C are all real. The forcing F and hence the

C solution X have both real and imaginery parts.

C

C CALLED FROM:

C DIRSOL4

C

C SUBROUTINES CALLED:

C INVERSE1: inversion of one 2x2 matrix;

C INVERSE2: inversion of a number of 2x2 matrices;

C MULTAB: multiplication of two matrices;

C MULTABK: multiplication of a number of two matrix set;

C MULTAV : multiplication of a matrix and a vector.

C MULTAVK : multiplication of a number of a matrix-vector sets.

C

C VARIABLES :

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

INPUT:

AM-array storing the coefficient matrix A as in (i);

BM-array storing the coefficient matrix B as in (i);

CM-array storing the coefficient matrix C as in (i);

DNl-array storing the real part of the vector of the forcing

term F in (i);

DN2-array storing the imaginery part of the vector of the forcing

term F in (I) ;

JIM-number of total grid points including the two boundary points.

OUTPUT:

FNl-array storing the real part of the solution vector X of (i);

FN2-array storing the imaginery part of the solution vector X of (i);

LOCAL:

ALPHA-array storing intermediate coefficient matrix alpha

(Lindzen and Kuo, 1969);

BETAI-real part of vector beta (Lindzen and Kuo, 1969);

BETA2-imaginery part of vector beta (Lindzen and Kuo, 1969);

VRKl-working array;

VRK2-worklng array;

VRK3-worklng array;

VRK4-working array;

VRK5-working array;

VRK6-working array;

WRKl-working array;

WRK2-working array;

WRK3-worklng array;

WRK4-working array;

23
DIMENSION DNI(LEN, 2,JM), DN2(LEN, 2,JM), FNI(LEN, 2,JM)

#, FN2 (LEN, 2, JM)

D IMENS I ON

DIMENSION

DIMENSION

#,

#.

#,

WRKI(:, :, :) = AM(:,:,:,1)

CALL INVERSE2(LEN, WRKI, WRK2)

WRKI(:, :, :) = - WRK2(:, :, :)

WRK4(:, :, :) = BM(:,:,:,I)

CALL MULTABK(LEN, WRKI, WRK4, WR_K3, 2)

ALPHA(:,:,:,1) = WRK3(:, :, :)

VRKI(:,:) = DNI(:,:,I)

VRK2(:,:) = DN2(:,:,I)

CALL MULTAVK(LEN, WRK2, VRKI, VRK3, 2)

CALL MULTAVK (LEN, WRK2, VRK2. VRK4, 2 )

BETA1(:, :,I) = VRK3(:, :)

BETA2(:,:,1) = VRK4(:,:)

DO J = 2, JM-I

WRKI(:, :,:) = ALPHA(:, :,:,J-l)

WRK2(:,:,:) = AM(:,:,:,J)

CALL MULTABK(LEN, WRK2, WRKI, WRK3, 2)

WRKI(:, :, :) = WRK3(:, :, :) + BM(:, :,:,J)

CALL INVERSE2 (LEN,WRKI,WRK2)

WRKI(:, :, :) = CM(:, :, :,J)

CALL MULTABK(LEN, WRK2, WRKI, WRK3, 2)

ALPHA(:,:,:,J) = - WRK3(:, :, :)

WRKI(:, :, :) = AM(:, :, :,J)

VRKI(:,:) = BETA1(:, :,J-l)

VRK2(:,:) = BETA2(:, :,J-l)

CALL MULTAVK(LEN, WR.KI, VRKI, VRK3, 2)

CALL MULTAVK(LEN, WR/(I, VRK2, VRK4, 2)

VRK3(:,:) = DNI(:,:,J) - VRK3(:,:)

VRK4(:,:) = DN2(:,:.J) - VRK4(:,:)

CALL MULTAVK(LEN, WRK2, VRK3, VRKI, 2)

CALL MULTAVK(LEN, WRK2, VR/(4, VRK2. 2)

BETA1(:, :,J) = VRKI(:, :)

BETA2(:, :,J) = VRK2(:. :)

ENDDO

VRKI(:, :) = DNI(:,:,JM)

VRK2(:, :) = DN2(:, :,JM)

VRK3(:, :) = BETA1(:, :,JM-I)

VRK4(:,:) = BETA2(:,:,JM-I)

WRKI(:.:.:) = AM(:, :, :,JM)

AM(LEN, 2,2,JM), BM(LEN, 2,2,JM), CM(LEN, 2,2,JM)

ALPHA(LEN, 2,2, JM) , BETA1 (LEN, 2, JM) , BETA2 (LEN, 2, JM)

WRKI(LEN, 2,2), WRK2(LEN, 2,2), WRK3(LEN, 2,2)

WRK4 (LEN, 2,2)

VRKI(LEN, 2), VRK2(LEN, 2), VRK3 (LEN, 2), VRK4(LEN, 2)

VRK5 (LEN, 2) , VRK6 (LEN, 2)

Page 12

24
Pa_ 13

CALL MULTAVK(LEN, WRKI, VRK3, VRK5, 2)

CALL MULTAVK(LEN, WRKI, VRK4, VRK6, 2)

VRKI(:, :) = DNI(:,:,JM) - VRK5(:, :)

VRK2(:,:) = DN2(., :,JM) - VRK6(:, :)

WRKI(:,:,:) = ALPHA(:,:,:,JM-I)

WRK2(:, :, :) = AM(:, :, :,JM)

CALL MULTABK(LEN, WRK2, WRKI, WRK3, 2)

WRK4(:, :, :) = BM(:,:,:,JM) + WRK3(:, :, :)

CALL INVERSE2(LEN, WRK4, WRKI)

CALL MULTAVK(LEN, WRKI, VRKI, VRK3, 2)

CALL MULTAVK(LEN, WRKI, VRK2, VRK4, 2)

FNI(:,:,JM) = VRK3(:,:)

FN2(:,:,JM) = VRK4(:,:)

DO J = JM-I, I, -i

WRKI(:, :, :) = ALPHA(:,:, :,J)

VRKI(:, :) = FNI(:, :,J+l)

VRK2(:, :) = FN2(:, :,J+l)

CALL MULTAVK(LEN, WRKI, VRKI, VRK3, 2)

CALL MULTAVK(LEN, WRKI, VRK2, VRK4, 2)

FNI(:, :,J) = VRK3(:, :) + BETA1(:, :,J)

FN2(:,:,J) = VRK4(:,:) + BETA2(:,:,J)

ENDDO

RETURN

END

SUBROUTINE SOMEAN(DNI, FNI, AM, BM, CM, JM)

C Subroutine SOMEAN

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994

C

C PURPOSE: This routine solves the second-order ODE which originally results from

C the \del^4 equation on the sphere. But appplication is not so restricted

C as far as the discretized equation is of the following format.

C The unknown to be solved is a 2xl vector (Lindzen and Kuo, 1969, MWR).

C This routine is used here to solve for the zonal mean, equivalent to FFT

C coefficients (k = 0) in the \del^4 direct solver.

C

C The discretized form of the ODE is:

C

C A(j) X(j-l) + B(j) X(j) + C(j) X(j÷l) = F(j) ...... (I)

C (for j = 2, 3, 4 ....... JM-I)

C

C The boundary conditions are:

C

C A(1) X(1) ÷ B(1) X(2) = F(1)

C

C and A(JM) X(JM-I) + B(JM) X(JM) = F(JM)

C

C The coefficient A, B and C are all real.

C solution X also real.

C

C CALLED FROM:

C DIRSOL4

C

C SUBROUTINES CALLED:

C INVERSE: matrix inversion;

C MULTAB: multiplication of two matrix;

C MULTAV: multiplication of a matrix and a vector.

C

C VARIABLES:

........ (2.1)

........ (2.2)

The forcing F and hence the

25
Pagl l&

C

C

C

C

C

C

c

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

INPUT:

AM-array storing the coefficient matrix A as in (i);

BM-array storing the coefficient matrix B as in (i);

CM-array storing the coefficient matrix C as in (I);

DNl-array storing the vector of the forcing term F in (I);

JM-number of total grid points including the two boundary points.

OUTPUT :

FNl-array storing the solution vector X of (1);

LOCAL:

ALPHA-array storing intermediate coefficient matrix alpha

(Lindzen and Kuo, 1969);

BETAI-real part of vector beta (Lindzen and Kuo, 1969);

VRKl-working array;

VRK2-working array;

VRK3-working array;

VRK4-working array;

VRK5-working array;

VRK6-working array;

WRKl-working array;

WRK2-working array;

WRK3-working array;

C WRK4 -working array;

DIMENSION AM(2,2,JM), BM(2,2,JM), CM(2,2,JM)

#, DNI(2,JM), FNI(2,JM), ALPHA(2,2.JM). BETAI(2.JM)

DIMENSION WRKI(2,2), WRK2(2,2), WRK3(2,2), WRK4(2,2]

DIMENSION VRKI{2), VRK2(2), VRK3(2), VRK4(2), VRK5(2)

#, VRK6 (2)

WRKI(:, :) = AM(:, :,i)

CALL INVERSE1 (WRKI, WRK2)

WRKI(:,:) = - WRK2(:,:)

WRK4(:.:) = BM(:,:,I)

CALL MULTAB(WRKI, WRK4,

ALPHA(:,:,1) = WRK3(:,:)

WRK3, 2)

VRKI(:) = DNI(:,I)

CALL MULTAV(WRK2. VRK1. VRK3. 2)

BETA1(:,1) = VRK3(:)

DO J = 2, JM-I

WRKI(:,:) = ALPHA(:, :,J-l)

WRK2(:.:) = AM(:,:,J)

CALL MULTAB (WRK2, WRKI, WRK3 ,2 )

WRKI(:.:) = WRK3(:,:) + BM(:,:.J)

CALL INVERSEI(WRKI.WRK2)

WRKI(:,:) = CM(:,:,J)

CALL MULTAB(WRK2,WRKI,WRK3,2)

ALPHA(:, :,J) = - WRK3(:.:)

WRKI(:.:) = AM(:,:,J)

VRKI(:) = BETAI|:,J-I)

CALL MULTAV(WRKI, VRK1, VRK3, 2)

VRK3(:) = DNI(:,J) VRK3(:)

CALL MULTAV(WRK2. VRK3, VRKI, 2)

26
Pagl 15

BETAI(:,J) = VRKI(:)

ENDDO

VRKI(:) = DNI(:.JM)

VRK3(:) = BETAI(:,JM-I)

WRKI(:, :) = AM(:,:,JM)

CALL MULTAV(WRKI, VRK3, VRK5, 2)

VRKI(:) = DNI(:,JM) VRK5(:)

WRKI(:,:) = ALPHA(:,:.JM-I)

WRK2(:,:) = AM(:,:,JM)

CALL MULTAB (WRK2, WRKI, WRK3, 2 )

WRK4(:,:) = BM(:,:,JM) + WRK3(:,:)

CALL INVERSE1 (WRK4, WRKI)

CALL MULTAV(WRKI, VRKI, VRK3, 2)

FNI(:,JM) = V'R.K3(:)

DO J = JIM-l, I, -I

WRKI(:.:) = ALPHA(:,:,J)

VRKI(:) = FNI(:,J+I)

CALL MULTAV(WRKI, VRKI, VRK3, 2)

FNI(:,J) = VRK3(:) + BETAI(:,J)

ENDDO

RETURN

END

SUBROUTINE WRAP1 (IMA, IM. JM. A)

DIMENSION A( IMA, JM)

DOJ = i. JM

DO I = IM+I, IMA

A(I,J) = A(I-IM,J)

ENDDO

ENDDO

RETURN

END

SUBROUTINE WRAPD (IMA. IM, JM, A)

DIMENSION A( IMA, JM)

DOJ = I, JM

A(I,J) = A(IM+I,J)

DO I = IM+2, IMA

A(I,J) = A(I-IM, J)

ENDDO

ENDDO

RETURN

END

27

References

Bates,J. R., S. Moorthi and R. W. Higgins, 1993: A global multilevel atmospheric model using a vector semi-Lagrangian finite-difference scheme. Part I: Adiabatic formulation.

Mon. Wea. Rev., 121,244-263.

Bates, J. R., Y. Li, A. Brandt, S. F. McCormick and J. Ruge, 1994: A global shallow water

model based on the semi-Lagrangian advection of potential vorticity. Accepted for

publication in Q. J. Roy. Met. Soc.

Gravel, S., A. Staniforth and J. CStd, 1993: A stability analysis of a family of baroclinic

semi-Lagrangian forecast models. Mon. Wea. Rev., 121,815-824.

Jakimow, G., E. Yakimiw and A. Robert, 1992: An implicit formulation for horizontal

diffusion in gridpoint models. Mon. Wea. Rev., 120, 124-130.

Lindzen, R. S. and H.-L. Kuo, 1969: A reliable method for the numerical integration of a

large class of ordinary and partial differential equations. Mon. Wea. Rev., 97,732-

734.

Lynch, P. and Huang, X.-Y., 1992: Initialization of the HIRLAM model using a digital

filter. Mon. Wea. Rev., 120, 1019-1034.

McDonald, A., 1994: Using second, fourth and sixth order implicit horizontal diffusion

to control noise in three dimensional semi-Lagrangian, semi-implicit, limited area,

gridpoint models of the primitive equations. Unpublished Note (available from Irish

Meteorological Service, Glasnevin Hill, Dublin 9, Ireland).

McDonald, A. and J. R. Bates, 1989: Semi-Lagrangian integration of a gridpoint shallow

water model on the sphere. Mon. Wea. Rev., 117, 130-137.

Mesinger, F. and A. Arakawa, 1976: Numerical methods used in atmospheric models, Vol.1.

GARP Publications Series No.17., World Meteorological Organization, Geneva, Switzerland.

Moorthi, S. and R. W. Higgins, 1992: A fast direct solver for a class of two-dimensional

seperable elliptic equations on the sphere. NASA Tech. Memo. 104567.

Moorthi, S. and R. W. Higgins, 1993: Application of fast Fourier transformations to the

direct solution of a class of two-dimensional separable elliptic equations on the sphere.

Mon. Wea. Rev., 121,290-296.

Moorthi, S., R. W. Higgins and J. R. Bates, 1994: A global multilevel atmospheric model

using a vector semi-Lagrangian finite-difference scheme. Part 2: Version with physics.

29

PAGE tANK NOT FILMEL1
Accepted for publication in Mon. Wea. Rev.

Moorthi, S., 1994: Numerical weather prediction experiments with a semi-Lagrangian semiimplicit grid-point general circulation model at NMC. lOth Conference on Numerical

Weather Prediction, July 18-22, 1994, Portland, Oregon.

Raymond, W. H. and A. Garder, 1991: A review of recursive and implicit filters. Mon.

Wea. Rev., 119, 477-495.

Ritchie, H., C. Temperton, A. Simmons, M. Hortal, T. Davies, D. Dent and M. Hamrud,

1994: Implementation of the semi-Lagrangian method in a high resolution version of

the ECMWF forecast model. Mon. Wea. Rev. (in press).

Semazzi, F. H. M. and P. Dekker, 1994: Optimal accuracy in semi-Lagrangian models. Part

I. Linear formulation. Mon. Wea. Rev. (in press).

Shapiro, R., 1971: The use of linear filtering as a parameterization of atmospheric diffusion.

J. Atmos. Sci., 28, 523-531.

Shapiro, R., 1979: Linear filtering on the surface of a sphere. AFGL-TR-79-0263, Environ.

Res. Paper No. 683, Air Force Geophysics Laboratory, Hanscom AFB, MA.

Shuman, F. G., 1957: Numerical methods in weather prediction: II. Smoothing and filtering.

Mon. Wea. Rev., 85, 357-361.

Tanguay, M., E. Yakimiw, H. Ritchie and A. Robert, 1992: Advantages of spatial averaging

in semi-implicit semi-Lagrangian schemes. Mon. Wea. Rev., 120_ 113-123.

Williamson, D. L. and J. G. Olson, 1994: Climate simulations with a semi-Lagrangian

version of the NCAR CCM2. Mon. Wea. Rev., 122, 1594-1610.

3O
1.00

0.80

0.60

0.40

020

000

i °

:°

,i

•°i-"

,;o ,

: • : o

!. .;;

,' i

.... i-.-;:---4 ........ i ......... _ ........ _......... _........ _......... _........ ÷ ........ i ........ * ........ _......... _........ _......... _........ _......

.......... 2nd-order (implicit)

.......... 2nd-order (semi-implicit)

4th-order (implicit)

4th-order (semi-im fliciO ..............

i i i i

E i : i

i i i :

, i i _ i i J i i

2 4 6 8 10 12 14 16 18 20

m

Figure 1: The damping rates as a function of wavenumber for the second and the fourth

order diffusion in both implicit and explicit formulations.

31
(a)

', ,'}'-,'_,:.- - .'r- -.:,".,'-',:_ _' ' :/If_

1 , , _ _ _ l , , i i ,'

.._ _' ._ s _'_ _ _ J _,s " Jl _ I I | I X V._

, , _'.- _-___'-_-_-:---_ ,, : \\"q

\ _ " _ _ _I--- - -- t- -- "--- "_" -- _ _ ' / s I

........... _ ........ -4-- -",---:.--_d_ ---- .... 4---_-z. .... _'----, ....... 4------_--_--'

, I I

i ! '

............ _..................... ±................................. _ .......... _........... q-.......... -i-.......... _......................... i i i

I I I

I I I I

i : ---4----':__ ! i : i , . i - -_--i- - -_-..... i- _ - ,,-

.... _ .......... _'T'r"--r ...... 4-.... I".... =--I-= ° ..... _"'/"T'T'}"_'+'t'T"

i ,, ',_-"-"" "-","_- -_- --:-j'- ".- ._'.",

I \ ',.L',._"---.._____../././_ ,,. _ -_." _--b-_------,,--: _.--..-,---.:- .,.. "., / _____-'-----z,,__ ----'_-_------_'_ : "- _ -i-'-'--Z-._---_.=-_.--4_-----_'-__-----_ _ -

(u)

........ _--t--_ .... _I..... _...q....t. I ...... ' --.L ....... , ...... 4-.... I. ," _, " :- /...-,. ,T, , ' r_', r_ ___-'_--->÷'H'H-t""

I ' ' _.x[._ _.-:_- .._ ..-_, _ / . _ ' I ,11 I I_ k.._ _"M_"x_-'i'-_J.J"J/_Jr'J I I I I

r -"--_'V-- _---=--'---_-_-.-7--%-_ --_-_.--_---_-_-_..... _----_ ..... L...._....'2____2___.____ __/ ....

, ... I, ,,. ....... ." _" -Is " A,_I v_L__IZ2_ .....

.... :.....,.,,,,,,.

F-° "-_r_ -_.--7_.__.-.__.____j ..... -_...._r._.#...r_-_..._._.._ .Y_c'.__._.._._,.=C..%.,,._...._.._...t

I I , _f___-"__,"_'%_ '%_-_ _ _ I , ", ,i, ; ", " -_: ," - "! s_'_'- -:-" _", ':_ _ -

--H-t _LL-LLLI._..L_._..,L._J.d.i:'_.__._._..'_.L.'__'I._.-._L._.._" ".,.._

' . 111,_,.-_,_.'2., "'.'!"'.,_,,'_, _'"1,

Figure 2: (a) Spherical harmonic wave number one; (b) Same as (a), but with 2Ax waves

embedded; (c) After applying the semi-implicit formulation of fourth order diffusion to (b);

(d) After applying the implicit formulation to (b); (e) After applying both the semi-implicit

and the implicit formulations to (b).

32
(c)

I, - ;i-_--- =---.-- _:---- _- -.- r- :.."/ \_._____/___.._/

...... "J'-- - * - l -_ _- --'-- --- l -- - _:l--'- =,:7T,._-'-_; - -- - 4 "_ ...... t ......... |-___,'_ ..... I

" _ _ - -_- - - - - - - + " - - _ - / :._-_,,' ,__ I

--'---b---'-- i / i i I

i t, i i i d

'_: :-oi i _._._, i oie -'vl

' I ': I I 1 , ,_ I , /

...._ ,.. _ .......... .+...._-:.__....-__..-.-__-.=._....,...,-,--__.__

i_____, ' ' _ ' , .-''_-,_--+--='" -4" •

1,/;_1-__;'t,.'-_\ _ -.--.-=--_-=- . .

. _L.... L -------'f----, .,._n.._,_-_,__-___....L....:..

_,',_,';j,. :_'-_'_._'.t".,',.,." : ,'

(d)

/If / _ --I.l_-- _!% %%T\11i ._.

/',','7"_.;-"._--- ':---_-,:;-'.',.'.',','___JJ) /

I _ '_ • - - ...... _ :- " - _" • • \ _._k._'"._.--._. --JJZJ I

1-"...... _-- -.i_-_---_---i._-_ -_-- -"_----=-_--_----'- ...... T-----,'- .... _- '_ - ---_.._--'_--:.2 -- _......

/ -..:_ --t----_-'_--_-" '.- _ :' ,--¢'-'-___.----kS'

: I I I I I I I I I

'. ', l : I _ I I l I I

................ .....,....,.__,.--..--.-,......--,_..,..-- .'...>........

I//__f/f_,____k_k_/I,,,2, 't,.',:., .-_-.. "_,'' i.J,.'.',-'.,-.I

r1V__-T-TTFTT-IT; i: -c-'CJ-C_",,-Tt ..... ::_--:: ""]-_"7""rTT';'.'*'!-r'l

Figure 2 (continued)

33
(e)

Ii ] _ -7-.., -._-=-2---L_k ,. I

"r_-r'_"i--Lr----_ _I_-..... -,-------_- ...... _---,----r----I--,--,----L-

...... _ -- - - _ - _- - --':'-- --I.---- _---_.-=:---- --_-_ -+,----'-...... .'I'--:; .......

I i t i

i

I I I l I /

i __----4---L__ i

......... _- +. ','.. _<._-....,-.4_. _. ,.+.._. ,.__..-....,,.. _._ -..._...,,,...._.-.. _._._....

I I t j. t Is t--t- t _ _ k ) l

"'r'r'T'k'_i"C _'"CTk".... z'_-'= ...... 3"'7"'7-7"i"+'+'_"r'1

I k _ _+._+ .,, _. 4 _- --l- - "_ k .. ,, sl l l l l [

I

- _ --------L=--------'=-_------_- ----= _ - " /

Figure 2 (continued)

34
(a)

(b)

............ !............ i----_ ........... t----_-i ......... •.......... i.......... _'_-,,:: .........

Figure 3: (a) Geopotential at 50 days for the integration with the initial condition of McDonald and Bates (1989) and a value of 0.03 for the uncentering parameter(contours every 60.0

meters); (b) Same as (a) but with the fourth order diffusion applied to geopotential(contours

every 60.0 meters); (c) Same as (a) but with the fourth order diffusion applied to divergence

(contours every 60.0 meters); (d) Same as (a) but with the fourth order diffusion applied

to both geopotential and divergence(contours every 60.0 meters).

35
(c)

(d)

Figure 3 (continued)

36
(a)

,_ _." ', _ ',_ ' _ ._i_

._-'= ._ ___............. ,

.......I-.... ""

"'"_':_"::_ ',__'-_" _ _,_-21_ u'_:_'_-_ ., i_:'_°

..-_. ........ _....__._.: .......... _..-_--_-----O_-o-i---,_ ....... .i---_---,-':_-_,,. ....._..... $--_--.

Figure 4: (a) Geopotential of the ECMWF analysis of 0000 UTC. 1 December, 1992 (contours every 60.0 meters); (b) Geopotential after 3 days of integration with a value of 0.05

{'or the uncentering parameter(contours every 60.0 meters); (c) Same as (b) but with a value

of 0.015 for the uncentering parameter and the the fourth order diffusion applied to both

geopotential and divergence (contours every 60.0 meters).

37
(b)

(c)

i ' , | : , , : I _, ___ , i__.'._ _._--r.....___'_ .....

,. ' ._ -.i : I ' '

,. ' ..... [ _..

i e • o

:; I /'._'_[ _'*_¢,_ I . I I 1 ". I ...-12"I

- '_ l_ , "_ _ i i , \J/,_ _ , .....,' ,

'. ' , . ' , ,_ , , , J I ;" I

................ t....... "_'_----_, ,'I-.--._.... _ ........ I.......... t......... -'- ..... " ............

Figure 4 (continued)

38
_N"

60N

(a)

EQ

30S

60$

gO$,

0 60E 120[ 180 12"0W 60W 0

_K3N

GON.

(b)

EQ.

30S.

605

90S

0 601[ 120 r 180 t20W 60W 0

Figure 5:50 mb temperature field from NASA/GLA 3-D semi-Lagrangian semi-implicit

GCM without the application of the horizontal diffusion. The value for the off-centering

parameter is 0.05. (a) 3-day forecast; (b) 5-day forecast; (c) 7-day forecast; (d) 10-day

forecast. (Contour interval is 3 degrees).

39
90N

90S

0 60( 120( 180 12011/

90N

Figure 5 (continued)

T

12or

4O

(c)
.207_

(a)

60S

30S'

50S

Figure 6:50 mb temperature field from NASA/GLA 3-D semi-Lagrangian semi-implicit

GCM with the application of the horizontal diffusion. The value for the off-centering parameter is 0.05. (a) 3-day forecast; (b) 5-day forecast; (c) 7-day forecast; (d) 10-day forecast.

(Contour interval is 3 degrees).

41
30N

(0

.210

gOS

Figure 6 (continued)

42

(c)

(d)
Previous Volumes in this Series

Volume 1

September 1994

Documentation of the Goddard Earth Observing System

(GEOS) general circulation model - Version 1

L.L. Takacs, A. Molod, and T. Wang


REPORT DOCUMENTATION PAGE For., Approve_

OMB No. 0704-0188

Public reportingburden for this collection of information is estimated to average 1 hour per response, includingthe timefor reviewing instructions, searching existing data sources.

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

October 1994 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

Technical Report Series on Global Modeling and Data Assimilation

Volume 2 - Direct Solution of the Implicit Formulation of Fourth

Order Horizontal Diffusion for Gridpoint Models on the Sphere

6. AUTHOR(S)

Yong Li, S. Moorthi, and J. Ray Bates

Max J. Suarez, Editor

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS (ES)

Laboratory for Atmospheres

Data Assimilation Office

Goddard Space Flight Center

Greenbelt, Maryland 20771

9. SPONSORING / MONITORING ADGENCY NAME(S) AND ADDRESS (ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

C - NAS5-32332

Code 910

8. PEFORMING ORGANIZATION

REPORT NUMBER

94B00143

10. SPONSORING / MONITORING

ADGENCY REPORT NUMBER

TM- 104606, Vol. 2

11. SUPPLEMENTARY NOTES

Y. Li: General Sciences Corporation, Laurel, Maryland (at GSFC Data Assimilation Office,

Greenbelt, Maryland); S. Moorthi, National Meteorological Center, Camp Springs, Maryland;

M. Suarez and J. Bates: Goddard Space Flight Center I Greenbelt1 Mar_'land

12a. DISTRIBUTION I AVAILABILITY STATMENT 12b. DISTRIBUTION CODE

Unclassified - Unlimited

Subject Category 46

This publication is available from the NASA Center for AeroSpace

Information, 800 Elkridge Landing Road, Linthicum Heights, MD

21090-2934_ (301)621-0390.

13. ABSTRACT (Maximum 200 words)

High order horizontal diffusion of the form KV 2 is widely used in spectral models as a means of

preventing energy accumulation at the shortest resolved scales. In the spectral context, an implicit

formation of such diffusion is trivial to implement. The present note describes an efficient method

of implementing implicit high order diffusion in global finite difference models. The method expresses the high order diffusion equation as a sequence of equations involving V 2. The solution is

obtained by combining fast Fourier transforms in longitude with a finite difference solver for the

second order ordinary differential equation in latitude.

The implicit diffusion routine is suitable for use in any finite difference global model that uses a

regular latitude/longitude grid. The absence of a restriction on the timestep makes it particularly

suitable for use in semi-Lagrangian models. The scale selectivity of the high order diffusion gives it

an advantage over the uncentering method that has been used to control computational noise in twotime-level semi-Lagrangian models.

14. SUBJECT TERMS

GEOS, General Circulation Model, GCM, DAO, Semi-Lagrangian, Implicit,

Diffusion

17. SECURITY CLASSIRCATION 18. SECURITY CLASSlRCATION 19. SECURITY CLASSlRCATION

OF REPORT OF THIS PAGE OF ABSTRACT

Unclassified Unclassified Unclassified

NSN 7540-01-280-5500

15. NUMBER OF PAGES

54

16. PRICE CODE

20. UMITATION OF ABSTRACT

UL

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39.18