NASA Technical Memorandum 104606, Vol. 5

Technical Report Series on

Global Modeling and Data Assimilation

Max J. Suarez, Editor

Goddard Space Flight Center

Greenbelt, Maryland

Volume 5

Documentation of the

ARIES/GEOS Dynamical Core:

Version 2

Max J. Suarez

Goddard Space Flight Center

Greenbelt, Maryland

Lawrence L. Takacs

General Sciences Corporation

Laurel, Maryland

National Aeronautics and

Space Administration

Goddard Space Flight Center

Greenbelt, Maryland

1995
This publication is available from the NASA Center for AeroSpace Information,

800 Elkridge Landing Road, Linthicum Heights, MD 21090-2934, (301) 621-0390.
Abstract

A detailed description of the numerical formulation of Version 2 of the ARIES/GEOS

"dynamical core" is presented. This code is a nearly "plug-compatible" dynamics for use

in atmospheric general circulation models (GCMs). It is a finite difference model on a

staggered latittMe-longitttde C-grid. It uses second-order differences for all terms except

the advection of vorticity by the rotational part of the flow, which is done at fourthorder accuracy. This dynamical core is currently being used in the climate (ARIES)

and data assimilation (GEOS) GCMs at Goddard.

I? Sl NOT FU.Mr 

.,°

111

d.

• - • L. :, %4 .,_

Contents

List of Figures

1 Introduction

vii

1

2 Continuous Equations 3

3 Conservation Properties

3.1 Global Mass Integrals ..............................

3.2 Total Energy ...................................

3.3 The VerticMly Integrated Pressure Gradient Force ..............

4 Horizontal Coordinate Systems

4.1 Latitude-Longitude Coordinates .........................

4.2 Equations in Component Form .........................

4.3 Coordinate Rotations ..............................

5

5

5

7

7

7

8

9

5 Vertical Differencing 11

5.1 The Vertical Grid ................................. 11

5.2 The Discrete Equations ............................. 11

5.3 Energy Conservation and the Discrete w .................... 15

5.4 The Vertically Integrated Pressure Gradient Force .............. 17

6 Horizontal Differencing

6.1 The Horizontal Grid ...............................

6.2 The Continuity Equation ............................

6.3 The Thermodynamic Equation .........................

18

18

21

21

:,PREC,II_I)II_W_ PA IDL,ANK i_K)T PIL,SINED PAGE -\_:_ INTENI_C_ALLYBLN_K
6.4 The Tracer Equations ..............................

6.5 The Momentum Equation ............................

6.5.1

6.5.2

6.5.3

6.5.4

The discrete equations in component form ...............

Fourth-order vorticity advection: the choice of c_,fl, ...........

Energy conservation and the momentum equation ..........

Instabilities of the discretized vector-invariant form ..........

22

23

23

25

28

32

7 Polar Filters 35

8 Time Differencing 36

8.1 The Brown-Campana Scheme .......................... 37

Using Subroutine DYCORE 38

9.1 The Argument List ................................ 38

9.2 Comments ..................................... 39

9

10 Final Remarks 4O

vi
List of Figures

The geographical and computational coordinate axes .............

The vertical grid .................................

Position of the variables on the C-grid and the indexing convention .....

Stencils of _s used to define the c,, fl, etc ....................

Stencil showing the indexing of the c_, fl, % 5 factors, which appear in the

energy conserving form of the inertial terms ...................

Stencil showing the indexing of the u and # factors, which appear in the

energy conserving form of the inertial terms ...................

10

12

19

28

31

32

vii

1 Introduction

Tile need to use modular coding techniques in atmospheric general circulation models

(GCMs) has been recognized for some time. The two main benefits of this approach are

that it makes codes much easier to maintain and modify, particularly when the work must

be done by seve.ral people, and that it makes it easier to exchange and compare codes. Modular coding practices are thus being used to solne extent by most atnlospheric modelers; but

most of the effort has gone into developing modular (:odes for the physical parameterizalions, such its radiation and cumulus conw_ction. Kalnay et al (1989) proposed a set of rules

for coding physical parameterizations aimed at facilitating the exchange of codes. If these

"phlg-compatible" rules are followed in coding both the GCMs and the parameterizations,

codes can l)e "unphlgged" from one model and "plugged" into another with little effort.

In this report we introduce the notion of a "dynamical core," which attempts to extend the

ide.as of t_lug-conlt_atible parameterizati_ms to the coding of the dynamics. Since the maimer

ill which the eq, lations of motion are discretized affects so flmdamentally the organization

of both tile code and the data in a GCM, it is not practical t,_ make the dynamics codes

truly phlg-conlpatible. They can, however, be made sutIiciently modular that once the

GCM's data structure is inodified to accotnmodate them the computational routines can

be e`asily ported.

The dynamical core we will be describing consists of a set of subroutines that compute

the time tendencies of winds, temperatures, surface pressure, and an arbitrary number of

tracers. The core is invoked by calling a single subroutine, and all data is passed to the

core through that subroutine's argument list. Inputs to tile core include: the dilnensions of

the grid, an array defining the vertical distribution of the layers, some physical parameters,

the state, variables at two time levels, the time step, and the time tendencies of each state

variable due to other processes. The core updates tile time tedencies to include the effects

of the dynamics. It does not update the state variables, or perform any temporal or spatial

damping that may be needed to control non-linear conlputational instability.

Ideally one would want the core to be independent of the time differencing, but this is

not possible except for the simplest explicit differencing. By passing the time step, two

time levels of state variables, and tendencies due to other processes, the core can compute

tendencies for other time differencing methods, such ,as the semi-implicit scheme commonly

used in spectral models. The version we describe here is for explicit time differencing only,

but it can produce "economical explicit" tendencies (Brown and Campana 1978) for use

with a leap-frog scheme.

All state variables and tendencies in the argument list are on the core's horizontal and

vertical grids in our case, an unstaggered or Lorenz grid in standard sigma coordinates in

the vertical and an Arakawa C-grid in the horizontal. Although the type of discretization

is determined by the dynamical core and the GCM's data structure must conform to it--

other things, such ,as the resolution or the placement of the layers in the vertical, need not

be. In our dynamical core, these are controlled through the argument list and are thus
completely arbitrary. (The resolution can even be altered from one call to another within

a single program.) Thus, the core can bc adapted to many applications by simply altering

the call, rather than modifying the source and rebuilding the binaries. This should simplify

the use of the core for those not interested in tinkering with the source code. Oil machines

for which this dynamic allocation results in an inefficient code, this feature can be cmsily

discarded and the array sizes fixed at compile time.

The differencing used in the ARIES/GEOS core follows closely that which h,'Ls been in use

in the UCLA GCM since around 1980. Both the ARIES/GEOS core and the UCLA model

arc ba,sed on a C-grid, they use the same spherical grid, they both use Arakawa and Suarez

(1983) vertical differem:ing, and they use similar forms of the momentum equation. There

are, however, a number of important differences. The UCLA model, for example, use.s a

modified sigma coordinate, while the ARIES/GEOS core uses standard sigma. Rather than

detailing these differences, we will present a complete description of the ARIES/GEOS core.

This doemnentation is for Version 2 of the core. Version 1 was very similar, differing

primarily in that it used a second-order scheme for advection in the momentum equation.

Version 1 hms been used for a number of applications at Goddard. In particular, it is the

dynamics of the GEOS-1 GCM (Takacs et al. 1994), which w,_s used for the re-anMysis of

the 1985-1990 period recently completed at the Data Assimilation Office (Schubert ct al.

1993).

Two additional features are new in this version of the core. First, the code has been

generalized so that the latitude-longitude coordinates on which the horizontal grid is based

need not coincide with geographical latitude and longitude. By displacing the pole of

the computational coordinate from that of the geographical coordinate, the more unifi)rm

grid near the computational equator can bc brought to other regions of interest. This

device h,'_s already been used successfiflly to model the polar stratosphere. Second, the map

factors on the grid have been generalized to allow grids that are not uniform in latitude

or longitude. Latitudinal stretching of the grid, with or without pole rotation, is already

implemented; longitudinal stretching requires modifying the polar filters, a change that will

not be available until the next version. The hope is that with these modifications w_,,can

easily convert the global models into regional, high-resolution models.

The presentation follows the general form and notation of Arakawa and Lamb (1977): repeating their discussion of the continuous equations in sections 2 and 3. Latitude-longitude

coordinates and rotations between the computational and geographical coordinates are discussed in section 4. The vertical differencing is presented in section 5 and the horizontal

differencing, including a new fourth-order scheme for advection of vorticity is presented in

section 6. Section 7 describes the polar filters and section 8 the interaction of the dynmnical

core with the time differencing. Section 9 describes the FORTRAN interface.
2 Continuous Equations

We use a a coordinate defined by

P --PT

- , (1)

?r

where p is the pressure, 7r = PS -PT, PS is the surface pressure, and PT is a constant

prescribed pressure at the top of the model atmosphere. With PT = 0 this coordinate

reduces to the conventional cl coordinate proposed by Phillips (1957). The code is written

in this more general a coordinate for backward compatibility with earlier wwsions of our

GCMs; in the current versions, however, we take PT = 0.

With this vertical coordinate, the continuity equation becomes

o_ _ -v.. (_v) o(_a), (2)

Ot Oct

where v is the horizontal velocity vector. Integrating (2) and applying the boundary conditions _ = 0 at p = PT and p = PS, we obtain the forms used in tim modeh

07r Ot _ j_01 V_ - (Trv) dc_ (3)

and

(_r/7) = -a_ - V_- (7_v) dc_. (4)

The equation of state for an ideal gas is c_ = RT/p, where a is the specific volume, T is

the temperature, and R is the gas constant. The following alternative forms will be used

below:

" N- N (5)

where 0 = TIP is the potential temperature, P = (p/po) '_, t,_= R/cp, % is the specific heat

at, constant pressure, and P0 is a reference pressure. In obtaining the forms in (5) we have

used 0g = x_p mad the relation

(o2) (o2) 7r (7"

(6)

For the time being we are ignoring virtual effects.

The hydrostatic equation is

0q)

0p
where _5 is the geopotential. Using (5) this can be written:

(02)= -rra = -Tr q's

rr (7 a 7r

From (7) we obtain

7r

which is the form used by Arakawa and Suarez (1983) for vertical diseretization.

(7)

(8)

The thermodynamic equation is written in flux form to facilitate the derivation of a 0-

conserving differencing scheme:

0(_o) o(_ao) _Q

ot - v,.(_v0) 0_ +q,/,, (9)

where Q is the diabatic heating pet" unit ,nass.

In addition to the equations of nlotion, the core. computes tendencies for an arbitrary number

of atmospheric constituents: such as water vapor and ozone. These are also written in flllX

ft H'In:

O(rrq (_:))

- O(_rirq(_:)) _r$ (_,), (10) v_. (_vq(*')) +

Ot Oa '

where q(_') is the specific mass of the kth constituent, and $(_') is its source per unit mass

of air.

The momentum equation is written in "vector-invariant" form, as in Sadourny (1975)

and Arakawa and Lamb (1981), to facilitate the derivation of an energy- and enstrophyconserving differencing scheme:

Ov Ov .q OT

0t - -(f + () k x v - ?*0a V.(q_ + K) - q,O V_P - ---,TrOc* (11)

= -,lkx (Trv)-/_ Ov V,(_5+K)-%0 (Op) _ Vrr---- o07- (12)

a 7r GQO _

where f is the Coriolis parameter, k is the unit vector in the vertical,

(=(V, xv).k

is the vertical component of the wn'tieity along a surfaces,

is an external potential vorticity,

(.f + ()

7r

1

K = _(v. v)

is the kinetic energy pet" unit mass, g is the acceleration of gravity, and T is the horizontal

frictional stress vector.

4
The advective form of the momentunl equation is obtained by using the identity

1

_V(v-v) = (v.V)v - (v x v) x v

in (11):

Ov _ .f k x v - (v. V_)v - d_0v V_ - Q,O (0/') _ V_r .qOT . (13) Ot _ 7r Oa

3 Conservation Properties

3.1 Global Mass Integrals

Integrating (2) over the entire atmosphere we obtain a statement of conservation of total

atmospheric m_ss,

0; _ds = O,

where S represents the surface of the sphere. Similarly integrating (9) and (10), we obtain

equations fi)r the global mass-weighted integrals of 0 and q(a'):

and

0 1 1

3.2 Total Energy

To discuss conservation of total energy, we first write an equation for the temperature.

Multiplying (9) by cvP we obtain, after some manipulation:

O(uc,,T)

Ot -- -V_. (rrvcvT) O(_rcpT)oa

+_revO -_-+v.V_P+_ _ _ +rrQ. (14)

The first two terms on the right-hand side of (14) vanish when integrated over the entire

domain. The first term on the second line represents the conversion to kinetic energy. Using

(6), this term can be written ,as:

OP VTr) +_r&] (15)
The thermodynamic equation can then be written in the more conventional form:

O(_rcpT)

Ot - V_. (_rv%T) O(_r_cpT) 0a + (_) + _Q' (16)

where w is the vertical p velocity:

w=lS=r6+o (°i+v'Vr ), (17)

and a is given by (5).

Next we obtain the kinetic energy equation by dotting Qrv) with (12). The first term on

the right-hand side of (12) is an "inertial force," which does no work since it acts normal

to the motion. The remaining terms give:

O(TrK) -V,. (TrvK) O(_6K) (07r O(r&) Ot - Oa +K.-_-+V_-(Trv)+ Oa ]

op 07 (18)

The last term on the first line of the right-hand side vanishes from conservation of mmss

(2). Terms on the second line represent the work done by the pressure-gradient force the

conversion from potential energy. We next consider the first of these, -(Trv• V_O):

-(_v. V_) = -V_-(_v_) - • (07r -5/+ 0(Tr&)) 0----g-

= -V_-(_) 0(_)0o _0_ + _°b-g_" 0_

=-vo. (_v_)o(_)o_ o(_) O_o_o_+_°_(_ +_°r)

where we have used (2), (7), and the identity

o(_e) 0¢ (20)

- Oa a Oa"

The kinetic energy equation (18) can now be rewritten ms

O(_K) V_. (_v_) Ot - -V_. QrvK) O(_riTK)oa O(rh_)Oa

o,oo, Oc_ Ot _r--a -_-_ _ [7r/r+a -_-+v. VTr)]-gv. Oa" (21)

Finally, using (17) and (5), we obtain a kinetic energy equation in which the conversion

term is written in the same form as in (16):

O(TrK) = -V_. (TvvK) O(Tr#K) O(_ri_¢) Ot 0_ V_. (_rv_) - - Oa

07"

O(aO) O_r (Trwoz) - gv . . (22)

Oa Ot Oa
Summing (16) and (22) and integrating over the entire domain, we obtain an equation for

the rate of change of the total energy of the atmosphere:

0 ! l

g

where _s is the surface geopotential.

--V • --d,s_

7r .q

(23)

3.3 The Vertically Integrated Pressure Gradient Force

The vertical integral of the pressure gradient force is irrotational in the absence of spatial

variations of the surface geopotential. In pressure coordinates,

Ps Vl,_dp = V c_dp - 62sVps ,

J PT L_ PT

= V (#P - _s) dp + (Ps - PT) VgPS,

LJ PT

In the last two forms, the curl of the first term is zero and the second term vanishes if (I)s

is constant.

To satisfy this property in the discrete model, we will need to obtain this form starting from

the pressure gradient force written in terms of gradients along a surfaces.

1 l 0P

f0 7rVpq)da = if0 7r[V,,g2+cpO(--_),VTr] dcr,

_- V[fo'(Tr,i_)da]_V_rfo' [_ +a-_a]Oa_ da,

= V[fo'(_rq))da ] -(I)sVTr,

= V [/o' _r(_- _s) da ]+ _rV_s, (25)

where we have used (7) and (20).

4 Horizontal Coordinate Systems

4.1 Latitude-Longitude Coordinates

In section 2, the continuous equations were written in vector notation. In this section we

rewrite thc equations as used in the codc, in component notation in spherical coordinates.
Let ¢ be the latitude and A the longitude. The differential elements of arc length along

latitude circles and meridians are:

dx = a cos ¢ dh (26)

d.v = a d¢, (27)

where a is the radius of the globe. The differential element of surface area on the sphere is:

The gradient of a scalar field, ¢, is

ds = a 2 cos ¢ dh de. (28)

0¢ _0¢

V¢ - + - -- (29)

a cos ¢ Oh a 0¢'

where A and ¢ are the horizontal unit vectors in longitude and latitude.

4.2 Equations in Component Form

Let u and v denote the zonal and meridional components of the horizontal velocity v in

the (h, ¢) coordinate system. To write the scalar conservation equations for mass, 0, and

constituent concentrations in component form, we require only the form of the horizontal

divergence of their mass-weighted fluxes:

V. lrvq -- 1 (Oruq 0(cos ¢ 7rvq) ) acos¢\ 0A + 0¢ ' (30)

where q can be replaced by 1 for use in (2), by 0 for llSe in (9), and by q(_') for use in (10).

The zonal and meridional components of the momentum eqnation (11) are:

0_ - (f + Or- _0_ acos¢/ 0-X +%0 _ _-_ 7rOck'

0vOt _ 0v 1 [0(_+K) (OP) 0__¢] .007" (32) ('f + ()u - 6 0a o, [ -0¢ + Q,O -_ _ 7r Oct'

where f = 2f_ sin ¢ is the Coriolis parameter and f_ is the rotation rate of the coordinate

system relative to the inertial frame. The relative vorticity, _, is

1 (0v 0(cosCu)) (33) (=V_xV-acos¢ Oh 0¢

The quantity (f + () is regarded ms the absolute vorticity, since .f is the curl of solid-body

rotation at rate l-h

f = v x (_a cos ¢_) -

1 O(f_a cos2¢)

a cos ¢ 0¢ -- 2f_ sin ¢. (34)
Eqs(31)and(32)canbewritten in themoreconventional advective form by combining the

relative vorticity and K terms:

1 OK v {0v 0(cos (bu) \

£'U

a cos ¢ OA a cos ¢ \ OA 0¢ ]

1 / Ou Ov "_

acos4,t?-5-i+ 7

,, O(cosCu) u &,,

a cos ¢ 04' a cos ¢ OA

tan ¢ u Ou v Ou

-- II'U

a ,zcos¢OA aO¢

-(11. l OK _ 1,, (Or O(cos 4)11,)) _ _1 (u Ou. + v Ov'_ a O¢ acos¢ OA 0¢ a -_ 0¢]

_ tan 4)u2 u 0v v Ov

a a cos ¢ OA a 04'

4.3 Coordinate Rotations

In the code, we allow the spherical coordinate system used for the computations to be displaced relative to the geographical latitude-longitude coordinates. Let (A, 05) denote longitude and latitude on the computational coordinate system and let (_, q_)denote geographical

longitude and latitude. The relation between the two coordinate systems is fiflly determined

by specifying the coordinates of the geographical north pole in the computational system,

which we denote by (ANp, CNP), and by a third parameter, A0, which represents a rotation

about the geographical pole. The relation between geographical and computational axes is

shown schematically in figure 1. Using these parameters, the two systems are related by:

or

sin 4_

cos(,i+ _0)cos$

sin(£ + ;_0)cos

= cos CNp cos 05cos()_ -- /\Np) + sin CNp sin ¢, (35a)

= sinCspcos¢cos(A-- ASp) - cos Cse sin ¢, (35b)

= cos¢sin(A- ANp), (35c)

sin 4) = cos 05NP cos _ COS(_ nt- )_0 -- 71") + sin 4)Np sin ¢, (36a)

cos(A+_r-ANp)COS¢ = sinCNpcos¢cos(_+A0--Tr)--cosCNpsin¢, (36b)

sin(A + 7r - ANp)cos¢ = eosCsin(_ + A0 -- 7r). (36c)

Note that these transformations can be inverted by making the exchanges:

(A,¢) _ (_,(_), (37a)

ANp ¢==> 7r-- A0. (37b)

The code assumes that all input variables are on the computational coordinate system, and

that the input u and v wind components are also along the computational coordinates. It

makes the same assumptions for the output variables. With these assumptions, the only

9
_NP

x

Y

Figure 1: The geographical and computational coordinate axes.

effect on the code of displacing the geographical pole away from the computational pole is

on the form of the Coriolis parameter:

f = 2f_sin¢ = 2_ (cosCNp cos¢cos(A -- ASp) + sinCNP sine). (38)

Although there is no other reference to quantities in geographical coordinates within the

dynamical core, it is generally necessary to transform between the two grids when initializing

a run or recording results in geographical coordinates. In addition to interpolation, this

requires transforming the velocity components between the two systems• If we let X denote

the local angle between the computational and geographical coordinates, the two sets of

velocity components are related by

u = cosx_,- sinxv, (39a)

v = sinxu. + cosxv, (39b)

and

fi, = cosxu + sinxv, (40a)

10
= cos X v - sin X u. (40b)

The angle X may be obtained from the relations:

cos ¢ cos X = cos ¢ sin ¢NP -- cos(A -- )_NP) cos CNP sin ¢, (41a)

cosCsinx = cosCNp sin(A -- )tNp ). (41b)

A subroutine package that implements these transformations is being developed and will

be documented in a future report in the series.

5 Vertical Differencing

5.1 The Vertical Grid

We use a Lorenz grid in the vertical, with both winds and temperatures defined at the

same levels. The atmosphere between a = 0 and a = 1 is divided into LM layers by

LM - 1 surfaces of specified constant al, l = 2, LM, so that al = 0 corresponds to the top

of the model atmosphere, and aLM+I = 1 corresponds to the lower boundary. At the LM

layers we define the velocity, vl, the potential temperature, 01, and the specific masses of

all *_race constituents, ql k). The vertical a-velocity, &t, is defined at the interfaces between

the layers and at the top and bottom surface. This arrangement of the variables is shown

schematically in figure 2.

Note that both layers and interface levels are numbered from top to bottom and that, to

maintain as close a correspondence with the code as possible, both are indexed by integers,

rather than using the more common convention of assigning the interfaces half-integer values. Level al thus corresponds to the top of layer I. Where there might be some confusion as

to whether a quantity is defined at the layers or at the interfaces, we will mark the interface

quantities with a hat (e.g., /5). Vertical differences will thus be written ms

(6h)l = ht+l - ht

at the layers, and A

( a)t = al - al-1

at the interfaces.

5.2 The Discrete Equations

Using (1), we define the pressures at the top, bottom, and interface levels as:

Pl = PT A- _ral, for l = 1, LM + 1.

1The forms of the equations used in the code will be denoted by boldface numbers.

(42)

11
iin_fTfiI_fflfTm

u, v, O, q

2 u, v, O, q

3 U, v, O, q

LM-2 U,V, O, q

LM-I u, v, O, q

LM U, v, O, q

Figure 2: The vertical grid

LM-2

LM-I

LM

LM+!

The values of the Exner function at the interfaces are obtained directly from (42):

which can be written

Pl = (Pl_ _ , for l= 1,LM + 1, (43)

\p0 /

for PT = 0. The latter form is much more economical to use since the first factor on the

right-hand side is independent of l and the second factor, 75t, depends only on a. The

code takes advantage of such economies when PT = 0 and uses the general form only when

PT 7L 0. We will thus write the PT = 0 forms when appropriate.

Layer temperatures Tl are defined from the potential temperatures a.s Tt = OIPl, where Pl,

12
theExnerfimction at thelayers,is givenby

1 [a(p/5)] , for I= 1, LM. (44)

P/- l+g [ _Sp Jl

This is the form proposed by Phillips (1974). For PT = 0 this becomes

--lr'_ [a(a_+l)] =Tr_Pt, (44a) /'t p_ (1 + _,) _------g-- t

where, again, 7_t depends only on a.

The vertical finite differencing of the thermodynamic and hydrostatic equations and of the

pressure gradient force follows Arakawa and Suarez (1983). In their scheme, the potential

temperature "interpolated" to the interfaces is

01+1 _--- (P/+I -- [-_l+l)Ol+l + (/_/+1 -- Pl)Ol for 1 = 1, LM - 1. (45)

PI+ I - Pl

As we show below, to conserve energy this form has to be used in the hydrostatic equation

for the thickness between layers and in the vertical advection term in the thermodynamic

equation. When PT = 0, (45) becomes:

"Pt+I -- _[_/-I-I _)/+1 - 7)I Or, (46)

0/+1 -- "_--//+; 7._g 0`+1 + _'+1 g

so that the interpolation depends only on c,, not pressure.

The hydrostatic equation is

¢I)LM ---- _S + CpOLM(&M+I -- PLM), (47)

Ot = ¢l+1 + c,01+l(Pt+l -/_),

= ¢t+1 + cp[(Pt+l -/_l+l)0t+l + (/_t+l -- Pt)Ot], for l = 1, LM - 1. (48)

For PT = 0 we write this as:

C K ^

gPI,M = (I)s + .pOLM 7r (PLM+I -- "PLM), (47a)

(I)l = (I)1+ 1 + CpTr_[(P/+l - 75t+1)01+1 + (75t+1 -- Pl)Ot]. (48a)

The continuity equation is

c%r (6(r&)_ (49)

Ot- V.(_rv_)-\ _a ]l'

The equation for the tendency of _r is obtained by summing (49) over all layers:

07r LM

- E v. (_v_) (_)_, (50) Ot

t=l

13
and the verticala-velocity at the interfaces is given by

oq r l

(r&)l+l = --al+l-0-_- -- Y_. V. (rve_) (Sa)f., fro" l = 1, LM - 1. (51)

_=1

The thermodynamic equation is

(5(r&{})'_ nQt for I = 1, LM, (52)

Ot t

where/_l is given by (45) and Pt by (44). We note that using (45) as the interpolation of 0 for

vertical advection gives up conservation of 02 or some other fllnction of 0. This is the price

Arakawa-Suarez pays to use a local hydrostatic equation. The scheme of Arakawa (1972)

satisfies all conservation properties of Arakawa-Suarez plus conservation of a fllnction of 8,

but it uses a non-local hydrostatic equation, in which one of the geopotentials depends on

a weighted sum of all the temperatures in the column.

For advection of the trace constituents, we are currently using the square-conserving form:

(5(r_-q(k)) _X r,5[/¢) for I 1, LM, (53) °(rq}_)) - -v. (rv,q}k))- \ _ / + = Ot l

where

;+5 = -_ _,qt+l + q ,

for l = 1, LM - 1. (54)

The form of the momentum equation is:

OVl

Ot

1 1

= -r/t k x (rv/) [,5t(vt -- V/-1) -'t- ,_t+l(vt+l - vt)]

2 (_o)t

OPlvT"f-'q('Q") for l = 1. LM, - V(¢9 + K)t - cpOt--_ r _ t'

where

OPt = at+,(Pl+l - Pl) + c*l(Pt - Pi) for l = 1, LM,

Or r(_)l '

,-_smay be verified by differentiating (44). With PT = 0, (56) becomes:

OPt _ _t+, (_5_+, _ Pl) + _(Pt - P_)

0_ r (Sa)l

where, again, the second factor depends only on or.

(55)

(56)

(56a)

14
5.3 Energy Conservation and the Discrete ._

The derivationandjustificationof theverticaldifferencingschemepresentedaboveisgiven

in ArakawaandSuarez(1983).Herewewill merelyshowthat thesedifferential-difference

equationsareenergyconserving. We will also obtain a discrete form of w fl'om the energy

conversion term.

We start, ms in section 2, by writing the thermodynamic equation in terms of the temperature at the layers, Tl = PlOl. Multiplying (52) by cpPl we have:

0(_c,,_) { _(_c,, 0))

= -V.(,vtc,,Tt)- \ G ]_

+ cpOtTr [ Ot + vt. VPI + 7rQl, (57a)

-[- _ [ (7r_)/+1(20/+1 -- PlOl+,)- (Tr[-Y)l(r l - I-)lOl)]

- u_r -_ +vt'VTr , +Tr@ for I = 1, LM. (57t))

Here the quantity 2b hm_ been introduced to put the equation in flux form; this allows llS

to separate the enO.rgy conversion term. Comparing (57l)) with (16), we see that the boxed

terms of (57) correspond to (_rc0c_)t.

To form the kinetic energy equation we take the dot product of (55) and _rvl:

01rKl

Ot - V.(TrvtKt)+Kl +V.(Trvt)+ \ 5a ]lJ

1 1

2 (t_O')l [(71"_)/+IVl+l 'Vl -- (TrO)/VI. VI_I]

- V. (_rvt_;) + _;V. (_vt) - c.ot_r _ OPt O--Tvl. VTr- .qv,.('7-)

1

1 1

-V" (TrVlKt) 2 (Sa)t [(TrS")t+lVt+l "vt - (Trb)tv/. vt-l]

- V. (rv_t) - _ -5[ + \-UJ-a ]tJ _p tr-b-_vl. VTr- .qvl. _ t'

(58)

Now introducing the edge-level geopotentials, _, and using the identities:

(_(O'_))_ __ (Yl+l(_)/+l -- (_l ) -t- (Yl( (_l -- (_l )

0_= \--TT-_j_ (_)_ '

(59)

15
whichis a finite differencedanalogof (20),and

wecanwrite:

(Tt'O)/+l(_l+l -- ¢_1) + (Trh)l(Ol -- _l)

(Sa)t

1 1

-- -V. (Trv/Kl) 2 (t_dr)l [(Trh)l+lV/+l " vl -- (Trh)/v/. Vt_l]

j _

(_rh)_+l(,i,_+_- o_)+ (_h)l(Ot - ¢_)

+

(6o)t

, (60)

+

_t+l('St+_ - ¢'1) + ot(,:I,t- ,51)O_-

(5a)l Ot

,_ OPl

- %_11rO---_-vl• VTr, for l = 1, LM, (61)

where we have also used the continuity equation (49).

Comparing (61) with (22), we see that the boxed terms in (61) correspond to the conversion

term, -(Trwa)t. Energy conservation requires that the forms of (Trwa)l in (57) and (61) be

identicM. Comparing quantities in the first box in each equation and noting that hi =

hI, M+l = 0, we see that we must require that

Cl,(T/+I -- PlOt+l) = Ol -- ¢/+1, for l = 1, LM - 1, (62a)

%(_-P_0_) = Ol-¢t, forl=2, LM. (62b)

Similarly, considering the second boxes in e_h equation we must require that

at+l(_t+l--_St)+al(Ol--¢Pt) _ OPt for/=l, LM (63)

(_)t = -c,0lTr 0--;-' "

Using (56), we see that this l_st is satisfied if the half thicknesses are given by the following

hydrostatic relations:

(_1+1 - lit) = -cvOt(Pt+l -/}), for l = 1, LM, (64a)

(Or - ¢_) = -cvOt(Pt - Pt), for l = 2, LM. (64b)

It is straightforward to verify that these forms of the half thicknesses are consistent with

(48), the hydrostatic relations between layers. Furthermore, substituting (64a,b) in the

2(LM - 1) relations (62a,b), the latter can be written:

cp(Pl+l -- Pl)Ol+l = _l -- _l+l, for l = l, LM - 1,

which are identical to (48), and

cprl= __ (Pl-[-Pl-l)Ol+-_(_l+Ol-l)--¢l, 1 for 1 = 2, LM, (65)

16
which define the LM - 1 interface temperatures we introduced in writing (57).

Using the half thicknesses (64a,b) we can rewrite (61) as:

07rKl 1 1

Ot - v. (_vtgt) 2 (_)t [(_ah+,vt+, -w- (_a)lvt. w-z]

-V.(_vtCt)- _ J_-\ _ ]_N-'qvt" _ t

(66)

where the boxed term again corresponds to -(rwa)t. Using (56), this term can be written

`as:

cO OPl: (Tr0")/+,(/_/+l - P/) -I- (Tr_)/(/_- Pl) (07r )] (_)t = "P_-g-#-_ _t+_(Pt+l - Pt) + _(Pt - it) + -_ + vt. v_

= rrT-- _rO+&t +vt.VTr (67)

_t 0_ -0_ "

Here tildes denote the weighted averages:

(_, _t) = (_' _)_+_(-_t+l - Pt) + (_0, _)j(P_ - _)

P/+l- /_l ' (68)

which may be thought of ,asinterpolations of these "edge" quantities to the layers. The term

in brackets on the right-hand side of (67) is a discrete form of wt. Clearly, the corresponding

c_t is:

cpOl OPt

OxI --

5t O_r

Note that the wt obtained from (67) is defined at the layers, not at the interface levels,

where at is defined.

Multiplying (57) and (61) by (_a)t, summing over all layers, and integrating over the horizontal, we obtain the verticMly discrete form of the total energy equation:

0fsLM .LM[ (ST)](_a)tds"

/=l

5.4 The Vertically Integrated Pressure Gradient Force

The discrete form of the vertically integrated pressure gradient force is:

LM LM

/=1 1=1

OP_Vlrl

+ cpOt--_ j (,_a)l

17
LM

/=1

Using (63), the identity (59) can be rewritten:

f_l \ Y ]l

Substituting this in (69) gives:

LM

l=l

which is in the same form a.s (25).

"_- CpOlTr _l _7 7r] ( 6o')l .

= V rCt)(6a)l - _sVTr

LM ]

(69)

6 Horizontal Differencing

6.1 The Horizontal Grid

The computational domain consists of a longitude sector extending from pole to pole, with

cyclic boundary conditions applied in longitude. The size of the sector is determined by

specifying the N-fold symmetry on the sphere. If Af = 1, the domain is the full globe.

Smaller sectors (Af > 1) can be specified only if the coordinate pole and the geographical

pole coincide.

The grid divides this domain into JM latitude belts of width (A¢)j, j = 1, JM, and IM

zonal sectors of width (A,k)i, i = 1, IM. We define the grid (hi, Cj) as follows:

and

71"

¢0 = -3; Cj = ¢j-1 + (A¢)j, for j = 1,JM; CJg = 2; (70)

7r 7r

)_1 = --_; )ti+l = hi + (A,_)i, for i = 1, IM; )_IM+I = _. (71)

Obviously, the (AA)_ and (A¢)j must be specified so that:

IM 271" JM

= and (zx¢ b =

i----1 j=l

The prognostic variables are located on an Arakawa C-grid. The pressure, temperature, and

Ml tracers are located at the points (hi, Cj), for all i and for j = 1, JM - 1, which excludes

18
the poles. We will referto these,as the "p-points." The "u-points," at which the zonal

wind components are defined, are located between the "p-points" and on the same latitude

circles, while "v-points" are located between "p-points" and on the same meridians. The

vorticity is defined at the "(-points," on the same latitude circles ,as v and on the same

meridians ,as u. In addition, the two pole points are considered both v- and (-points. This

arrangement is shown schematically in figure 3, which also shows the relation of the grid

intervals to the staggered variables and the indexing convention. Note that (AA)/is defined

\

"tti _1,1

\

\

\

(,

ti,l

/

/

/

0

(zx¢b

SOUTH POLE

(Aa)_

A

,-Pi,j-- ui,j-- Pi+ I,j--

""Vi,j ....... £i,j ....... Vi + 1 ,j""

I I

-Pid-l-- u/,,j-l-- Pi+l,j-l--

i : i

INTERIOR

q_JM

/

/

/

_6/- 1 ,J M -

M

\

\

\

7t _'_1 M - I

/ Pi JM-1 \

_-I,JM-, J (i,J_M- 1

/ .............. _KJfi=i......... \

/ P,,JM 2 \

/ Ui-l,JM-2 _ti,Jg_2

I

NORTH POLE

Figure 3: Position of the variables on the C-grid and the indexing convention.

at the u-points and (A¢)j at the v-points. Variables defined at p-points are indexed the

same ,as the u-point variables half a grid interval to the east, the v-point variables half a

grid interval to the south, and the _-point variables to the southeast. With this convention,

we define the averaging operators as:

w =

( )i,j l(ui "_ ?'i--1)j

At the p-points (iTJ)i,j = 1⁄2(vj + vj+l)i

{ (p-_)_,j =

At the u-points 7J (_)_,_

l (pi -t- Pi+l)j

1

and similarly for u and v at the (-points, and p and ( at the v-points.

operators, (Siu), (6ip), etc., are defined analogously.

The difference

19
We. will designate the grid spacings in arc length in longitude and latitude as Ax and

Ay, respectively. These will be defined at the four points, p, u, v, _, and denoted by the

appropriate superscript; thus, (AUx)i,j is the longitudinal grid spacing at u-point (i, j). For

the interior points we take:

(AUx)i,j = acosCj(A_)i , forj = 1, JM - 1, (72a)

(APz)i,j = (-_-_X)i,j , for j = 1, JM- 1, (72b)

(AVx)i,j = (_-"_2)i,j , for j = 2, JM - 1, (72c)

(A"y)j = a(A¢)j, for j = 1,JM, (73a)

(APy)j = (_-_J)j, for j = 2, JM - 2, (735)

(AUy)j = (APy)j, for j = 1,JM - 1. (73c)

The intervals (A"x)i,0, (A=x)i,jM, (APx)i,0, (APz)i,JM, (AUy)o, (Auy)JM, (APy)0, (APy)JM,

(AVy)t, and (A"y)jM do not appear in the equations. The remaining intervals are special

cases at the poles:

(AVx)i,1 = 0, (74a)

A"x" (74b) )i,JM = O,

(APy)_ = a[(A¢)l + 1⁄2(A¢h], (r4c)

1

(APy)jM_ 1 = a[(A¢)jM + _(A¢)JM_I]. (74d)

The grid area elements are:

2

(Ap)_,j = (-_V_/)_,j(APy)j,

2 _AUx, ,AU (Au)i,j = ( )i,j( Y)j,

2 V V

(h.)_,_ = (A x)_,j(A Vb,

forj = 2, JM- 2,

for j = 2, JM - 2,

forj = 1,JM,

for j = 2, JM - 1,

(75a)

(755)

(ZSc)

(75d)

and at the poles

2 2

(Ap)i, 1 : (mu)i, 1

A2 2

( p)i,JM-1 = (Au)i,JM-I

(i_)l

= (APx)_,_(A"y)I,

= p A v (A _)_,JM-_( .V)JM,

1 [i_M _(Ap2 ]

: -- )i,1 ,

2 i=1 .I

IM A 2

1 1 _( p)i,JM-I (_)J" = _ i_ _=_

(75e)

(ZSf)

(75g)

(75h)

2O
6.2 The Continuity Equation

We. use a simple centered second-order form of the continuity equation:

Ozri,j _ 1

[[_iUl + (_j'Ol ]i,j -- \ [_l(y / i,j,l '

where u* and v* arc:

u_,j, t = (A"y)j(-ff'ut)i,j, for j = 1,JM - 1,

Vi:j, l = (Avx)i,j ('_3Vl)i,j, for j = 1, JM.

(76)

(77)

(r8)

Summing (76) we have:

O-t -- --_--1 (A2p)i'j[_iu£

+ [)jV_]i,j(_O')i (79)

and

Ozri,j t 1 . .

(Tr[Y)i'j'l+l-_--O'l+l Ot _1= (A2)i,j" [6iu_ t " + 5jv_]i,j(Sa)_, for I = 1, LM - 1, (80)

which are the forms used in the code. The boxed term in (79) is the dynamical tendency

of _r computed by the core. Although there should be no other contributions to the Ir

tendency, the code adds the dynamical tendency to the _r increment,. Increments of the

surface pressure (if any) passed to the core are not included in the _t that appears in

(80); however, they can enter the calculation in the case of economical explicit differencing

described below.

We note that from (74a,b), vi; l = vi*,jM = O. Using these and cyclic conditions in longitude,

the global sum of the divergence vanishes, and we have:

_0 ___,(_ta) _ E _ri,j (Ap)_,j 2 = o. (81)

Ott , j

6.3 The Thermodynamic Equation

We are using the simple second-order, square-conserving (in the horizontal) form of the

potential temperature equation:

O(_rOl)i,j_ 1 " (3l(Tr&O)_ ] (*rQ'_ (82)

Ot (A2)i,j[_i(u*_)l"l-_j(v*-O2)lli'j--_ _la ]id,l [+\cPPl/i,J'

21
wherethe0 are given by (45). The boxed quantity in (82) is the dynamical tendency added

by the core to the _r0 increment. Note that it is a mass-weighted tendency. Taking the

mass- and area-weighted sum over the domain gives:

_ ( ,,)_,_= F_,(_o)t _2

ot t _ _ _ • \ cp/3] _,j(

p)i,j. (83)

×

Multiplying (82) by (A 2 cpPl)i,j and summing over the horizontal we have:

E E [0@rcpTt) o 0t3 0_1 2

i j

i j " " I i,j

In obtaining the first term on the second line, we have used the relations:

E o,_,b= - _ b,e,a, _ _,_= o, E o,r = E b,r, (85)

i i i i i

and similar relations for j with the polar boundary conditions properly accounted for. The

second term on the second line of (84) does not involve the horizontal differencing and

so can be manipulated like the analogous term in (57); then using (628) and taking the

mass-weighted sum over all layers we have:

0 E(_,a) ' E E(rrcpTl)i,J(A_)i,J

Ot

{

= E E + Z Ag),,,

t i j l _ j

"Tr_/+l(/_/+l- /3)-]-Tr(Tl(P l -- P/)] [(dPl_ 07r ((u*(_iP)li-] - (v'_jP)lJ) -

(86)

6.4 The Tracer Equations

We. are currently using the same horizontal advection scheme for tracers as for potential

temperature advection:

O(Trq}k))i,j_ ot 1 [ 5j(v.q-_j)l] (Sl(Tro'O(k))) (TrS}k))i,j, (a_)_,_ _(u'q(_=_))_+ - + (sz)

i,j \ 61_r i,j,l

where the 0(_) are given by (54). Summing over the domain we have

-o (_,),,;= E(_,o)E E @s?)) Ot i,j ( p)i,j

l i j l i j

(88)

22
6.5 The Momentum Equation

The differencingof the momentumequationis patternedafterthat of the UCLA model.

For the vertical advectionterm and the pressuregradientforceweusea straightforward

energy-conservingform,andforthe inertialterms,afourth-orderversionof theform of the

Sadournyschemeusedby BurridgeandHaseler(1977).Toderivethis fourth-orderscheme,

wefolloweda proceduresimilar to that usedby Takanoand Wurtcle (1982)to obtain a

fourth-orderversionof the Arakawaand Lamb (1981)scheme.Tile resultingschemeis

somewhatsimplerthan that of TakanoandWurtele,but retainsmost of the conservation

advantagesof the Arakawa-Lambscheme.

6.5.1 The discrete equations in component form

We begin by writing the equations in a general energy-conserving form similar to that used

by Arakawa and Lamb (1981):

Oui4,t _ 1 [

Ot ( iux )i,j -°Li'jvi*+ l'j+ l -]-/_i'jvi*'j+l -t- "yi,jvi: j n u [)i,jvi*+ l,j (89a)

-}-Vi,jIt i,j_ l -- Vi,j + l Ui,j + l -- f-i + l,jUi+ l,j -t- fi'j_l'i- l'J J 1

1 r .---------_-i. . "

1 [(=_)t(u,- ut_,)+ (_)t+l(U,+,- u_)],jj

_--P lit_iTr]

1

t" J_^"z'_,j _(Ot + K_) + cpO_

i,j (_r ),,j \ _t,r ] i,jS

forj=l,JM-1,

Ovi,j,l

Ot (89b)

for j = 2, JM - 1.

Here the a, fl, 7, _, e, _o, u, and # are linear combinations of neighboring potential vorticities.

The v and # terms were added to the original Arakawa-Lamb formulation by Takano and

Wurtele to expand the stencil enough to be able to require fourth-order accuracy. The

potential vorticity at interior points is defined as

_i,j,I = ('_')-lrA_J-i, j , for j = 2, JM - 1, (90)

where the relative vorticity is given by

23
1

(i,j,l = _ [(AVy)j(6iv)i,j,t- (6j(A"xut))i,j], forj = 2,JM - 1. (91)

At the two poles, the potential vorticity is given by:

(A_)l(fl -[- _l,l)

_]l,l = (90a) 1 1 IM

2I_ E,=l ('/r Ap2),, l

(A_)jM (fJM + (JMJ) (90b)

_JM,I = 1 1 IM

where

IM

1 1 y.]_(AUxut),,1 ' (91a)

(1,l = (A_) 11M,=I

IM

1 1 _(Aux u/),,JM_I" (91b)

(JM,t = (A_)jM IM ,=1

In the code, polar values of the vorticities are replicated at all longitudes.

The discrete form of the kinetic energy, K, is presented below.

The discrete form of the Coriolis parameter is obtained by noting that f is the "planetary"

vorticity; that is, the part of the absolute vorticity due to the rotation of the reference frame

fixed to the planet. The velocity components of the rotating frame along the geographical

coordinates are:

_p = f_a cos _, _p = O,

and from (39), the components along the computational coordinates are:

uv = f2a cosq_cosx, Vp= fla cos Csinx.

On the grid, we first define Up and Vp at the p-points (and at the poles) and then interpolate

to the u- and v-points:

i

(up)i,j = [_a cos ¢ cos X, for j = 1, JM - 1, (92a)

and

J

(vp)i,j = f_acosq_sinx , for j = 2, JM - 1. (92b)

Using the same discretization as for relative vorticity, we substitute these in (91) to obtain:

_a [ A"ySi(cos_)sin x) j 6_ _b )]'",,s forj=2, JM-1, fi,j -- - _(A"xcos_cosx i

24
fta F

[AVySi(cosCNPsin(A- ANp)) (@,,:

- 5j (A"x[sinCNp cos¢--cosCNpcos(A--ANp)'

•

sine]) ]i,j'

where we have used (41b). The form used in the code is

-fla {

.fi,j- (A_)i, j sin CNP 5j(AUx cos ¢)

--cosCNp [A"y6, sin(A - ANp)+ cos(A -- ASp) / 55(A"x sine)] _..,

(93)

J

which takes into account that, in latitude-longitude coordinates, (A - )_NP) is independent

of j and ¢ is independent of i. The right-hand side of (93) is a discrete form of

_ { 0cos2 ¢ FOsin(_ - ANp) .0cosCsin¢.] }, (94) sin CNP cos CSP + cos(A -- ANp) _-¢-

cos ¢ 0¢ [ 0A

which reduces to (38). At the poles, we have:

-f_a 1 IM

f' = (A_), IM E [(A"x),,, (sinCNp COSq_l

t---=l

-- COSCNP COS()_,-- _NP) sinq_l)] , (93a)

fla 1 IM FtAux_

fJM -- (A_)jMIM _,=! [t ..J,,Jg-l(sinCsp cosCJM-I

--cOSCN P cos()_,- _NP) sint_aM-1) ] ,

(93b)

which are obtained by substituting (923) in (913,b).

6.5.2 Fourth-order vorticity advection: the choice of a, fl,...

The procedure followed by Takano and Wurtele to obtain a fourth-order version of the

Arakawa-Lamb scheme was to 1) form an equation for the second-order difference form of

the vorticity, 2) assume that the flow is non-divergent and replace the advecting velocities

by second-order differences of the streamfunction, and 3) equate terms in the resulting

non-divergent vorticity equations with corresponding terms in the fourth-order Jacobian

derived by Arakawa (1966). This procedure gives a scheme that is fourth-order only in the

advection of a second-order vorticity by the non-divergent part of the flow. We will follow

the same procedure, but without requiring that in the shallow-water equations the potential

enstrophy be conserved for a general flow. Rather, wc will require enstrophy conservation

only for the case of non-divergent flow, ,'_sin the Sadourny scheme. We will refer to the

resulting scheme as a fourth-order Sadourny.

25
As pointedout by Arakawaand Lamb (1981),the second-orderSadournyschemecorrespondsto

aS " 1

_,3 = T_(r]i,j+l -t- _i,j -1- 71i+l,j+l), for j = 1, JM - 2, (95a)

_iS; " = l-_(_i,j+l "-{-_i,j -1-?}i--l,j+l), forj = 1,JM - 2, (95b)

"TiSj = l(r_i,j+l + rli,j + rli_m,j), for j = 2, JM - 1, (95e)

6S. = l(_i,j+ 1 + _?i,j+ r/i+l,j), for j = 2, JM - 1, (95d)

where the j ranges take into account that vi_,l = vi_JM = 0, and so some of the vahles

neighboring the poles arc not used. For the second-order Sadourny scheme,

and

eS. s

,o = _i,j = O, (95e)

uS. /s. = O. (95f) _,3 = # ,J

Requiring enstrophy conservation only for non-divergent flow allows us to eliminate the

eid and g)i,j terms. Dropping these terms is very attractive because we have found them

quite troublesome. In three-dimensional calculations with the GCM, we found that they are

responsible for nmch of the problem our implementation of the Arakawa-Lamb scheme had

near the poles. These problems were eliminated when we used the second-order Sadourny,

and so we used that scheme for Version 1 of the ARIES/GEOS dynamical core. We note

that the u and # terms, required for fourth-order accuracy, are similar in that they involve u

in the u equation and v in the v equation. Indeed we found a similar, though lesser, problem

when these terms were included. This problem, however, can be partially corrected, as we

will discuss at the end of section 6.5.4.

Following the procedure outlined above, Takano and Wurtele found that the following relations nmst be satisfied to obtain fourth-order accuracy:

1

o_i,j = Ci+l,j -k -_ [16r/i,j+l + 8(_/i+l,j+l -t- ?/i,j) - 3(_]i-l,j+l -1- r}i,j+2)

-{-(r]i+l,j- 1 -[- ?']i+2,j) -- (7/i-l,j -1- 7]i,j+l -[- _']i+l,j+2 -t- fli+2,j+l)]

1

-- _(t_ nt- (P)i+I.j, for j = 1, JM - 2, (96a)

1

[_i,j = -Ci,j -{- -_ [167"]i,j+1 -Jr- 8(?'/i-ld+l -1- _l}i,j) - 3(r/i+l,j+l -t- 7}i,j+2)

d-(7]i- l,j-1 -t- r]i-2,j ) - (7"]i+l,j Jr- _i,j-I -'}- ?']i-l,j+2 "at- r]i-2,j+l)]

1

--_(e+_o)i.j, for j= 1,JM-2, (96b)

1 [16r/i,j + 8(r/i_l

")'i ,j = Ci,j q- -_ ,j q- r]i,j + 1 ) - 3 (_i j - 1 -t- r]i + l,j )

-}-(r]i-2,j+l -t- _]i-l,j+2) -- (I]i-2,j nt- _']i-l,j-1 -1- _i,j+2 -t- _/i+l,j+l )]

26
[)i ,j

1

+-_(e+(p)i.j, for j = 2, JM- 1,

1

= -Ci+l,j + -_ [16_i,j + 807i,j+ l + rli+l,j ) - 3(rli-l,j + 7h,j-1)

-I-(_i+[,j+2 "1- 7]i+2,j+1) -- (7]i+2,j "_ 7]i+l,j-[ q- ?]i,j+2 + _'_i--l,j+l)]

1

+ _(_ + _o)i+,.j, for j = 2, JM - 1,

(96c)

(9Od)

1

tli,j -_ _(_Ti+J,j -Tli-,,j), forj = 1,JM, (97a)

1

#i,j -- _-4(r/i,j-I -_li,j+J), for j = 2, JM - 1, (97b)

where the quantities C, e, and _bmay be any combination of r/s. From (97a)_ ui,1 = u,,dM = 0,

and so terms involving 7_7,1 and uT,jM in (89a) disappear.

Takano and Wurtele, following the same procedure as Arakawa and Lamb, use the choice

of C, e, and 4) to obtain enstrophy conservation for a general flow. We follow Sadourny by

setting

f-i,j : _i,j : 0 (98)

and thus abandoning exact potential enstrophy conservation for divergent flows. This leaves

only C to be determined. Note that the choice of C will only affect c_, /3, % and _5; the

quantities u and t_ are the same for all fourth-order schemes.

We will choose Ci,j fairly arbitrarily, trying only to make the scheme m_ simple as possible.

Figure 4 shows a stencil of the rlS neighboring _Li,j and the corresponding cx,/_, 7, and/i. This

stencil covers all the _s that appear in equations (96a-d). For comparison, we also show the

stencils for the second-order Sadourny, as well ,as the Arakawa-Lamb and Takano-Wurtele

schemes. The second-order Sadourny uses only the vorticities within the d,a_hed rectangle in

computing c_, fl, ?, and 6 at (i, j). The Arakawa-Lamb and Takano-Wurtele schemes use the

vortieities within the hexagonal region only two more points than the Sadourny scheme.

Equations (96a-d) involve an additional ei.qht points: the four enclosed in solid squares,

which appear underlined in (96a,d), and the four enclosed in dashed squares, which appear

underlined in (96b,e).

We. first require that C be such that the scheme uses only the Arakawa-Lamb stencil. This

is satisfied if

1

Ci,j = di,j + _ [(rl_-l,j-I + _i-2,j) -- (_;-l,j+2 + rl_-2,j+,)

"1- (_]i,j+2 "l- _i+l,j+l) -- (7]i,j-I -_ "g]i+l,j)] , (99)

where Ci,j depends only on vorticities within the desired stencil. We then choose

1

d_,j = 2-_ [(rl,-,,j + r/;,j+,)- (s_-,,j+_ + v;,j)], (100)

27
j+2

j+l

j-1

J |

Lo.J

i----i

,_° | |

t.. o.i

?_] Vi f] Vi+ 1

Ci,j (fit u_) Ci+l,j

j

rl vi rl vi+l rl

D

D

i-2 /-1 i i+1 /+2

Figure 4: Stencils of rls used to define the cx, fl, etc.

for which the scheme can be written in the particularly simple form:

3 as 1 _g (_, _, 7, _)_,_= _ (_, _, 7, )_,j- _ (_, _, 7, )_,j,

where

(101)

c_¥. 1

,,j --_ -_(t]i+l,j -1- ?]i,j+2 + ?]i-l,j+l), for j = 1, JM - 2, (102a)

&,_j_ 112(rli,j+2 + _i-l,j +_?i+l,j+l), for j = 1,JM - 2, (10=b)

1

7i N,j : --_(rli,j_ 1 +r_i+l,j--J-tli_l,j+l) , forj=2, aM-1. (102c)

_N. 1

',3 ---- -i-_(?_i+l,j+l -Jr-tli,j-I -'F _i-l,j), for j = 2, JM - 1. (102d)

6.5.3 Energy conservation and the momentum equation

To form a kinetic energy equation we multiply (89a) by (A"x u_)i,j and (89t)) by (/kvy)j v_,j,l:

28
nc ItT,j, l [OG,jVi*4_|,j+ 1 nc fli,jVi,j+

* l "Jc "_i,jVi,j . -1- _i,jVT+ l,j ]l

-- 7zi'J'l _i'j+l?ti'j+l -- lJi'jTti'j-1 l

(ZX,,)i,j 2 1 r-------_ . i 1

(6_cJ)t 2

(103)

2 jo1⁄2,,?l [A_,,,(_'_;l

-- 1_i,j, l C_i_ l ,j-17ti_ l ,j_ 1 -t'- [_i,j-- l ll, i,j_ 1 "]- "_i,j_l,i, j "1- 1

-- vi,j,l .[:+#id v I,j -- #i-!,Jvi-l,j ,] l

(_a,,)_,, 2 • 1 rv----:-vj • J ]

-

(6l_)t 2

[(_),,,,(,,,-,,,-,1+ (_),+,,),(,,,+,-,,)j,,_

-v],j d 6j(_l+Kt)+cp_l-g-- owj oj_l1 , forj=2, JM-1.

arc J i,j

(lO4)

The global integral of the kinetic energy tendency at each level can now be obtained by

adding (103) and (104) and summing over M1 i and j. Doing this, we have on the left-hand

side:

= ,_JM-_E _,,Jo5 o [_L_"2 • + A_ T--2J]' 5"_ '

• j= 1 i,j

3M-I c

E E _ )2 ou,,j,, -At'-i'J_ri'J Ot '

i j=l

(105)

since (Av),,12. ____(A2)i,JM ---- 0. Here K _ is the discrete form of the kinetic energy per unit

mass:

K_'J'I - (A2p)i,j 2 + A_ _v_ , for j = 1, JM - 1. (106)

i,j

We can perform a similar manipulation under the summation of the right-hand side, since

in the u equation all reordering is done over i, which is cyclic, and in the v equation, terms

29
at the poles that appeax when the sums are reordered as averages at the p-points contain a

factor of either v* or A 2, which are both zero at j = 1 and j = JM.

i j 3

1 1 + (Tr&)t+lA_ul(ut+l ut) i)

- E, Ej

)J Vl) j) + ( (Tr_)lA2 Vl(Vl- Vl-I + (Tr_r)l+lA2vVl(Vl+' -- ]i,j

+ Ei _,j (_ + Kt) (Siut+ _jvt) - Cp t-_ + ..,

i j 2

c [ 07r i,j 1

+ E, Ej ),j l Ot +

+ E_ Ej

- E;E_

(Kt - Kt )i,j(_fiul + tfjvt )i,j

(,5_(7r0)_ ]

(_Siu_ + _jv_)i,j + \ _la / i,j,lJ

i - (Tr_Y)lA2 _l'lTZl-I i) (_o)l 2

-j j

(108)

The second line of (108) drops out using the continuity equation, and the third line vanishes

if we take

K = K _.

Note that the a, fl, 7, 6, v, and # terms do not appear in (108). The c_, fl, 7, _i terms cancel

between the u and v equations when summed over the globe, _ssuming an appropriate

treatment at the poles. This cancellation can be clarified by considering the stencil shown

in figures 5.

In figure 5, the a, fl, 7, and 6 are shown between the us and vs they affect. Each term appears

in the equation for the neighboring u multiplied by the neighboring v and in the neighboring

v equation multiplied by the neighboring u, but with opposite sign. For example, fli,j

appears in the ui,j equation multiplied by vi*j+l, and in the vi,j+l equation (obtained by

incrementing j by 1 in (89b)) multiplied by u_*,j. (Note that the convention we use here

and in the code is that the a, fl, 7, 6 are indexed the same as the u equation in which

they appear.) The v and # terms vanish when (103) and (104) are summed over i and j,

respectively. The stencil for the v and # terms is shown in figure 6. Note that v and # are

defined at the (-points.

3O
: I :

................... Vi,j+l ................... i

fli-l,j !o_i-l,j Ifli,j i ai,j

-- ui-1,j--Ki,j -- ui,j-

_ i - 1,j i6i- 1,j I ")'i ,j i _ i ,j

, 1

..................... v i,j ......................

t

Figure 5: Stencil showing the indexing of the a, fl, % 6 factors, which appear

in the energy conserving form of the inertial terms.

Multiplying (108) by (6ta)t, summing over all layers, and using (76), we obtain the following

form of the total kinetic energy equation:

Ot i,j

l i j

_ Z(_,_),_2 { [0_r (_i(_rd)'_] cO OPt (_+v_-_j)}, (109) (A_)¢, 57+ _ _ ],J +., '-N; ,,_

t i j

where 7) is the tot_al kinetic energy dissipation:

l __ \ _ta ]l \ _ta ]l J,,j

The first term in the curly brackets in (109) involves no horizontal differencing and can be

manipulated in the same way as in the horizontally continuous case discussed in section 5.3,

equations (61) through (66). Doing this gives

-- = __,(_r _pOtA_,)_,j

Ot

l i j l i j

{ (_rO)t+l(f'_+1-Pt)+(Trd)t(Pt--151) ] (dPl_ [-' Or uTS_r +vT_ j . (110) X

Notice that the factors in curly brackets in the conversion terms in (110) and (86) are not

identical, and thus the horizontally discrete form of the equations does. not conserve energy

exactly. The difference stems from the 0 \ a_ } Vcr term in the pressure-gradient force. If

this term were discretized as

g_iiP,

rather than

0 \ O_rJ _i_r,

31
"Oi- 1,j+l Vi+I,j+I

i

Ui-l,j+l

i

(v, #)i'-l,j+l vi,j + 1

I

7Ii,j + l

i

l

i

(12, #)i,j+l

Ui- 1,j -- Ki,j -- ui,j

: I "

?'i--l,j (lJ, #)i--l,j Vi,J (V, h)i,j Vi+l,J

',

tti,_- l

,,

U i- l,j- 1

Figure 6: Stencil showing the indexing of the u and # factors, which appear

in the energy conserving form of the inertial terms.

conservation would be exact. We have not done this because the conserving form does not

allows us to maintain the irrotationality constraint on the vertically integrated pressuregradient force. In any ca.se, we will be abandoning exact energy conservation in applying the

polar filter and other modifications described in the following sections. The non-conserving

form of the equations also allows us to simplify the implementation of the economical explicit

scheme described below.

Comparing (110) with (66), we see that the horizontally discrete form of (Trwc_)i,j,t analogous

to (67) is

(111)

where (7r6, 5)i,j,l are given by (68) evaluated at each grid point. This is the form used to

evaluate the w returned by the code.

6.5.4 Instabilities of the discretized vector-invariant form

Hollingsworth-Kiillber 9 instability.

The discrete form of the equations presented above suffer from the computational instability

discussed by Hollingsworth et al. (1983). This instability is mssociated with the horizontal

differencing of the inertial terms and is common in schemes written in vector-invariant form.

It manifests itself as noise with short meridionM scales and long zonal scales in regions of

strong zonal flow.

32
Hollingsworthet al. (1983)presenta heuristicargunlentsuggestingthat aspurioussource

of momentumis producingtheinstability. Toillustrate thisargumentweconsiderthe onedimensionalshallow-waterequationswritten in vector-invariantform. Ignoringsphe.ricity

we tlave:

at (f + ()v = o,

Ov Oh 0 1 1 2

0--7,+ (f + () u = -9 Oy O*.l(-2u2 + -_v ).

Expanding terms in the v equation, we obtain the advective form:

Ov Ou Oh 01⁄2u 2 Ov

0--7.+ 'f u - u-_v = -.q &.J &.J ,, &.J"

The. argument of Hollingsworth et al. is that the underlined terms do not cancel in the

discrete form of the equations, resulting in a spurious source of momentunl. The following

linear system illustrates the effect of having such a term in the equations:

Ot f v = O,

Ov Oh Ou

+ 'f _" = -q ov - U ov'

Oh Ov

-- -- U--,

Ot Oy '

where the linearized error arising from non-cancellation is mssumed to be of the form U 0,,

The dispersion relation for this system is:

a2 = .f2 + gill 2 +, fUl, (112)

where a is the frequency, l is the wavenumber in y, and , = x/Z] -. Note that unstable

solutions occur when f and U are non-zero.

The simplest way of obtaining cancellation is to modify the form of the kinetic energy that

appears on the right-hand side of the momentum equation. Hollingsworth et el. present

such a modification for the second-order Sadourny scheme. We will follow a similar appro_h

that is currently being used in the UCLA model to remove the instability from the ArakawaLamb scheme. Rather than requiring exact cancellation, this approach requires only that

the terms cancel when the equations are linearized in cartesian coordinates. Although a

much weaker constraint than exact cancellation, this removes the instability with a smaller

modification of the kinetic energy term, and so maintains better conservation.

For one-dimensional flow,

10u 2

u(-

20y "

We wish the discrete form of the inertial terms to satisfy this relation when linearized about

a uniform zonal flow. In this case, the discrete form of u( at latitude j is:

1

U (o_j_l + _j-i +Tj + _Sj)-- U -i-_(10(j + (j+l + (j-l), (113)

33
since

cU-l = flj-I = (2Q -t- _j-1) -- _--_(_j-I 1 -t- ij -I- _j+l),

= 15_j-_+ _j-_(j+t,

1

1 5 l(j___.

---- "i-_jq_ 1 q- _-_(j --

It is worth noting that the quantity (_j-1 + flj-I + 7j + _ij) is independent of the choice of

C; thus that choice could not be uscd to optimize the form of K. Substituting

in (113), we have:

V (o_j_ l q- f_j- l -_-')'j -t- _j)

(j : --(_tj --U j-l)

1 (10(uj -?_j-1) -'1- _/,j + "/_j-1

: -U _-_ 7tj+l - - _tj-2),

= -U l((uj+l + lOuj + -1Lj_I) - (?l,j -'1-lOuj-I + uj_2)).

Since we wish to require that

v (-j-t +Zj-_ + 7_+ _j) = -(gj - gj_l),

we see that the linearized form of the kinetic energy must be

Kj = U _j + _(_,_-_ + uj+_) .

This is the linearization of

Kj=-_ 115_ u+ 62_ ] ,

where

1

_j = 2(?Zj-1 "1- 7/'j+l).

Inserting the map factors and applying this scheme in both directions, we have:

K=5_K,,+IK .',

6 6

where K _ is given by (106) and

Ki,j,l ., -- 21 A_ _ + A_ -2v , for J = 2,JM - 2,

(AI,)',J i,j

where

(Zi,j, l = _(?ti,j+l, l -t- Ui,j-l,l),

1

75i,j,l = -_ ( _Ji + l ,j,l + Vi- l ,j,l )"

(114)

(115)

(116)

(117)

34
Thesecorrectionsare applied at the interior points, j = 2, JM - 2. At the poles the

scheme for the inertial terms differs significantly and so a different modification of K is

required. However, we have not derived such a scheme and, for simplicity, have been using

the conserving form:

Ki,l,l = K_[I,t and Ki,jM-J,l = K':i,JM_IJ. (118)

It seems likely that at lemst some of the problems remaining at the poles are due to this

crude treatment.

The u term. at the poles.

As discussed above, we found that the u term given by (97a) tends to produce zonal noise

near the poles. This is similar to the problem we had encountered earlier with the e and

¢ terms in Arakawa-Lamb scheme, and seems to be related to the fact that these inertial

terms, although conservative, are not acting normal to the flow. The problem, however, is

confined to the poles, and the culprit seems to be the u factor at j = 2 and JM - 1. We

are circumventing this problem by replacing these us by interpolations between interior and

pole values (note that u vanishes at the pole):

(Apy) 

Pi,2 = _'i,3 (APy)l -k- (APy)2' (119a)

(APy)JM-I

(l19b)

---- b'i,JM_ 2 (Apy)jM_ 1 -[- (APy)jM_ 2 '

Vi,JM- 1

since t)i,l = t*i,JM :=a O.

7 Polar Filters

In the dynamical core, we apply a polar Fourier filter to the tendencies of all prognostic

variables. The purpose of the polar filter is to avoid linear computational instability due to

the convergence of the meridians near the poles. The filter acts poleward of 45 ° latitude, and

its strength is gradually increased towards the pole by increasing the number of affected

zonal wavenumbers and the amount by which they are damped. The polar filter is also

applied to the diagnostic w.

Let _bi, j denote a single level of any of the tendencies to be filtered. Its zonal Fourier

expansion is:

IM/2

¢i,J = E _m,jexp(--zrnAi),

ra=O

where Cm,j is the complex amplitude of the ruth zonal wavenumber and z = v/AT The

filtered tendency, _bi,j, is

IM/2

¢i,j = E era,j@.,j exp(-*mAi), (120)

rn=0

35
wherethe filter coefficientsaregivenby:

n

\cos(¢450)

(121)

The rationale for the form of the filtering coefficients is the same as in Arakawa and Lamb

(1977). Crudely speaking, for n = 1 this filtering "slows-down" the zonal propagation of

each zonal component enough to satisfy the same stability criterion ,as is required by the

shortest wave at a particular latitude in this case 45 °. Arakawa and Lamb (1977) use the

meridional grid size, rather than the zonal grid size at 45 ° latitude, and they use the filter

more selectively, treating only the gravity waves by smoothing only the mass divergence and

the pressure-gradient term, rather than the entire tendencies. We found that applying it to

the entire tendencies made the model more stable without adversely affecting the results.

Also, if the filter is applied more selectively, one must be careful of how it interacts with the

inertial instabilities discussed in section 6.5.4. We. have also found it useflfl to take n = 2,

which provides additional filtering near the pole and seems to enhance the stability.

Note that in this form the polar filter is only appropriate for grids with equal spacing

in longitude. Modifying the filter for more general grids should be straightforward; for

example, one can first interpolate in the zonal direction to a uniform grid as fine as the

finest spacing on the non-uniform grid and then apply the same scheme. Another possibility

is to design a local filter that approximates the Fourier filter for a uniform grid. We are

currently developing such schemes.

8 Time Differencing

As discussed in the introduction, the task of the dynamical core is limited to computing

the dynamical contributions to the tendencies of all state variables. With time continuous,

these depend only on the instantaneous state. For some explicit time differencing schemes,

such a.s the leap-frog t)r the Euler-backward, these tendencies depend on the state at a

single time level and are independent of the particular time-differencing chosen. But for

multi-level explicit schemes and for all implicit schemes, the tendencies depend on the time

differencing and on the size of the time step.

In the ARIES-GEOS core, the argument list includes two time levels of the state variables

and the time step. This argument list can accommodate the most commonly used, explicit

or implicit, timc-diffcrcncing scheme. In Version 2 we support only explicit schemes. When

the specified time step is less than or equal to zero, the core returns the explicit, single-level

tendencies. When a positive time step is specified, the core returns "economical explicit"

tendencies, using the scheme of Brown and Campana (1978), ms described below. We plan

to inchlde semi-implicit differencing in future versions of the core.

36
8.1 The Brown-Campana Scheme

Brown and Campana (1978) proposed an explicit sctmme to use in conjunction with leapfrog differencing which relaxes somewhat the instability condition fi)r gravity waves. The

scheme does not affect the advection or inertial terms. The idea is to average the pressure

gradient fi_rce over three time levels. This can be accomplished very simply when using a

standard leap-frog scheme by updating the mass field (potential temperatures and surface

pressure in our case) first., then computing the pressure gradient at the three levels and

taking an approt)riate average. An even simpler strate.gy, which is equivalent when the

pressure gradient force is linearized, is to average the three time lcve.ls of the mass field and

only compute the pressure gradient once, using the averaged values. We take the latter

approach, simplifying it still fllrther by using the averaged mass field only for selected parts

of the pressure gradient calculat.ion.

We finite difference the moment, urn equation in time as follows:

v "+'-v .... l [ 0_r'_"

At =- -vq_" - _c,,O--_ ) Vrc* +( all other terms )', (122)

where the '-)* are obtained fl'om the hydrostatic equation:

q'[x, = _s + q,0,.x,t(/;,.M+l -- PLM) '_, (123)

^ l

_" = _I)7+1 + %[(P_+1 - Pl+l)"Ol+l + (/5l+1 - Pl)'"_']_ for l = 1, LM - 1, (124)

and

_* = _'. (125)

Here t,he ( yt denotes the three time level average

W' = _E(vr '_+l + vr"-1) + (1 - 2cz_7)rr '_. (126)

As pointed out by Brown and Campana (1978), the averaging parameter, aE, must be

chosen careflflly to maximize the time step. The best, value depends on the strength of the

Robert (1966) time filter used to control the splitting of solutions with the leap-frog. We use

c_E = .2475, which is near optimal for a time-filtering parameter of .05. We emphasize that

the Robert filter is not done by the dynamical core, only the Brown-Campana smoothing

of the pressure gradient that enters the momentum tendencies returned by the core.

Another point to note is that the dynamical core will use a completely updated mass field

at time level n + 1 only if the tendencies due to all other (non-dynamical) processes are

input. However, it is not essential to do this to obtain the time step improvement, since

what matters is the treatment of the gravity waves.

37
9 Using Subroutine DYCORE

The dynamical core is accessed by a call to FORTRAN subroutine DYCORE. Inputs to

the routine include the state v_iables, various parameters that specify the horizontal and

verticM grids, and some physical constants. The outputs of the routine are the updated

time tendencies of all state variables. The state variables themselves are not modified.

A typical GCM application would use the following procedure at each time step:

• Compute the "physics" tendencies, storing them in UOI, VOI, POI, QOI (PII would

normally be zero, since sources or sinks of mass are usually ignored).

• Call DYCORE, which would add the dynamical contributions to the tendencies.

• Add the contribution of any horizontal diffusion.

* Perform the time stepping, including any time filtering that may be needed, to obtain

updated values of the state variables.

9.1 The Argument List

The following is a listing of the SUBROUTINE declaration and the "banner" that is inchlded

at the top of the DYCORE source code. The banner describes the argument list and the

toutine s storage requirements.

SUBROUTINE DYCDRE ( IM,JM,JDIM,LM,SIG,PTDP,KM,DT,

OMEGA, CP, RGAS, AE,

PHS ,PKH,

PIB ,UOB ,VOB ,POB ,QOB,

PIM ,UOM ,VOM ,POM ,QOM,

PII, UOI, V01, P01, QOI,

0MG, VOR, DIAG )

C INPUT ARGUMENT DESCRIPTION

C IM ....

C JM ....

C JDIM ..

C LM ....

C SIG ...

C PTOP ..

C KM ....

C DT ....

Number of Grid Intervals in Zonal Direction

Number of Grid Intervals in Meridional Direction

Meridional (Second) Dimension of Input Fields

Number of Vertical Levels

(LM+I): Sigma at Interfaces. SIG(1)=O; SIG(LM+I)=I

Model Top Pressure at Sigma = 0.0

Number of Scalars, Including H2D, but not There.

Time-Step fron n-I to n+l (Seconds)

OMEGA..

CP .....

RGAS...

AE .....

Rotation rate (tad/see)

Specific heat at constant pressure (J/(kg K))

Gas content (J/(kg K))

'Earth' radius (meters)

PHS ... (IM,JDIM): Surface Geopotential (m * m/see**2)

38
PEH .. (IM,JDIM,LM+I): (P/POO)**KAPPA (Not used if PTOP:O)

PIB

UOB

VOB

POB

QOB

C PIM

C UOM

C VOM

C POM

c QOM

.. (IM,JDIM): Mass (Psurf-Ptop) mb at Current Time-Level

.. (IM,JDIM,LM): Zonal Wind m/s at Current Time-Level

.. (IM,JDIM,LM): Meridional Wind m/s at Current Time-Level

... (IM,JDIM,LM): Potential Temp. K at Current Time-Level

.. (IM,JDIM,LM,KM): Scalar Fields at Current Time-Level

.. (IM,JDIM): Mass (Psurf-Ptop) at Previous Time-Level

.. (IM,JDIM,LM): Zonal Wind at Previous Time-Level

.. (IM,JDIM,LM): Meridional Wind at Previous Time-Level

.. (IM,JDIM,LM): Potential Temperature at Previous Time-Level

.. (IM,JDIM,LM,KM): Scaler Fields at Previous Time-Level

C OUTPUT ARGUMENT DESCRIPTION-- Tendencies are in per second.

PII ... (IM,JDIM): Updated Surface Pressure Time-Tendency

UOI ... (IM,JDIM,LM): Updated Zonal Wind Time-Tendency

VOI ... (IM,JDIM,LM): Updated Meridional Wind Time-Tendency

POI ... (IM,JDIM,LM): Updated PI-Weighted Theta Time-Tendency

QOI ... (IM,JDIM,LM,KM): Updated PI-Weighted Scalar Time-Tendency

OMG ... (IM,JDIM,LM): Omega Diagnostic (mb/sec)

VOR ... (IM,JDIM,LM): Vorticity Diagnostic (I/sec)

DIAG .. Logical (On/Off) Flag for Diagnostics

C NOTES:

C (i) JDIM is to be used for DIMENSION Purposes ONLY

C Vectorization will be performed over IM*JM.

C (2) The Vertical Layers are numbered from TOP(1) to BOTTOM(LM).

C (3) All Time-Tendencies are INCREMENTED (bumped).

C The Momentum Time-Tendencies ARE NOT mass-weighted.

C ]'he Potential Temperature and Scalar Time-Tendencies ARE

C mass-weighted (by PI).

C (4) JM is 180 degrees divided by the meridional grid size.

C (5) UXX(I,J) are located half a grid interval EAST of PXX(I,J).

C VXX(I,J) are located half a grid interval SOUTH of PXX(I,J).

C (6) If PTOP>O, the PKH MUST be defined.

C (7) The previous time level fields (PIM,UOM,etc) are used for the

C economical explicit calculation done in conjunction with

C leap-frog steps. If you are not doing leap-frog or do not

C wish to have economical explicit tendencies, pass DT<=O.

C SPACE REQUIREMENTS:

C (I) Takes IM*JM*I9+4*JM+2*IM words from the heap for STATIC storage;

C these are kept throughout the run.

C (2) Takes IM*JM*(LM+25) ÷ 3*LM + I words from the heap for DYNAMIC

C storage when PTOP=O; for PTOP!=O, add IM*JM*LM words.

C All of this storage is freed before returning.

9.2 Comments

• The argument JDIM is provided in case the calling program ll,_ the state variables

organized by level rather than variable, or has extra latitudes.

39
With PTOP set to zero, the code economizes by raising only the surface pressure to

the _:. If PTOP is non-zero the code expects the (P)_ at the interfaces to be passed

-poin the array PKH. This is done to avoid costly recomputations of this factor, which is

usually needed outside the dynamical core. If PTOP is zero, the input variable PKH

is completely ignored; it need not, be defined, or even allocated as an array. In this

case we are not concerned about duplicating the exponentiation.

Tile code sets the grid according to the values of IM and JM. This is done on the first

call or whenever IM or JM change from the previous call, and the grid is saved in a

dynamically allocated static space. The vertical grid is inexpensive to compute and

so it is reset oll each call.

In the current version the time step is used only for computing Brown-Campana

tendencies. It is Mso used as a flag for these tendencies, which are only computed

when DT is greater than zero. If DT is less than or equal to zero, the code returns an

explicit tendency computed with the PIB, UOB, etc. state variables, and the PIM,

U()M, etc. are ignored.

Access to the grid rotation or irregular grid features is not provided through the

argument list in this version. The location of the pole and the zonal symmetry can be

oh;raged by altering the appropriate parameters in subroutine GRIDH. An irregular

grid in 1;ttitude can be specified by modifying the definition of the ¢I>j, also in routine

GRIDH. Irregular grids in longitude are not supported in this version.

A message passing version is also available.

10 Final Remarks

We have presented a detailed description of the Enlerian grid-point dynanfics currently being

used as the "dynamical core" of global atmospheric modeling at Goddard. This dynamical

core is written in a nlodular fornl that allows easy ilnt)lementatioIl , assuming the GCM's

variables are already stored on an Arakawa C-grid. 111 addition to a detailed discussion of

the numerics, inchlding the introduction of a new fourth-order sdlenle, we have described

the usage _f the routine..

The next w_'rsion of the core will inchlde:

• A flmrt, h-order scheme for the horizontal advection of potential temperature and tracers

• _. flfll telescoped grid in latitude and longitude

• An improved trcatnlent of the poles in the monlentum equation

4O
All of these changes are in fairly advanced stages of development.

The FORTRAN source code is available over the network. To obtain access to it, or for

filrther information, send e-mail to:

Max.J.Suarez(@gsfc.n,xsa.gov or Lawrence.L.Takacs(_gsfc.n,'_sa.gov

41

References

Arakawa, A., 1966: Computational design for long-term nmnerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I. J. Cornput.

Phys., 1, 119-143.

Arakawa, A., 1972: Design of the UCLA generM circulation model. Tech. Report No.Z

Dept. Meteor., UCLA, Los Angeles CA 90024. Mon. Wea. Rev., 109, 18-36.

Arakawa, A. and V.R. Lamb, 1977: Computational design of the basic dynamical processes

of the UCLA general circulation model, in Methods m Cornputation.al Physics, J

Chang ed., Academic Press, 17, 173-265.

Arakawa, A. and V.R. Lamb, 1981: A potential enstrophy and energy conserving scheme

for the shallow water equations. Mon. Wea. Rev., 109, 18-36.

Arakawa, A. and M.J. Suarez, 1983: Vertical differencing of the primitive equations in sigma

coordinates. Mon. Wea. Rev., 111, 34-45.

Brown, J.A. and K. Campana, 1978: An economical time-differencing system for numerical

weather prediction., Mon. Wea. Rev, 106, 1125-1136.

Burridge, D.M. and J. Haseler, 1977: A model for medium range weather forec`a_ting

adiabatic formulation, Tech. Report. No. 4, European. Cen.ter for Medium Ran.qe

Weather Forecasts, Brachnell, Berkshire, UK.

Hollingsworth, A., P. K£11berg, V. Renner, and D.M. Burridge, 1983: An internal symmetric

computational instability, Quart. J. R. Met Soc. , 109, 417-428.

Kalnay, E., M. Kanamitsu, J. Pfaendtner, J. Sela, M. Suarez, J. Staekpole, J. rlSmcillo, L.

Umscheid, and D. Williamson, 1989: Rules for the interchange of physical parameterizations, Bull. Am Met. Soc., 70, 620-622.

Phillips, N.A., 1957: A coordinate system having some special advantages for numerical

foreca.sting. J. Meteor., 14, 184-185.

Phillips, N.A., 1974: Application of Arakawa's energy conserving layer model to operational numerical weather prediction. Office Note 104, National meteorological Center, NWS/NOAA 40pp.

Robert , N.A., 1966: The integration of a low-order spectral form of the primitive meteorological; equations. J. Meteor. Soe. Japan, 44,237-245.

Sadourny, R., 1975: The dynamics of finite difference models of the shallow water equations,

J. Atmos. Sci., 32, 680-689.

Schubert, S. D., J. Pfaendtner and R. Rood, 1993: An ,assimilated data set for Earth Science

applications, Bull. Am Met. Soc., 74, 2331-2342

Takacs, L. L. , A. Molod, and T. Wang, 1994: Goddard Earth Observing System (GEOS)

43

•i .'A,):;
General Circulation Model (GCM) Version 1. NASA Technical Memorandum 104606

Volume 1, Goddard Space Flight Center, Greenbelt, MD 20771, 97 pp.

Takano, K. and M.G. Wurtele, 1982: A fourth-order energy and potential enstrophy conserving difference scheme, Air Force Geophysics Laboratory Report AFGL-TR-82-0205,

(NTIS AD-A126626), AFGL Hanscom AFB, Mmssachusetts, 01731, 85pp.

44
Previous Volumes in This Series

Volume 1

September 1994

Documentation of the Goddard Earth Observing System

(GEOS) general circulation model - Version 1

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

Volume 2

October 1994

Direct, solution of the implicit fornmlation of fourth order

horizontal diffusion for gridpoint models on the sphere

Y. Li, S. Moorthi, and J.R. Bates

Volume 3

December 1994

An efficient, thermal inDared radiation parameterization for

use in general circulation models

M.-D. Chou and M.J. Suarez

Volume 4

January 1995

Documentation of the Goddard Earth Observing System

(GEOS) Data Assimilation System - Version 1

James Pfaendtner, Stephen Bloom, David Lamieh,

Michael Seablom, Meta Sienkiewiez, James Stobie,

and Arlindo da Silva

45


REPORT DOCUMENTATION PAGE Form Approve_

OMB No. 0704-0188

Public reporting burden for this collectionof informalion is estimated to average 1 hour per response, includingIhe timefor reviewing instructions, searching existing data sources,

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

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

April 1995 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

Technical Report Series on Global Modelingand Data Assimilation

Volume 5 - Documentation of the ARIES/GEOS Dynamical Core:

Version 2 C - NAS5-32332

6. AUTHOR(S) RTOP 578-41-45-20

Max J. Suarez and Lawrence L. Takacs

7. PERFORMINGORGANIZATIONNAME(S)ANDADDRESS(ES)

Laboratory for Atmospheres

Climate and Radiation Branch

Data Assimilation Office

Goddard Space Flight Center

Greenbelt, Maryland 2077 l

9.SPONSORING/ MONITORINGADGENCYNAME(S)AN'DADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

Code 913

8. PEFORMING ORGANIZATION

REPORT NUMBER

95B00069

10, SPONSORING/MONITORING

ADGENCY REPORT NUMBER

TM- 104606, Vol. 5

11. SUPPLEMENTARYNOTES

Suarez: Goddard Space Flight Center, Greenbelt, Maryland

Takacs: General Sciences Corporation, Laurel, Maryland

12a. DISTRIBUTION I AVAILABILITY STATMENT

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)

12b. DISTRIBUTION CODE

A detailed description of the numerical formulation of Version 2 of the ARIES/GEOS "dynamical

core" is presented. This code is a nearly "plug-compatible" dynamics for use in atmospheric general

circulation models (GCMs). It is a finite difference model on a staggered latitude-longitude C-grid.

It uses second-order differences for all terms except the advection of vorticity by the rotation part of

the flow, which is done at fourth-order accuracy. This dynamical core is currently being used in the

climate (ARIES) and data assimilation (GEOS) GCMs at Goddard.

14. SUBJECTTERMS

Atmospheric dynamics, general circulation models, GCM, finite differences,

numerical methods, primitive equations.

17. SECURITY CLASSIFICATION

OF REPORT

Unclassified

18. SECURITY CLASSIFICATION

OF THIS PAGE

Unclassified

NSN 7540-01-280-5500

19. SECURITY CLASSIFICATION

OF ABSTRACT

15. NUMBER OF PAGES

58

16. PRICE CODE

20. LIMITATION OF ABSTRACT

UL Unclassified

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39.18

298-102