5

,?

• i  ̧

NASA Technical Memorandum 104606, Vol. 1

Technical Report Series on

Global Modeling and Data Assimilation

Max J. Suarez, Editor

Goddard Space Flight Center

Greenbelt, Maryland

Volume 1

Documentation of the

Goddard Earth Observing System (GEOS)

General Circulation Model - Version 1

Lawrence L. Takacs

Andrea Molod

Tina Wang

General Sciences Corporation

Laurel, Maryland

•i:/

National Aeronautics and

Space Administration

Goddard Space Flight Center

Greenbelt, Maryland

1994
"06E0- I _9 (I 0_) '17_6_-0601 g GIN 'slqg!oH tuna!qlu!q 'pgo_I gut.pug] ogP!J:'II_[ 008

'uoDmu.iojui oagdgoJo v Joj JOlUOD VSVN oql tuoJj Olqel!gae S! uoggo!Iqnd s!q_L

, [>' ,

<'.<(, '!

i/i <i_

/i %11

k

i

<<
Abstract

This report is a documentation of Version 1 of the Goddard Earth Observing System

(GEOS) General Circulation Model (GCM) developed by NASA's Data Assimilation

Office and the Climate and Radiation Branch. The report is separated into three parts.

Part I describes the overall features of the GCM, including numerical schemes for the

hydrodynamics, sub-grid scale physical parameterizations, and boundary conditions.

Part II provides a comprehensive description of the diagnostics available within the

model, including the methodology for their selection and retrieval. Part III concludes the

report with a User's Guide, which includes the location of frozen versions of the GEOS

system, instructions on how to create user-defined applications for running the GEOS

system, and an overview of established production applications and output utilities

available to the communitiy.

:•,! _.i ̧

(,

°°.

111

PAGE CLANK NOT FN.MED

i,

!i i__

_'i _.

Contents

List of Figures

I GEOS-1 GCM: A Descriptive Documention

1 Introduction and Model Lineage

1.1 Relevant Model Documentation .........................

2 Atmospheric Dynamics

2.1 Horizontal and Vertical Discretization .....................

2.2 Time Integration Scheme ............................

2.3 Smoothing / Filling ...............................

Vlll

1

1

1

2

4

5

7

3 Atmospheric Physics 9

3.1 Moist Convective Processes ........................... 9

Sub-grid and Large-scale Convection .................. 9

Cloud Formation ............................. 11

3.2 Radiation ..................................... 12

3.2.1 Shortwave Radiation ............................ 13

3.2.2 Longwave Radiation ........................... 13

3.2.3 Cloud-Radiation Interaction ....................... 14

3.3 Turbulence .................................... 15

z

• :::i ̧

i ,_: i_

3.3.1 Planetary Boundary Layer ........................ 16

4 Boundary Conditions 16
i!il:_i

4.1 Orography ..................................... 16

4.2 Land-Sea-Ice Surface Type ........................... 17

4.3 Sea Surface Temperature ............................ 17

4.4 Sea Surface Roughness .............................. 21

4.5 Albedo ....................................... 21

4.6 Sea Ice ....................................... 21

4.7 Snow Cover .................................... 21

4.8 Land Surface ................................... 21

II GEOS-1 GCM: Diagnostics 23

5 Introduction 23

6 Diagnostic Overview 23

7 GEOS GCM Diagnostic Menu 26

8 Diagnostic Description 30

III GEOS-1 GCM: User's Guide 65

9 Introduction 65

10 GEOS Superstructure 65

10.1 GEOS Applicutions ................................ 65

10.2 Subroutines GCMDRV & RESTART ...................... 67

vi
11 GEOS Production 71

11.1ProductionNamelists .............................. 73

11.2PhoenixFormat.................................. 77

11.3GeopotentialHeightsonSigmaSurfaces.................... 81

11.4Interpolationfrom Sigmato PressureSurfaces................. 83

11.5Sea-LevelPressure................................ 88

12 GEOS Unix Script 90

12.1Creatingthe GEOSGCM Namelist....................... 90

12.2Creatingthe DIAGSIZESubroutine ...................... 91

12.3RemotecopyfrozenGEOSLibrary ....................... 92

12.4CompileUserApplicationand Output Routines................ 93

12.5LoadandRun GEOSSimulation ........................ 94

References 97

/!i_(i

vii
. /

List of Figures

1

2

3

Sigma-level pressures used in the 20-level GEOS-1 GCM ........... 6

Shapiro filter response function used in the 2° x 2.5 ° GEOS-1 GCM ..... 8

Comparison between the Lanczos and mth-order Shapiro filter response functions for m = 2, 4, and 8 ............................. 18

GEOS-1 GCM Land/Water Mask and Topography used at 4 ° x 5° resolution. 19

GEOS-1 GCM Land/Water Mask and Topography used at 2° x 2.5 ° resolution. 20

Vertical placement and index notation for sigma levels in the GEOS GCM . 83

Index notation for vertical interpolation of Geopotential Heights ...... 84

viii
i ̧ i:i

.)

_i  ̧ :i/

:i:_  ̧/•/ ̧!

:r

)(i:::)

Part I

GEOS-1 GCM: A Descriptive

Documention

1 Introduction and Model Lineage

The Goddard Earth Observing System-1 (GEOS-1) General Circulation Model (GCM)

was developed by the Data Assimilation Office (DAO) at the Goddard Laboratory for

Atmospheres (GLA), in collaboration with the Climate and Radiation Branch, for use in

the system being developed to analyze EOS data. The GEOS Data Assimilation System

(DAS), currently using the GEOS-1 GCM in conjuction with an Optimal Interpolation

(OI) analysis scheme, has been used to produce a multi-year global atmospheric data set

for climate research (Schubert et al., 1993). The GEOS-1 GCM has also been used to

produce multiple 10-year climate simulations as part of the DAO's participation in the

Atmospheric Model Intercomparison Project (AMIP)sponsored by the Program for Climate

Model Diagnostics and Intercomparison (PCMDI) (see Gates, 1992). The GEOS-1 DAS has

been used operationally to provide scientific flight guidance during NASA's participation in

the Measurements for Assessing the Effects of Stratospheric Aircraft (MAESA) experiment.

The earliest predecessor of the GEOS-1 GCM was developed in 1989 based on the "plugcompatible" concepts outlined in Kalnay et al. (1989), and subsequently improved in 1991

(Fox- Rabinovitz, et al., 1991, Helfand et al., 1991). The plug-compatibility of the physical parameterizations together with the plug-compatible "Dynamical Core" introduced by

Suarez and Takacs (1994) facilitate the development and testing of new algorithms. Together the DAO and the Climate and Radiation Branch at GLA have produced a library

of physical parameterizations and dynamical algorithms which can be utilized for various

GCM applications.

1.1 Relevant Model Documentation

A description of the overall features of the GEOS Data Assimilation System (DAS) used by

the DAO may be found in Schubert et al. (1993). A comprehensive documentation of the

Aries/GEOS Dynamical Core incorporating the horizontal and vertical discretization and

finite-difference schemes used is given in Suarez and Takacs (1994). The Relaxed ArakawaSchubert cumulus convectNe parameterization and the re-evaporation of falling rain are

based upon the works of Moorthi and Suarez (1992) and Sud and Molod (1988), while the

radiative processes are described by Harshvardhan et al. (1987). The planetary boundary
layer(PBL) andthe upperlevelturbulenceparameterizationarebasedonthe level2.5closuremodelof HelfandandLabraga(1988)andHelfandet al. (1991).Additionalpostscript

copiesof this document may be obtained through anonymous ftp from hera.gsfc.nasa.gov.

The path and filename of the document is:

pub/gcm/geos 1.0_gcm.doc.ps

2 Atmospheric Dynamics

The momentum equations used in the GEOS-1 GCM are written in the "vector invariant"

form, as in Sadourney (1975) and Arakawa and Lamb (1981), to facilitate the derivation

of the energy and potentiM enstrophy conserving differencing scheme. The thermodynamic

(potential temperature) and moisture (specific humidity) equations are written in flux form

to facilitate potential temperature and moisture conservation. A complete description of

the finite-difference scheme used can be found in Suarez and Takacs (1994).

The GEOS GCM uses a a coordinate defined by

P - PT

a -- , (1)

where p is the pressure, _ = p_ - PT, P_ 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 a coordinate proposed by Phillips (1957).

With this vertical coordinate, the continuity equation becomes

0r

--=ot -v_.(_v)- 0(_) oe ' (2)

where v is the horizontal velocity vector. Integrating (2) and assuming & = 0 at p = Ps and

P -_ PT, we obtain the forms used in the model:

07r 0t - f01 v_. (_v) de (3)

and

(_r&) = -e-070r- ]0°V_-(_rv) de. (4)

The equation of state for an ideal gas is _ = RT/p, where _ is the specific density, T is the

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

2
 ̧/' ̧5!

4

where 0 -= T/P is the potential temperature, P - (p/po) _, _ = R/cp, Cp is the specific heat

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

used OF = _E and the relation

p

_r O" G

For the time being virtual effects are neglected.

The hydrostatic equation is

0¢

Op

where _ is the geopotential. Using (5) and (6), this can be written:

0a- lra=-cp0_ _ _=-c'0 -JJa _"

From (7) we obtain

0¢

_-fi = -cp0,

which, following Arakawa and Suarez (1983), is the form used in the model.

(7)

(8)

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

conserving differencing scheme:

0(_0) 0(_0) Lq

Ot - -V_. (_rv0) Oa + cpP' (9)

where Q is the diabatic heating per unit mass.

In addition to the equations of motion, the Aries/GEOS Dynamical Core computes tendencies for an arbitrary number of atmospheric contituents, such as water _rapor and ozone.

These are also written in flux form:

°(_q(k)) -v_.(_vq (k)) °(_q(k))

Ot - Oa + 7r8 (k), (10)

where q(k) is the specific mass of the kth constituent, and 8 (k) is its source per unit mass

of air.

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

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

0v 0v

Ot - -(I+()kxv-&oGG-Vo(@+K)-%0V_,P----

= -r]kx(Trv)-.0v --(dP) _-V_(¢+K)-c_0_/_

g OT

0_,' (11)

g OT

VTr (12)

7r OG'
)i _ _i,,) _

ii

where

(/+ 0

_r

is an "externM" potential vorticity, f is the Coriolis parameter, k is the unit vector in the

vertical,

(--V_xv

is the vertical component of the vorticity along a surfaces,

is the kinetic energy per unit mass, g is the acceleration of gravity, and 7" is the horizontal

frictional stress.

2.1 Horizontal and Vertical Discretization

L

/

The GEOS-1 GCM is constructed in the horizontal using the staggered Arakawa C-grid

and employs Version 1 of the Aries/GEOS Dynamical Core for the finite-differencing algorithm. The Aries/GEOS Dynamical Core is a plug-compatible dynamics module which

is used at both the DAO in the GEOS GCM as well as at the Climate and Radiation

Branch in the Aries GCM. The Aries GCM has been extensively used in many climate and

coupled ocean/atmosphere simulations (eg. Schubert et al. 1993, Higgins and Schubert,

1993). Version 1 of the Aries/GEOS Dynamical Core is the second-order energy and potential enstrophy conserving scheme of Sadourny described by Burridge and Haseler (1977),

while Version 2, to be used in subsequent versions of the GEOS GCM, is the fourth-order

formulation described in Suarez and Takacs, 1994. This scheme conserves total energy

and potential enstrophy for the non-divergent component of the flow in the shallow water

equations.

The Aries/GEOS Dynamical Core uses a Lorenz or unstaggered vertical grid in generalized

sigma coordinates. The vertical differencing scheme is that of Arakawa and Suarez (1983)

which ensures that:

The pressure gradient force generates no circulation of vertically integrated momentum along a contour of surface topography

The finite-difference analogues of the energy-conversion term have

the same form in the kinetic energy and thermodynamic equations

The global mass integral of the potential temperature is conserved

under adiabatic processes

The hydrostatic equation for the lowest thickness has a local form

The hydrostatic equation is exact for vertically isentropic atmospheres

:i!i5 /i

,•_ • ..... iL
iii ,!:!

,i  ̧ _.,:

:i:  ̧- • .

_, :_i' I

i _

' i

"!]i: _ ' , :

_ i• /

• i _

/

The pressure-gradient force is exact for three-dimensionally isentropic atmospheres

The GEOS-1 AMIP simulation was run using a global 4° x 5°latitude/longitude grid in

the horizontal, together with 20 sigma levels in the vertical. The GEOS-1 DAS 5-year

re-analysis was run using 2° x 2.5 ° horizontal resolution. The vertical distribution of the

sigma levels is such as to provide enhanced resolution in the planetary boundary layer and

at the jet level. The model top pressure was prescribed to be 10 mb with a rigid lid. For

a surface pressure of 1000 mb, the lowest sigma level is at 994 mb while the highest sigma

level is at 19 mb. There are 5 sigma levels below 800 mb, and 7 sigma levels above 200 rob.

The sigma-level pressure used in the 20-level GEOS-1 GCM are shown in Figure 1.

2.2 Time Integration Scheme

The GEOS GCM has the ability to use the Matsuno time integration scheme or the Leapfrog

time integration scheme together with an Asselin (1972) time filter. The GEOS GCM multiyear simulations have all been run using the Leapfrog time scheme together with the Asselin

averaging parameter equal to 0.05. The GEOS-DAS 5-year Re-analysis used the Matsuno

time scheme. The GEOS GCM employs a somewhat unique method for incorporating

adjustments due to diabatic processes (ie, moist convection, radiation and turbulence) and

the analysis increments during an assimilation. At every model time step, all prognostic

fields are updated due to both dynamical and sub-grid scale diabatic processes, shown

schematically using the Leapfrog time scheme for an arbitrary prognostic field q as:

where

= qn-1+ 2Z t(0qn 

\_ )tota_ (13)

_ (Oqn-l_

-_ J total -- \ O--t ] Dynamics _ -_ ] Shapiro Filter

-_ (°-_ 0Moist Processes +(o0Turbulence

+ (o_ 0Long_aveRadiation+_,o•cos¢z (o_ 0S_ort_aveRadiation

[_Oq

+

_-_ ) Analysis Increment (14)

with So defined as the solar contant, Ra as the earth-sun distance in Astronomical Units,

and cos¢_ as the cosine of the zenith angle. The Dynamics and Shapiro Filter time tendencies are updated every model time step. The diabatic time tendencies are computed

5
(qtu) e Jnsse_d

§

§

!:

i , i , , i , , , _

(_qtu) gJnssgJd

_-----7_ : _ _-

(_quJ) e Jnssg Jd

ob

©

r..q

_>

6

gr_

i>

¢..o
• /

/

i:/ -

/

(updated) at a time step appropriate to the physical parameterizations, and are held constant between Physics calls. The time tendency for Moist Convection is updated every 10

minutes, and for Turbulence every 30 minutes. The time tendency for Longwave Radiation

is updated every 3 hours. Shortwave Radiation is updated once every three hours assuming

a normalized incident solar radiation, and adjusted at every model time step by the true

incident radiation. During GEOS-1 DAS assimilations, the Analysis Increment is updated

every synoptic time period, or 6 hours. By gradually incorporating the diabatic adjustments during the model integration, shocks and dynamical imbalances are greatly reduced

(cf. Bloom, et al. 1991).

2.3 Smoothing / Filling

An eighth- (sixteenth-) order Shapiro (1970) filter adjustment is computed for the winds,

potential temperature, and specific humidity in conjunction with the 2° x 2.5 ° (4 ° x 5° )

GEOS GCM in order to globally damp small-scale dispersive waves. In order to reduce

the dynamical imbalances generated by intermittent use of the full filter, a Shapiro filter

tendency for any quantity q is defined as:

Oq) qF_ q (15)

-_ Shapiro Filter -- T

where qF is the quantity after application of the full Shapiro filter, q is the unfiltered

quantity, and _- is an adjustable time scale. Thus, a fl'action /,t of the full Shapiro filter is

T

incorporated at each model timestep. In the GEOS-1 GCM, the time scale _-is equal to 1.5

hours, and has been chosen so as to remove the two-grid interval wave in about six hours.

The Shapiro filter response function used in the 2 ° x 2.5 ° GEOS-1 GCM is shown in Figure

2. The GEOS-1 GCM also uses a high-latitude Fourier filter to avoid linear instability due

to violation of the CFL condition for the Lamb wave and internal gravity waves near the

poles. This polar filter is only applied to the time-tendencies of the prognostic fields, ie.

winds, potential temperature, specific humidity, and surface pressure.

For the moisture equation, negative values computationally generated through advection

are filled by "borrowing" fl'om below while conserving the vertically integrated moisture,

ie.,

I pqdz] = pqdz

Jle+l ] final \Jle+l initial

Using the hydrostatic relation pSz = -_ = -_-Sa, we may write

g g

(16)
8th-Order Shapiro Filter Response Function

0°9"

0.8 .......................................

\

0.7 ........................................

0.6 .........................................

\

0.5

0.4...............,..........:-...................

...... ,...... ..... . . . . .

i fullShopiroGtter l&3Applicotions

fred[onQI ShopiroQllel 2&61_urs

FfQc[iona[ ShopiroQ[terl[imestep

\i ............

5 10 15 20 25 30 ,:35 40 45 50 55 60 65 70

Wovenumber

Figure 2: Shapiro filter response function used in the 2 ° x 2.5 ° GEOS-1 GCM.

8

• (' •
Iv)_rqda = Ir _rqda ) \dle+l final \dle+l initial

Approximating equaton (17) by

(TrqlAal + 7rql-lAffl-1)final = (TrqlAal + rql-lACrl-1)initial

(17)

(18)

r!i:-

( •

an expression for the updated moisture at level 1 is given by

r qlf,n_l = 7rql + 7rql-1 --7---I /\O.l_ 1 _ -- //(Trql-1 _)mO'l_ 1 "_

LXffl ] initial \ / final

(19)

Assuming that 7cql-1 is initially negative, we require that its final value is set to zero. Thus

7Cql_lfina_ = O,

_rqlf_naz = _rql + _rql-1---v---J AO.l_ 1 x (20)

ZA(Tl / initial

• i

This process is repeated until the lowest level is reached. If the resulting _rqnlav is negative,

the mass-weighted specific humidity in the lowest model level is simply set to zero.

3 Atmospheric Physics

3.1 Moist Convective Processes

3.1.1 Sub-grid and Large-scale Convection

Sub-grid scale penetrative and shMlow cumulus convection are parameterized using the Relaxed Arakawa Schubert (RAS) scheme of Moorthi and Suarez (1992), which is a simplified

Arakawa Schubert type scheme. The RAS scheme predicts mass fluxes from cloud types

which have different entrainment rates and levels of neutral bouyancy, depending on the

properties of the large scale environment.

The thermodynamic variables that are used in RAS to describe the grid scale vertical profile

are the dry static energy, s = cpT + gz, the moist static energy, h = cpT + gz + Lq, and the

saturation moist static energy, h* = %T + gz + Lq*, where the superscript • refers to the

saturation value.
: i_i_

The conceptualmodelbehindRAS depictsthe cloudas a rising plume,entrainingmass

from the environmentduringascent,anddetrainingmassat thelevelof neutralbuoyancy.

RAS assumesthat the cloudmassflux, _, normalizedby the cloudbasemassflux, is a

linearfunctionof height,expressedas:

O_(z) P ) cp _ )_ or - t?)_

Oz OP" g

where we have used the hydrostatic equation written in the form:

Oz _ % t_

OP _ g

The entrainment parameter, A, characterizes a particular cloud type based on its detrainment level, and is obtained by assuming that the level of detrainment is the level of neutral

buoyancy, ie., the level at which the moist static energy of the cloud, hc, is equal to the saturation moist static energy of the environment, h*. Following Moorthi and Suarez (1992),

A may be written as

hB - h_9 A=

_r p. O(h*D - h)dP,_'

g JPD

where the subscript B refers to cloud base, and the subscript D refers to the detrainment

level.

The convective instability is measured in terms of the cloud work function A, defined as the

rate of change of cumulus kinetic energy. The cloud work function is related to the buoyancy,

or the difference between the moist static energy in the cloud and in the environment:

A= 1+7 [ P'_ dP'_

where 7 is L _-T obtained from the Claussius Clapeyron equation, and the subscript c refers

to the valuePinside the cloud.

To determine the cloud base mass flux, the rate of change of A in time due to dissipation by

the clouds is assumed to approximately balance the rate of change of A due to the generation

by the large scale. This is the quasi-equilibrium assumption, and results in an expression

for mB:

m B -- K

where K is the cloud kernel, defined as the rate of change of the cloud work function

per unit cloud base mass flux, and is obtained by differentiating the expression for A in

time. The rate of change of A due to the generation by the large scale is written as the

difference between the current A(t + At) and its value after the previous convective time

step A(t), divided by the time step. Since the convective parameterization is designed to

nearly neutralize the instability, A(t) is approximated as some (near zero) critical Ac_.it.

10
7,,

/ 7,  ̧

• > :::

The predicted convective mass fluxes are used to solve grid-scale temperature and moisture budget equations to determine the impact of convection on the large scale fields of

temperature (through latent heating and compensating subsidence) and moisture (through

precipitation and detrainment):

O-_tc = aRc-c_g_op mB Os

and

wheree = andP = (p/po).

O-_ttc = n mB r Oh Os

As an approximation to a full interaction between the different allowable cloud types, many

clouds are simulated frequently, each modifying the large scale environment some fraction c_

of the total adjustment. The parameterization thereby "relaxes" the large scale environment

towards neutrality.

In addition to the RAS cumulus convection scheme, the GEOS-1 GCM employs a Kesslertype scheme for the re-evaporation of falling rain (Sud and Molod, 1988), which correspondingly adjusts the impact on the large scale environment by the factor R. The scheme

accounts for the rainfall intensity, the drop size distribution, and the temperature, pressure

and relative humidity of the surrounding air. The moisture deficit in any model layer into

which the rain may re-evaporate is a free parameter.

Due to the increased vertical resolution in the Planetary Boundary Layer (PBL), the lowest

two model layers are averaged to provide the sub-cloud layer for RAS (nominally 50 mb

thick). Each time RAS is invoked (every ten simulated minutes), the possiblity for shallow

convection is checked for the two layers just above cloud base. RAS also randomly chooses

10 other cloud-top levels for the possibility of convection, from just above cloud base to the

model top layer.

Supersaturation or large-scale convection is defined in the GEOS GCM whenever the specific

humidity in any grid-box exceeds its supersaturation value. The large-scale precipitation

scheme rains at supersaturation, and re-evaporates during descent to partially saturate

lower layers in a process that accounts for some simple microphysics.

3.1.2 Cloud Formation

Convective and large-scale cloudiness which is used for cloud-radiative interactions are determined diagnostically as part of the cumulus and large-scale parameterizations. The

convective and large-scale cloud fractions are combined into two separate arrays for use in

11
the shortwaveandlongwaveradiationpackages,onearrayfor randomoverlap(CLRO)and

onearray for maximumoverlap(CLMO) cloudiness.

If convectionoccurswith a cloud-toppressurelessthan 400mb, a cloudfractionequalto

1.0is assignedinto CLMO at the detrainmentlevel. Thesetall convectivetowersaresaid

to bemaximumlyoverlapped,ie. theyareradiativelycorrelatedin the vertical.If shallow

convectionoccurswith a detrainmentlevelwithin 1or 2levelsabovethecloudbase,acloud

fractionequalto 0.5is assignedinto CLRO,againat the detrainmentlevel. Theseshallow

cloudsarenot correlatedin theverticalandaresaidto berandomlyoverlapped.Nocloud

fractionsareprescribedto convectionwhosecloud-topisbetween400mband2levelsabove

cloud-base.

Supersaturationor largescalecloudinessis definedwheneverthe largescaleprecipitation

schemedeterminesthat thegrid boxat anylevelbecomesupersaturated.In orderto ensure

that at any instant the total cloudfractionbe lessthan or equalto 1.0,supersaturation

cloudsareonlyprescribedwhentherearenodeepconvectiveclouds.Undersuchconditions,

the grid box is assumedto be instantaneouslyfully coveredwith randomlyoverlapped

clouds,with a cloudfractionof 1.0beingassignedinto CLRO.

3.2 Radiation

/r. I

The parameterization of radiation includes both shortwave radiation and longwave radiation

in the GEOS-1 GCM. A single model "sponge" layer is inserted into the GEOS-1 GCM above

PTop with a temperature equal to the GCM's first level temperature. This is done in order

to avoid exposing the first GCM level to a zero downward thermal flux condition at PTop.

Radiative fluxes are calculated at each model edge-level in both up and down directions.

The heating rates/cooling rates are then obtained from the vertical divergence of the net

radiative fluxes.

The net flux is

F = F _ - F _

where F is the net flux, F [ is the upward flux and F l is the downward flux.

The heating rate due to the divergence of the radiative flux is given by

OpcpT OF

Ot Oz

or

OT g OF

Ot Cp_rOa

where g is the accelation due to gravity and Cp is the heat capacity of air at constant

pressure.

12
..... !

./

The infrared and solar radiation parameterizations employed in the GEOS-1 GCM follow

closely those described by Harshvardhan et al. (1987). The time tendency for Longwave

Radiation is updated every 3 hours. The time tendency for Shortwave Radiation is updated

once every three hours assuming a normMized incident solar radiation, and subsequently

modified at every model time step by the true incident radiation. For the AMIP simulations,

the GEOS-1 GCM used the AMIP-prescribed solar constant value of 1365 W/m 2 and a C02

mixing ratio of 345 ppm. For the 5-year GEOS-DAS re-analysis, a solar constant value of

1380 W/m 2 was used together with a C02 mixing ratio value of 330 ppm. For the ozone

mixing ratio, monthly mean zonally averaged climatological values specified as a function of

latitude and height (Rosenfield, et al., 1987) are linearly interpolated to the current time.

3.2.1 Shortwave Radiation

The parameterization of heating due to shortwave radiation is an extension of that developed

by Lacis and Hansen (1974). The documentation of the computational aspects of the scheme

are given in Davis (1982). Harshvardhan improved the treatment of the cloud scattering

by allowing it to depend on the solar zenith angle. Two major absorbers, H20 and 03,

are considered in the scheme. Radiative heating by water vapor absorption is modeled in

the 0.7-0.4 pm region. The total absoption of the ozone includes the absoption of visible

radiation by the Chappuis bands (450-750 nm) and the absorption of ultraviolet radiation

by the Hartly (240-280 nm) and Huggins bands (280-360 nm). For a cloudy atmosphere, the

Meadow-Weaver (1980) modified Eddington approximation is adopted to solve the radiative

transfer equation in a scattering medium. The multiple-scattering computation with the

k-distribution functions is used to calculate the cloudy-sky water absorption.

3.2.2 Longwave Radiation

_'• •S ̧ / !

,_ i ̧

The parameterization of longwave radiation employs a wide band model. Four broadband

transmissions have been used in the model. The parameterization of the H20, containing

the two water bands at 15 #m and 9.6 #m, is calculated using the method of Chou (1984)

based on the water vapor line and water vapor e-type absorption approximations. The C02

band employs the scheme of Chou and Peng (1983) which seperates the band wing and

band center scaled paths. The 03 parameterization approach applies the formulation of

Rosenfield et al. (1987), which modifies the line width used by Rodgers (1968) to include

the affects of Doppler broadening.

13
3.2.3 Cloud-Radiation Interaction

The cloud fraction values produced by the Moist Convective Processes are used for cloudradiation interactions as a means for evaluating clear line-of-site probabilities and effective

optical thicknesses.

If we define the time-averaged random and maximum overlapped cloudiness as CLRO and

CLMO respectively, then the probability of clear sky associated with random overlapped

clouds at any level is (1-CLRO) while the probability of clear sky associated with maximum

overlapped clouds at any level is (1-CLMO). The total clear sky probability is given by

(1-CLRO)*(1-CLMO), thus the total cloud fraction at each level may be obtained by 1-(1-

CLRO)*(1-CLMO).

At any given level, we may define the clear line-of-site probability by appropriately accounting for the maximum and random overlap cloudiness. The clear line-of-site probability is

defined to be equal to the product of the clear line-of-site probabilities associated with

random and maximum overlap cloudiness. The clear line-of-site probability C(p, pl) associated with maximum overlap clouds, from the current pressure p to the model top pressure,

pl = Prop, or the model surface pressure, pt = Psurf, is simply 1.0 minus the largest maximum

overlap cloud value along the line-of-site, ie.

1- MAX p' (CLMOp)

Thus, even in the time-averaged sense it is assumed that the maximum overlap clouds are

correlated in the vertical. The clear line-of-site probability associated with random overlap

clouds is defined to be the product of the clear sky probabilities at each level along the

line-of-site, ie.

p!

I'I (1-CLROp)

P

The total cloud fraction at a given level associated with a line- of-site calculation is given

by

: ?

• r'

p!

1 - (1 - MAX;' [C'LMOp]) II (1 - CLROp)

P

The GEOS-1 GCM cloud optical thicknesses are specified based on cloud type and temperature. The "maximum overlap" clouds are assigned an optical thickness of 0.16 km -1

14
i//i')

, i, _ •

• /

• i

for the shortwave radiative feedback. The "random overlap" clouds are assigned an optical

thickness based on an empirical relation between local temperature and optical thickness.

The relation gives a thicker cloud when the temperature is higher and is expressed:

0

(T - 190.66) 2 , 2x10 -6

_ = 6.95x10 -3 * T - 1.82

0.08

for T < 190.66

for 190.66 < T _< 263.16

for 263.16 < T _< 273.38

for 273.38 < T

The relation for the range 190.66 < T _< 263.16 is from Platt and Harshvardahn, 1988. The

optical thickness for the longwave radiative feedback is assumed to be 75% of these values.

The entire Moist Convective Processes Module is called with a frequency of 10 minutes.

The cloud fraction values are time-averaged over the period between Radiation calls (every

3 hours). Therefore, in a time-averaged sense, both random overlap and maximum overlap

cloudiness can exist in a given grid-box.

3.3 Turbulence

The GEOS GCM turbulence parameterization consists of elements which handle vertical

diffusion (Helfand and Labraga, 1988) and surface fluxes of heat, moisture and momentum

(Helfand, et al, 1991, and Helfand and Schubert, 1994). The parameterization is invoked

every 30 minutes, and employs a backward-implicit iterative time scheme with an internal

time step of 5 minutes. The vertical regime is divided into a free atmosphere, a surface

layer, and a viscous sub-layer above the surface roughness elements. The turbulent eddy

fluxes are calculated using a variety of methods depending on the vertical location in the

atmosphere.

Turbulent eddy fluxes of momentum, heat and moisture in the surface layer are calculated

using stability-dependant bulk formulae based on Monin-Obukhov similarity functions. For

an unstable surface layer, the chosen stability functions are the KEYPS function (Panofsky,

1973) for momentum, and its generalization for heat and moisture. The function for heat and

moisture assures non-vanishing heat and moisture fluxes as the wind speed approaches zero.

For a stable surface layer, the stability functions are those of Clarke (1970), slightly modified

for the momemtum flux. The moisture flux also depends on a specified evapotranspiration

coefficient, set to unity over oceans and dependant on the climatological ground wetness over

land. The gradients in the viscous sublayer are based on Yaglom and Kader (1974). The

surface roughness length over oceans is computed as an interpolation between the functions

of Large and Pond (1981) for high winds and of Kondo (1975) for weak winds, and over

land is specified from the climatology of Dorman and Sellers (1989).

15
Abovethe surfacelayer,turbulent fluxesof momentum,heatandmoisturearecalculated

by the Level2.5Mellor-Yamadatypeclosureschemeof ttelfandandLabraga(1988)5which

predictsturbulent kinetic energyand determinesthe eddy transfercoefficientsusedfor

verticaldiffusion.

3.3.1 Planetary Boundary Layer

The Planetary Boundary Layer (PBL) height is diagnosed by the parameterization as the

level at which the turbulent kinetic energy is reduced to a tenth of its surface value, and

generally occurs within the first two to four sigma levels above the surface (sigma equal

0.994 to about sigma equal 0.875).

4 Boundary Conditions

" i_

The conditions at the lower boundary which are not explicitly predicted during execution

must be specified from observed fields. Boundary condition data sets are available at present

in monthly mean increments, at the model's 4° x 5° and 2 ° x 2.5 ° resolutions, for either

climatological or yearly varying conditions (beginning from 1979). Monthly mean values

are interpolated during each model timestep to the current time. Future model versions

will incorporate boundary conditions at higher spatial (1 ° x 1°) and temporal resolutions.

4.1 Orography

Surface geopotential heights are provided from an averaging of the Navy 10 minute by 10

minute dataset supplied by the National Center for Atmospheric Research (NCAR) to the

model's grid resolution. For the 2° x 2.5 ° GEOS-1 GCM, the averaged topography was then

passed through a Lanczos (1966) filter in both dimensions which removes the smallest scales

while inhibiting Gibbs phenomena. Negative values in the topography resulting from the

filtering procedure were not filled.

In one dimension, we may define a cyclic function in x as:

N

a0

f(x) = --f + (ak cos(kx) q- bksin(kx)) (21)

k=l

_ 2_r

where N = _ and IM is the total number of points in the x direction. Defining Ax - i_,

we may define the average of f(x) over a 2Ax region as:

/

16
1 {x+Ax f(x')dx'

f(x) - 2Ax J_-Ax

Using equation (21) in equation (22) and integrating, we may write:

a0 fi[

f(x) = --_ + _ [ak k / x-Ax - bk -_ /x-Ax] k--1

or

(22)

(23)

f(x) = --_a°+ _N sin(kAx)_ (akcos(kx)+ bksin(kx))

k=l

(24)

Thus, the Fourier wave amplitudes are simply modified by the Lanczos filter response function _ This may be compared with an ruth-order Shapiro (1970) filter response kAx "

function, defined as 1 - sinm( k_ shown in Figure 3.

4.2 Land-Sea-Ice Surface Type

The determination of the land or sea category of surface type was made from NCAR's.

10 minute by 10 minute Navy topography dataset, which includes information about the

percentage of water-cover at any point. The data were averaged to the model's 4 ° x 5 ° and

2 ° x 2.5 ° grid resolutions. Any grid-box whose averaged water percentage was >_ 60% was

defined as a water point. Permanent ice or glacier points for the 2° x 2.5 ° grid were identified

from the climatology of Parkinson's satellite-derived fields (NASA SP 459: Arctic Sea Ice,

1973-1976, and NASA SP 489: Antartic Sea Ice, 1973-1976). The Land/Water masks and

topography used for the 4° x 5° and 2 ° x 2.5 ° GEOS-1 GCM grids are shown in Figure 4

and Figure 5, respectively.

4.3 Sea Surface Temperature

Climatological sea surface temperatures for the 4 ° x 5 ° and 2 ° x 2.5 ° grid resolutions are

interpolated from the data of Reynolds, 1988. Yearly varying monthly mean data for the 4 °

x 5 ° grid are from the AMIP specified monthly mean fields (Gates, 1992), and for the 2° x

2,5 ° grid are interpolated from Reynolds' (1988) data. Monthly mean data are interpolated

linearly to the current time at every GCM timestep.

17
1

0.9

0.8

0.7.

0.6,

0.5

0.4

0.3

0.2-

0.1

0

-0.1

Filter Response Functions

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

i : i\i i "Xi. i ! i .... i ..... i

........::....._:.....::.....i i\_ .......X,i' ....................

: ) \; i :_

LORCZOS

2nd-Order Shopiro

4th-Order Shopiro

8th-Order Shapiro

5 1'0 1'5 2'0 2'5 3'0 3'5 4:0

I

: ..... • "1

I

4'5 _o _5 _ _ 7'0

Wavenumber

Figure 3: Comparison between the Lanczos and ruth-order Shapiro filter response functions

for m = 2, 4, and 8.

18
: ii

i _ ;' ,J,

• L _,

i'

i

,) )

z

o

V

eQ.

o

L..

o

Q.

0

I..-

"0

c

o

0

6.

0

"0

e..

0

..J

&

"0

LO

x

V

o

o

lad

o

iiiiii

iiiiii_

iiliiiiiiiii i

iliiiiii_i_i _iiiiiiiiiii_I f ili_i_i

c_

-0o

ILl

.o

Figure 4:GEOS-1 GCM Land/Water Mask and Topography used at 4 ° x 5°resolution.

Light grey shading denotes land surfaces, while dark grey denotes permanent glacier.

•%: •, •!

ii ̧ ̧_,•!:!

i i!__

/ i• •

19
!

2

"v

L--

O

...t

&

ur)

04

x

04

v

W

{.9

Figure 5:GEOS-1 GCM Land/Water Mask and Topography used at 2° x 2.5 ° resolution.

Light grey shading denotes land surfaces, while dark grey denotes permanent glacier.

2O
i ̧:)i

_ _ i_ • '

i_i_

.L

4.4 Sea Surface Roughness

Ocean surface roughness is computed iteratively as a function of wind stress in the Turbulence Scheme.

4.5 Albedo

Monthly mean albedos over land are data of Posey and Clapp (1964), with modifications as

described in Kitzmiller (1979). Over oceans the albedo is specified as 0.07, and over ice (sea

ice or glaciers) 0.80. Monthly mean values are linearly interpolated to the current time.

4.6 Sea Ice

Monthly mean sea ice extent were interpolated to the model's 4 ° x 5° and 2° x 2.5 ° grid resolutions from the AMIP monthly mean fields. Monthly mean data are linearly interpolated

to the current time. Any sea ice is assumed to be three meters thick, and the conduction

through sea ice is accounted for in the Turbulence Parameterization as part of the surface

energy budget.

4.7 Snow Cover

Snow cover is assumed to occur when the ground temperature over the land is less than

freezing in conjunction with climatological albedo values greater than 0.4. Surface conductive characteristics are adjusted accordingly.

4.8 Land Surface

Monthly varying climatological roughness lengths are specified for each land surface vegetation type from Dorman and Sellers (1989). A gridded surface vegetation type designation

is from Sud, et M. (1990). Monthly mean values are interpolated linearly to the model's

current time. The soil wetness climatology and yearly varying fields are from the monthly

estimates of Schemm, et al. (1992) based on the procedure developed by Mintz and Serafini

(1984) using an inverted single layer "bucket" model together with observed fields of surface

air temperature and precipitation. Values are interpolated linearly to the model's current

time. Ground temperature over land is predicted from a surface energy balance equation.

(

21

'_ ,/if

/

i •

Part II

GEOS-1 GCM: Diagnostics

5 Introduction

This section of the documentation describes the Diagnostics Utilities available within the

Goddard Earth Observing Systems (GEOS) General Circulation Model (GCM). In addition

to a description on how to set and extract diagnostic quantities, this document also provides

a comprehensive list of all available diagnostic quantities and a short description of how

they are computed. It should be noted that this document is not intended to be a complete

documentation of the various physical parameterization schemes used in the GEOS GCM,

and the reader should refer to original publications for further insight.

6 Diagnostic Overview

A large selection of model diagnostics is available in the GEOS GCM. At the time of this

writing there are 119 different diagnostic quantities which can be enabled for an experiment. As a matter of philosophy, no diagnostic is enabled as default, thus each user must

specify the exact diagnostic information required for an experiment. This is accomplished

by enabling the specific diagnostic of interest cataloged in the Diagnostic Menu (see Section

7). The Diagnostic Menu is a hard-wired enumeration of diagnostic quantities available

within the GEOS GCM. Diagnostics are internally referred to by their associated number

in the Diagnostic Menu. Once a diagnostic is enabled, the GEOS GCM will continually

increment an array specifically allocated for that diagnostic whenever the associated process

for the diagnostic is computed. Separate arrays are used both for the diagnostic quantity

and its diagnostic counter which records how many times each diagnostic quantity has been

computed.

There are several utilities within the GEOS GCM available to users to enable, disable,

clear, and retrieve model diagnostics, and may be called from any user-supplied application

and/or output routine. The available utilities and the CALL sequences are listed below.

SETDIAG: This subroutine enables a diagnostic from the Diagnostic Menu, meaning that

space is _llocated for the diagnostic and the model routines will increment the diagnostic

value during execution. This routine is useful when called from either user application

routines or user output routines, and is the underlying interface between the user and the

desired diagnostic. The diagnostic is referenced by its diagnostic number from the menu,

and its calling sequence is given by:

PAGE B,L.AI'4IK NOT F MPD

23
where

CALL SETDIAG (NUM)

NUM = Diagnostic number from menu

GETDIAG: This subroutine retrieves the vMue of a model diagnostic. This routine is

particulary useful when called from a user output routine, although it can be called from an

application routine as well. This routine returns the time-averaged value of the diagnostic

by dividing the current accumulated diagnostic value by its corresponding counter. This

routine does not change the value of the diagnostic itself, that is, it does not replace the

diagnostic with its time-average. The calling sequence for this routine is givin by:

where

CALL GETDIAG (LEV,NUM,QTMP,UNDEF)

LEV

NUM

QTMP

UNDEF

= Sigma Level at which the diagnostic is desired

= Diagnostic number from menu

= Time-Averaged Diagnostic Output

= Fill value to be used when diagnostic is undefined

CLRDIAG: This subroutine initializes the values of model diagnostics to zero, and is

particularly useful when called from user output routines to re-initialize diagnostics during

the run. The calling sequence is:

where

CALL CLRDIAG (NUM)

NUM = Diagnostic number from menu

ZAPDIAG: This entry into subroutine SETDIAG disables model diagnostics, meaning

that the diagnostic is no longer available to the user. The memory previously allocated to

the diagnostic is released when ZAPDIAG is invoked. The calling sequence is given by:

where

CALL ZAPDIAG (NUM)

NUM = Diagnostic number from menu

DIAGSIZE: We end this section with a discussion on the manner in which computer

memory is allocated for GEOS diagnostics. All GEOS GCM diagnostic quantities are

24
storedin the singlediagnosticarray QDIAG whichis locatedin the DIAG COMMON,

havingthe form:

COMMON/DIAG/NDIAG_MAX,QDIAG(IM,JM,1)

i

where NDIAG_MAX is an Integer variable which should be set equal to the number of

enabled diagnostics, and QDIAG is a three-dimensional array. The first two-dimensions

of QDIAG correspond to the horizontal dimension of a given diagnostic, while the third

dimension of QDIAG is used to identify specific diagnostic types. In order to minimize the

maximum memory requirement used by the model, the GEOS GCM libraries are compiled

with room for only one horizontal diagnostic array, as shown in the above example. In

order for the User to enable more than 1 two-dimensional diagnostic, the size of the DIAG

COMMON must be expanded to accomodate the desired diagnostics. This can be accomplished by including the COMMON DIAG in the User's Application program, with the

QDIAG array properly dimensioned. Another method, currently used by DAO production

routines, involves the use of the DIAGSIZE subroutine. Using the DIAGSIZE subroutine

is a convenient method to allocate just enough permanent memory during a run to be

used for diagnostic storage. All production DAO applications contain a call to subroutine

DIAGSIZE, which expands the size of the DIAG COMMON to include the total number

of surface and upper-air diagnostics the User has enabled. An example of the DIAGSIZE

subroutine is shown below:

subroutine diagsize

parameter ( IM = 144 )

parameter ( JM = 91 )

parameter ( NLAY = 20 )

parameter ( ndiags = 11 )

parameter ( ndiagu = 4 )

common /diag/ ndiag_max,qdiag(im,jm,ndiags+ndiagu*nlay)

ndiag_max = ndiags + ndiagu*nlay

return

end

25
7 GEOS GCM Diagnostic Menu

N NAME UNITS LEVELS DESCRIPTION

1 UFLUX Newton/m 2 1

2 VFLUX Newton/m 2 1

3 HFLUX Watts/m 2 1

4 EFLUX Watts/rn 2 1

5 QICE Watts/rn 2 1

6 RADLWG Watts/rn 2 1

7 RADSWG Watts/rn 2 1

8 RI deg/day NLAY

9 CT m/sec 1

10 CU m/sec 1

11 ET m2/sec NLAY

12 EU m2/sec NLAY

13 TURBU rn/sec/day NLAY

14 TURBV reseeday NLAY

15 TURBT deg/day NLAY

16 TURBQ g/kg/day NLAY

17 MOISTT deg/day NLAY

18 MOISTQ g/kg/day NLAY

19 RADLW deg/day NLAY

20 RADSW deg/day NLAY

21 PREACC rnm/day 1

22 PRECON ramday 1

23 TUFLUX Newton/m 2 NLAY

24 TVFLUX Newton/m 2 NLAY

25 TTFLUX Watts/m 2 NLAY

26 TQFLUX Watts/m 2 NLAY

27 CN rn / sec 1

28 WINDS m/sec 1

29 DTSRF deg 1

30 TG deg 1

31 TS deg 1

32 DTG deg 1

Surface U-Wind Stress on the atmosphere

Surface V-Wind Stress on the atmosphere

Surface Flux of Sensible Heat

Surface Flux of Latent Heat

Heat Conduction through Sea-Ice

Net upward LW flux at the ground

Net downward SW flux at the ground

Richardson Number

Surface Drag coefficient for T and Q

Surface Drag coefficient for U and V

Diffusivity coefficient for T and Q

Diffusivity coefficient for U and V

U-Momentum Changes due to Turbulence

V-Momentum Changes due to Turbulence

Temperature Changes due to Turbulence

Specific Humidity Changes due to Turbulence

Temperature Changes due to Moist Processes

Specific Humidity Changes due to Moist Processes

Net Longwave heating rate for each level

Net Shortwave heating rate for each level

Total Precipitation

Convective Precipitation

Turbulent Flux of U-Momentum

Turbulent Flux of V-Momentum

Turbulent Flux of Sensible Heat

Turbulent Flux of Latent Heat

Neutral Drag Coefficient

Surface Wind Speed

Air/Surface virtual temperature difference

Ground temperature

Surface air temperature (Adiabatic from lowest model layer)

Ground temperature adjustment

26
N NAME UNITS LEVELS DESCRIPTION

r

i

33 QG g/kg 1

34 QS g/kg 1

35 TGRLW deg 1

36 ST4 Watts/m 2 1

37 OLR Watts/m 2 1

38 OLRCLR Watts/m 2 1

39 LWGCLR Watts/m 2 1

40 LWCLR deg/day NLAY

41 TLW deg NLAY

42 SHLW g/g NLAY

43 OZLW g/g NLAY

44 CLMOLW 0- 1 NLAY

45 CLROLW 0- 1 NLAY

46 CLMOSW 0- 1 NLAY

47 CLROSW 0- 1 NLAY

48 RADSWT Watts/m 2 1

49

50

51

52 EVAP ramday 1

53 DPDT mm/day 1

54 ANALU m/sec/sec NLAY

55 ANALV m/sec/sec NLAY

56 ANALT mbdcg/mb_/sec NLAY

57 ANALQ mb,g/g/sec NLAY

58 OMEGA mb/day NLAY

27

Ground specific humidity

Saturation surface specific humidity

Instantaneous ground temperature used as input to the Longwave radiation subroutine

Upward Longwave flux at the ground (aT 4)

Net upward Longwave flux at the top of the

model

Net upward clearsky Longwave flux at the top

of the model

Net upward clearsky Longwave flux at the

ground

Net clearsky Longwave heating rate for each

level

Instantaneous temperature used as input to

the Longwave radiation subroutine

Instantaneous specific humidity used as input

to the Longwave radiation subroutine

Instantaneous ozone used as input to the

Longwave radiation subroutine

Maximum overlap cloud fraction used in the

Longwave radiation subroutine

Random overlap cloud fraction used in the

Longwave radiation subroutine

Maximum overlap cloud fraction used in the

Shortwave radiation subroutine

Random overlap cloud fraction used in the

Shortwave radiation subroutine

Incident Shortwave radiation at the top of the

atmosphere

Disabled

Disabled

Disabled

Surface evaporation

Total surface pressure tendency

Analysis u - Wind increment

Analysis v - Wind increment

Analysis _r_ increment, t_ =

Analysis _rq increment

Vertical (Omega) velocity
N NAME UNITS LEVELS DESCRIPTION

DUDT reseeday NLAY

DVDT m/sec/day NLAY

DTDT deg/day NLAY

DQDT g/kg/day NLAY

VORT sec -1 NLAY

59

60

61

62

63

64

65

66

67 USTAR

68 Z0

69 FRQTRB

70 PBL

71 SWCLR

72 OSR

73 OSRCLR

m/sec 1

m 1

0- 1 NLAY-1

mb 1

deg/day NLAY

Watts/m 2 1

Watts/m 2 1

74 CLDMAS kgm/sec NLAY

75 UAVE m/sec NLAY

76 VAVE m/sec NLAY

77 TAVE deg NLAY

78 QAVE g/g NLAY

79

80 PAVE mb 1

81 QQAVE (m/sec) 2 NLAY

82 SWGCLR Watts/m 2 1

83 ANALP mb/sec 1

84 SDIAG1 1

85 SDIAG2 1

86 UDIAG1 NLAY

87 UDIAG2 NLAY

88 DIABU m/sec/day NLAY

89 DIABV m/sec/day NLAY

90 DIABT deg/day NLAY

91 .DIABQ g/kg/day NLAY

Total U-Wind tendency

Total V-Wind tendency

Total Temperature tendency

Total Specific Humidity tendency

Relative Vorticity

Disabled

Disabled

Disabled

Surface USTAR wind

Surface roughness

Frequency of Turbulence

Planetary Boundary Layer depth

Net clearsky Shortwave heating rate for each

level

Net downward Shortwave flux at the top of

the model

Net downward clearsky Shortwave flux at the

top of the model

Convective cloud mass flux

Time-averaged u - Wind

Time-averaged

Time-averaged

Time-averaged

Disabled

Time-averaged

Time-averaged

Net downward

ground

v - Wind

Temperature

Specific Humidity

Psur] -- Prop

T u r bul ent K inet ic E ner g y

clearsky Shortwave flux at the

Analysis SurfacePressure increment

User-Defined Surface Diagnostic-1

User-Defined Surface Diagnostic-2

User-Defined Upper-Air Diagnostic-1

User-Defined Upper-Air Diagnostic-2

Total Diabatic forcing on u - Wind

Total Diabatic forcing on v - Wind

Total Diabatic forcing on Temperature

Total Diabatic forcing on Specific Humidity

28
N NAME UNITS LEVELS DESCRIPTION

92

93

94

95

96

97

98

99

100

101

102

103

104

105

I06

107

108

109

110

111

112

113

114

115

116

117

118

119

VINTUQ m/sec.g/kg 1

VINTVQ m/sec.g/kg 1

VINTUT m/sec, deg 1

VINTVT m/sec, deg 1

CLDFRC 0- 1 1

QINT gm/cm 2 1

U2M m / sec 1

V2M m/sec 1

T2M deg 1

Q2M g/kg 1

U10M m/sec I

V10M m/sec 1

T10M deg 1

Q10M g/kg 1

DTRAIN kg • m/sec NLAY

QFILL g/kg/day NLAY

VINTQANA ramday 1

VINTQFIL mm/day 1

Disabled

Disabled

Disabled

Disabled

Disabled

Disabled

Disabled

Disabled

Disabled

Disabled

Vertically integrated u q

Vertically integrated v q

Vertically integrated u T

Vertically integrated v T

Total Cloud Fraction

Precipitable water

U-Wind at 2 meters

V-Wind at 2 meters

Temperature at 2 meters

Specific Humidity at 2 meters

U-Wind at 10 meters

V-Wind at 10 meters

Temperature at 10 meters

Specific Humidity at 10 meters

Convective Cloud Detrainment

Filling of negative specific humidity

Vertically integrated analysis specific humidity increment

Vertically integrated QFILL

%

/

29
8 Diagnostic Description

In this section we list and describe the diagnostic quantities available within the GEOS

GCM. The diagnostics are listed in the order that they appear in the Diagnostic Menu,

Section 7. In all cases, each diagnostic as currently archived on DAO production output

datasets is time-averaged over its diagnostic output fi'equency:

1 t:TTOT

DIAGNOSTIC-TTOT _ diag(t)

t----1

where TTOT NQDIAG NQDIAG is the output frequency of the diagnositc, and At is

-- At '

the timestep over which the diagnostic is updated. For further information on how to set

the diagnostic output frequency NQDIAG, please see Part III, A User's Guide. Also, note

that several Diagnostic Index locations are currently disabled due to changes which have

occured from previous GCM versions.

A word of caution: The Arakawa C-grid stencil used in the GEOS GCM employs a

"stretched" coordinate near the poles due to conservation requirements. As a result, a

well-defined vorticity is is located at the poles, while the first mass point is located one

A¢ away from the poles. Since all Physical (diabatic) processes of the model are computed

at mass points, all Physics diagnostics are only truely defined at latitudes -90 + A¢ to

90- A¢. While an attempt was made to "diagnose" some pole value diagnostics through

interpolation over the poles, others have been left unfilled.

1) UFLUX Surface Zonal Wind Stress on the Atmosphere (Newton/m 2)

The zonal wind stress is the turbulent flux of zonal momentum fi'om the surface. See section

3.3 for a description of the surface layer parameterization.

UFLUX -pCDWsu where : CD = C_

where p = the atmospheric density at the surface, CD is the surface drag coefficient, Cu is the

dimensionless surface exchange coefficient for momentum (see diagnostic number 10), Ws

is the magnitude of the surface layer wind, and u is the zonal wind in the lowest model layer.

2) VFLUX Surface Meridional Wind Stress on the Atmosphere (Newton/m 2)

The meridional wind stress is the turbulent flux of meridional momentum from the surface.

See section 3.3 for a description of the surface layer parameterization.

VFLUX = -pCDWsV where : CD = C2u

3O
wherep = the atmospheric density at the surface, CD is the surface drag coefficient, Cu is

the dimensionless surface exchange coefficient for momentum (see diagnostic number 10),

Ws is the magnitude of the surface layer wind, and v is the meridional wind in the lowest

model layer.

3) HFLUX Surface Flux of Sensible Heat (Watts/m 2)

The turbulent flux of sensible heat from the surface to the atmosphere is a function of the

gradient of virtual potential temperature and the eddy exchange coefficient:

HFLUX = PgpCpCHWs(Osurface - ONLAY) where : CH = CuCt

where p = the atmospheric density at the surface, cp is the specific heat of air, CH is the

dimensionless surface heat transfer coefficient, W_ is the magnitude of the surface layer

wind, Cu is the dimensionless surface exchange coefficient for momentum (see diagnostic

number 10), Ct is the dimensionless surface exchange coefficient for heat and moisture (see

diagnostic number 9), and _ is the potential temperature at the surface and at the bottom

a-level.

4) EFLUX Surface Flux of Latent Heat (Watts/m 2)

The turbulent flux of latent heat from the surface to the atmosphere is a function of the

gradient of moisture, the potential evapotranspiration fraction and the eddy exchange coefficient:

EFLUX = fl/_LCHWs(qsurface - qNLAY ) where : CH = CuCt

where p = the atmospheric density at the surface,/3 is the fraction of the potential evapotranspiration actually evaporated, L is the latent heat of evaporation, CH is the dimensionless surface heat transfer coefficient, W_ is the magnitude of the surface layer wind, C_ is

the dimensionless surface exchange coefficient for momentum (see diagnostic number 10),

Ct is the dimensionless surface exchange coefficient for heat and moisture (see diagnostic

number 9), and qsurface and qNLAY are the specific humidity at the surface and at the

bottom a-level, respectively.

5) QICE Heat Conduction Through Sea Ice (Watts/m 2)

Over sea ice there is an additional source of energy at the surface due to the heat conduction

from the relatively warm ocean through the sea ice. The heat conduction through sea ice

represents an additional energy source term for the ground temperature equation.

Cti

QICE = -7}-_(Ti - Tg)

_r/i

31
whereCti is the thermal conductivity of ice, Hi is the ice thickness, assumed to be 3 m

where sea ice is present, Ti is 273 degrees Kelvin, and Tg is the temperature of the sea ice.

NOTE: QICE is not available through model version 5.3, but is available in subsequent

versions.

6) RADLWG Net upward Longwave Flux at the surface (Watts/m 2)

_Net RADLWG = *LW, nlay+l

: I[ _LW,nlayT1 -- F_W, nlay+l

where nlay+l indicates the lowest model edge-level, or p = p_r/. FtLW is the upward Longwave flux and F_LW is the downward Longwave flux.

7) RADSWG Net downard shortwave Flux at the surface (Watts/m 2)

_Net RADSWG = l sw,_t_v+l

--_ _7 _SW,nlay__ 1 __ F_W, T nlayq_ 1

where nlay+l indicates the lowest model edge-level, or p = p_7.f. F_w is the downward

Shortwave flux and F_w is the upward Shortwave flux.

8) R_II Richardson Number (dimensionless)

The non-dimensional stability indicator is the ratio of the buoyancy to the shear:

RI=

g _ c _ °P_ 0-70z P Oz

(_au)2 + (Ova2 _zJ - Eo.u'_2 _azJ + (OzJ Ov'_2

, "!

i, _

• _ i  ̧.

where we used the hydrostatic equation:

0¢

Op _ - cpOo

Negative values indicate unstable buoyancy AND shear, small positive values (< 0.4)

indicate dominantly unstable shear, and large positive values indicate dominantly stable

stratification.

32
9) CT Surface Exchange Coefficient for Temperature and Moisture (dimensionless)

The surface exchange coefficient is obtained from the similarity functions for the stability

dependant flux profile relationships:

CT =_(w--0,) (w,qi) k

u, O u, q (¢h+ Cg)

where _bh is the surface layer non-dimensional temperature change and Cg is the viscous

sublayer non-dimensional temperature or moisture change:

and 0.2!(hou* )1/2 Cg = vl/2 - ho_j u._j

and: ho = 30zo with a maximum value over land of 0.01

Ch is the similarity function of _, which expresses the stability dependance of the temperature and moisture gradients, specified differently for stable and unstable layers according

to Helfand and Schubert, 1993. k is the Von Karman constant, _ is the non-dimensional

stability parameter, Pr is the Prandtl number for Mr, v is the molecular viscosity, z0 is the

surface roughness length, u. is the surface stress velocity (see diagnostic number 67), and

the subscript ref refers to a reference value.

10) CU Surface Exchange Coefficient for Momentum (dimensionless)

The surface exchange coefficient is obtained from the similarity functions for the stability

dependant flux profile relationships:

u. k

CU- -

w_ ¢_

where '_bmis the surface layer non-dimensional wind shear:

Cm is the similarity function of _', which expresses the stability dependance of the temperature and moisture gradients, specified differently for stable and unstable layers according

to Helfand and Schubert, 1993. k is the Von Karman constant, _ is the non-dimensional

stability parameter, u. is the surface stress velocity (see diagnostic number 67), and W_ is

the magnitude of the surface layer wind.

11) ET Diffusivity Coefficient for Temperature and Moisture (m2/sec)

..... 33
In the level 2.5 version of the Mellor-Yamada (1974) hierarchy, the turbulent heat or moisture flux for the atmosphere above the surface layer can be expressed as a turbulent diffusion

coefficient Kh times the negative of the gradient of potential temperature or moisture. In

the Helfand and Labraga (1988) adaptation of this closure, Kh takes the form:

(w'8_v) _ { q_SH(GM, GH) for decaying turbulence ET = Kh -- _ - _ g SH(GM_, GH_) for growing turbulence

where q is the turbulent velocity, or x/2* turbulent kinetic energy, qe is the turbulence

velocity derived from the more simple level 2.0 model, which describes equilibrium turbulence, _ is the master length scale related to the layer depth, SH is a function of GH and

GM, the dimensionless buoyancy and wind shear parameters, respectively, or a function of

GH_ and GM_, the equilibrium dimensionless buoyancy and wind shear parameters. Both

GH and GM, and their equilibrium values GH_ and GM_, are functions of the Richardson

number.

For the detailed equations and derivations of the modified level 2.5 closure scheme, see

Helfand and Labraga, 1988.

In the surface layer, ET is the exchange coefficient for heat and moisture, in units of m/see,

given by:

ETNLAY = Ct * u. = CHWs

where Ct is the dimensionless exchange coefficient for heat and moisture from the surface

layer similarity functions (see diagnostic number 9), u. is the surface friction velocity (see

diagnostic number 67), CH is the heat transfer coefficient, and W_ is the magnitude of the

surface layer wind.

12) EU Diffusivity Coefficient for Momentum (m2/sec)

In the level 2.5 version of the Mellor-Yamada (1974) hierarchy, the turbulent heat momentum flux for the atmosphere above the surface layer can be expressed as a turbulent diffusion

coefficient /(m times the negative of the gradient of the u-wind. In the Helfand and Labraga

(1988) adaptation of this closure, K,_ takes the form:

EU = Km (u'w 0 { q_.SM(GM, GH) for decaying turbulence -- O__U_U= q2---_SM(GM_,GH_ ) for growing turbulence

Oz qe

where q is the turbulent velocity, or x/2* turbulent kinetic energy, qe is the turbulence

velocity derived from the more simple level 2.0 model, which describes equilibrium turbulence, _ is the master length scale related to the layer depth, SM is a function of GH and

GM, the dimensionless buoyancy and wind shear parameters, respectively, or a function of

GH. and GM_, the equilibrium dimensionless buoyancy and wind shear parameters. Both

%,

34

• : : ,: rl • •
GH and GM, and their equilibrium values GH_ and GM_, are functions of the Richardson

number.

For the detailed equations and derivations of the modified level 2.5 closure scheme, see

Helfand and Labraga, 1988.

In the surface layer, EU is the exchange coefficient for momentum, in units of m/sec, given

by:

EUNLAY = Cu * u. = CoWs

where Cu is the dimensionless exchange coefficient for momentum from the surface layer

similarity functions (see diagnostic number 10), u. is the surface friction velocity (see diagnostic number 67), Co is the surface drag coefficient, and W_ is the magnitude of the

surface layer wind.

13) TURBU Zonal U-Momentum changes due to Turbulence (m/sec/day)

The tendency of U-Momentum due to turbulence is written:

TURBU = Ou0-7, rb= 0 (-u'w') = 0 ._. Ou

The Helfand and Labraga level 2.5 scheme models the turbulent flux of u-momentum in

terms of Kin, and the equation has the form of a diffusion equation.

14) TURBV Meridional V-Momentum changes due to Turbulence (m/sec/day)

The tendency of V-Momentum due to turbulence is written:

TURBV = Ov 0 0 . Ov -O-tt_,,.b= (--v'w') = _z (hm_z)

The Helfand and Labraga level 2.5 scheme models the turbulent flux of v-momentum in

terms of K,_, and the equation has the form of a diffusion equation.

15) TURBT Temperature changes due to Turbulence (deg/day)

The tendency of temperature due to turbulence is written:

TURBT- OTOt _ p,_O0O---tturb = p,_ O (_w,O, )

= o (i(hOo ) -5-;

The Helfand and Labraga level 2.5 scheme models the turbulent flux of temperature in

terms of Kh, and the equation has the form of a diffusion equation.

35
16) TURBQ Specific Humidity changes due to Turbulence (g/kg/day)

The tendency of specific humidity due to turbulence is written:

Oq 0 _ 0 . Oq

TURBQ = -_turb = _zz(--w q) = -_z(Kh-_z)

The Helfand and Labraga level 2.5 scheme models the turbulent flux of temperature in

terms of Kh, and the equation has the form of a diffusion equation.

17) MOISTT Temperature Changes Due to Moist Processes (deg/day)

where:

and

OT

MOISTT = -_-c

= R_ a rs and

c • i

08

Fs = g_p

+

The subscript c refers to convective processes, while the subscript Is refers to large scale

precipitation processes, or supersaturation rain. The summation refers to contributions

from each cloud type called by RAS. The dry static energy is given as s, the convective

cloud base mass flux is given as roB, and the cloud entrainment is given as y, which are

explicitly defined in Section 3.1.1, the description of the convective parameterization. The

fractional adjustment, or relaxation parameter, for each cloud type is given as a, while R

is the rain re-evaporation adjustment.

18) MOISTQ Specific Humidity Changes Due to Moist Processes (g/kg/day)

where:

and

MOISTQ = 0_tt c+ 0_tt ts

O-_ttc= RE (a_-_ -(Fh- Fs)) and Oq =(q._q) i i Ot is

Os Oh

Fs = g_-z- and Ph = g71

ap Op

The subscript c refers to convective processes, while the subscript Is refers to large scale

precipitation processes, or supersaturation rain. The summation refers to contributions

36
,L

ii I:

.... ii_

from each cloud type called by RAS. The dry static energy is given as 8, the moist static

energy is given as h, the convective cloud base mass flux is given as taB, and the cloud

entrainment is given as y, which are explicitly defined in Section 3.1.1, the description of

the convective parameterization. The fractional adjustment, or relaxation parameter, for

each cloud type is given as (_, while R is the rain re-evaporation adjustment.

19) RADLW Heating Rate due to Longwave Radiation (deg/day)

The net longwave heating rate is calculated as the vertical divergence of the net terrestrial

radiative fluxes. Both the clear-sky and cloudy-sky longwave fluxes are computed within the

longwave routine. The subroutine calculates the clear-sky flux, F_ rsky, first. For a given

cloud fraction, the clear line-of-sight probability C(p, pt) is computed from the current level

pressure p to the model top pressure, pJ = prop, and the model surface pressure, pJ = Ps_r],

for the upward and downward radiative fluxes. (see Section 3.2.3). The cloudy-sky flux is

then obtained as:

FLW = C(p,p') "_LW_dearsky

Finally, the net longwave heating rate is calculated as the vertical divergence of the net

terrestrial radiative fluxes:

OpcpT _ 0 _NET

Ot -----_Z _'LW '

or

RADLW g 0 _NET

-- go. r LW • Cp?r

where g is the accelation due to gravity, % is the heat capacity of air at constant pressure,

and

ELNETW -_ FTLW FILW

J

20) RADSW Heating Rate due to Shortwave Radiation (deg/day)

The net Shortwave heating rate is calculated as the vertical divergence of the net solar

radiative fluxes. The clear-sky and cloudy-sky shortwave fluxes are calculated separately.

For the clear-sky case, the shortwave fluxes and heating rates are computed with both

CLMO (maximum overlap cloud fraction) and CLRO (random overlap cloud fraction) set to

zero (see Section 3.2.3). The shortwave routine is then called a second time, for the cloudysky case, with the true time-averaged cloud fractions CLMO and CLRO being used. in all

cases, a normalized incident shortwave flux is used as input at the top of the atmosphere.

37
The heatingrate dueto ShortwaveRadiationundercloudyskiesis definedas:

OpcpT 0 NET

Ot -- _zF(Cloudy)sw . RADSWT,

or

RADSW - g _aF(cloudy)NET. RADSWT.

CpT_

where g is the accelation due to gravity, Cp is the heat capacity of air at constant pressure, RADSWT is the true incident shortwave radiation at the top of the atmosphere (See

Diagnostic #48), and

21) PREACC Total (Large-scale + Convective) Accumulated Precipltion (ramday)

For a change in specific humidity due to moist processes, Aqmoist, the vertical integral or

total precipitable amount is given by:

PREACC :i,o pAqmoistdz = -- /,o,Aqraoist ? = /01 - Aqmoistda

asurS gsur/ g

A precipitation rate is defined as the vertically integrated moisture adjustment per Moist

Processes time step, scaled to mm/day.

22) PRECON Convective Precipition (ramday)

For a change in specific humidity due to sub-grid scale cumulus convective processes, Aqc_m,

the vertical integral or total precipitable amount is given by:

PRECON = /top pAqcumdz = /top dp TCjfol - Aqcum-- = - Aq_u_da

asurf a_urf g g

A precipitation rate is defined as the vertically integrated moisture adjustment per Moist

Processes time step, scaled to ramday.

38
23) TUFLUX Turbulent Flux of U-Momentum (Newton/m 2)

The turbulent flux of u-momentum is calculated for diagnostic purpose8 only from the

eddy coefficient for momentum:

TUFLUX = p(u'w 7) = p(-K._z )

where p is the air density, and Km is the eddy coefficient.

24) TVFLUX Turbulent Flux of V-Momentum (Newton/m 2)

The turbulent flux of v-momentum is calculated for diagnostic purpose8 only from the

eddy coefficient for momentum:

TVFLUX = p(v'w') = p(-Km_z )

where p is the air density, and Km is the eddy coefficient.

25) TTFLUX Turbulent Flux of Sensible Heat (Watts/m 2)

The turbulent flux of sensible heat is calculated for diagnostic purpose8 only from the

eddy coefficient for heat and moisture:

TTFLUX = cppP_(w'O ') = cppP_(-Kh _ )

where p is the air density, and Kh is the eddy coefficient.

26) TQFLUX Turbulent Flux of Latent Heat (Watts/m 2)

The turbulent flux of latent heat is calculated for diagnostic purposes only from the eddy

coefficient for heat and moisture:

TQFLUX = Lp(w'q') = Lp(-I(h_z )

where p is the air density, and Kh is the eddy coefficient.

39
27) CN Neutral Drag Coefficient (dimensionless)

Thedragcoefficientfor momentumobtainedby assuminga neutrally stablesurfacelayer:

k

CN= 1-.(1⁄4--U

where k is the Von Karman constant, h is the height of the surface layer, and z0 is the

surface roughness.

NOTE: CN is not available through model version 5.3, but is available in subsequent versions.

28) WINDS Surface Wind Speed (meter/sec)

The surface wind speed is calculated for the last internal turbulence time step:

WINDS • = CU2NLAY -_- VNLAY 2

where the subscript NLAY refers to the lowest model a - level.

29) DTSRF Air/Surface Virtual Temperature Difference (deg K)

The air/surface virtual temperature difference measures the stability of the surface layer:

DTSRF = (t_vNLAY +I -- OVNLAY )Ps_ur]

where

Tg

OVNLAY+ 1 -- lp_r f

(1+.609qNLAY+I) and qNLAY+I = qNLAY +_(q (Tg, Ps)--qNLAY)

fl is the surface potential evapotranspiration coefficient (_ = 1 over oceans), q*(Tg, Ps) is the

saturation specific humidity at the ground temperature and surface pressure, level NLAY

refers to the lowest model level and level NLAY + i refers to the surface.

30) TG Ground Temperature (deg K)

The ground temperature equation is solved as part of the turbulence package using a backward implicit time differencing scheme:

c OTg

TGisobtainedfrom: g_ = Rs_ - Rl_ + Qice - H - LE

4O
whereRs_ is the net surface downward shortwave radiative flux, Rl_ is the net surface

upward longwave radiative flux, Qice is the heat conduction through sea ice, H is the upward

sensible heat flux, LE is the upward latent heat flux, and Cg is the total heat capacity of

the ground. Cg is obtained by solving a heat diffusion equation for the penetration of the

diurnal cycle into the ground (Blackadar, 1977), and is given by:

Cg A/_-C-_s I(0.386 + 0.536W + 0.15W2)2x10 -386400"

= •

Here, the thermal conductivity, _, is equal to 2x10 -3 t_y__ o;__,s_c the angular velocity of the

earth, w, is written as 86400 sec/day divided by 27r radians/day, and the expression for

C8, the heat capacity per unit volume at the surface, is a function of the ground wetness, W.

31) TS Surface Temperature (deg K)

The surface temperature estimate is made by assuming that the model's lowest layer is

well-mixed, and therefore that 8 is constant in that layer. The surface temperature is

therefore:

TS = _NLAY Psur]

32) DTG Surface Temperature Adjustment (deg K)

The change in surface temperature from one turbulence time step to the next, solved using

the Ground Temperature Equation (see diagnostic number 30) is calculated:

DTG = Tg n - Tg _-1

where superscript n refers to the new, updated time level, and the superscript n - 1 refers

to the value at the previous turbulence time level.

33) QG Ground Specific Humidity (g/kg)

The ground specific humidity is obtained by interpolating between the specific humidity at

the lowest model level and the specific humidity of a saturated ground. The interpolation

is performed using the potential evapotranspiration function:

QG = qNLAY+I = qNLAY + t_(q*(Ta, Ps) - qNLAY)

where j3 is the surface potential evapotranspiration coefficient (f_ = 1 over oceans), and

q*T( g, Ps) is the saturation specific humidity at the ground temperature and surface pres-

/

41
sure.

34) qS Saturation Surface Specific Humidity (g/kg)

The surface saturation specific humidity is the saturation specific humidity at the ground

temprature and surface pressure:

qs = q*(Ta, P,)

35) TGRLW Instantaneous ground temperature used as input to the Longwave

radiation subroutine (deg)

TGRLW = Tg(A, ¢, n)

where Tg is the model ground temperature at the current time step n.

36) ST4 Upward Longwave flux at the surface (Watts/m 2)

ST4 = aT 4

where a is the Stefan-Boltzmann constant and T is the temperature.

37) OLR Net upward Longwave flux at p = Ptop (Watts m2)

L.N ET OLR = "LW, top

where top indicates the top of the first model layer. In the GEOS-1 GCM, prop

= 10 mb.

38) OLRCLR Net upward clearsky Longwave flux at p =Ptop (Watts m2)

NET

OLRCLR = F(clearskY)LW, top

where top indicates the top of the first model layer. In the GEOS-1 GCM, prop = 10 rob.

39) LWGCLR Net upward clearsky Longwave flux at the surface (Watts/m 2)

Net

LWGCLR = F(clearsky)gw, nlay+l

= F(clearskY)_W, nlay+l- F(clearskY)_LW, nlay+l

42
wherenlay+l indicatesthe lowestmodeledge-level,or p = Psur]. F(clearsky)TLW is the

upward clearsky Longwave flux and the F(clearskY)lLW is the downward clearsky Longwave

flUX.

40) LWCLR Heating Rate due to Clearsky Longwave Radiation (deg/day)

The net longwave heating rate is calculated as the vertical divergence of the net terrestrial

radiative fluxes. Both the clear-sky and cloudy-sky longwave fluxes are computed within the

longwave routine. The subroutine calculates the clear-sky flux, F_ arsky, first. For a given

cloud fraction, the clear line-of-sight probability C(p, p_) is computed from the current level

pressure p to the model top pressure, p_ = Prop, and the model surface pressure, p_ = Psur.f,

for the upward and downward radiative fluxes. (see Section 3.2.3). The cloudy-sky flux is

then obtained as:

FLW = C(p,p') "ILW_d_arsky

Thus, LWCLR is defined as the net longwave heating rate due to the vertical divergence

of the clear-sky longwave radiative flux:

OpcpT Ft

= T,

Ot clearsky UZ

or

LWCLR - g 0

cpTrOa F(clearsky)NEwT"

where g is the accelation due to gravity, Cp is the heat capacity of air at constant pressure,

and

41) TLW Instantaneous temperature used as input to the Longwave radiation

subroutine (deg)

TLW = T(A, ¢, a, n)

where T is the model temperature at the current time step n.

, L,./

, _/i.,i}_(

• : v

42) SHLW Instantaneous specific humidity used as input to the Longwave radiation subroutine (kg/kg)

SHLW = q(A, ¢, a, n)

43
whereq is the model specific humidity at the current time step n.

43) OZLW Instantaneous ozone used as input to the Longwave radiation subroutine (kg/kg)

OZLW = OZ(A, ¢, a, n)

where OZ is the interpolated ozone data set from the climatological monthly mean zonally

averaged ozone data set.

44) CLMOLW Maximum Overlap cloud fraction used in LW Radiation (0 - 1)

CLMOLW is the time-averaged maximum overlap cloud fraction that has been filled by the

Relaxed Arakawa/Schubert Convection scheme and will be used in the Longwave Radiation

algorithm. These are convective clouds whose radiative characteristics are assumed to be

correlated in the vertical. For a complete description of cloud/radiative interactions, see

Section 3.2.3.

CLMOLW = CLMORAs, LW(_, ¢, a)

Note, in Version 1 of the GEOS GCM, shortwave and longwave cloud fields are computed

identically, ie.

CLMOLW = CLMOSW

45) CLROLW Random Overlap cloud fraction used in LW Radiation (0- 1)

CLROLW is the time-averaged random overlap cloud fraction that has been filled by

the Relaxed Arakawa/Schubert and Large-scale Convection schemes and will be used in the

Longwave Radiation algorithm. These are convective and large-scale clouds whose radiative

characteristics are not assumed to be correlated in the vertical. For a complete description

of cloud/radiative interactions, see Section 3.2.3.

CLROLW = C LRORAS, LargeScale,LW(._, ¢, 5r)

Note, in Version 1 of the GEOS GCM, shortwave and longwave cloud fields are computed

identically, ie.

CLROLW = CLROSW

46) CLMOSW Maximum Overlap cloud fraction used in SW Radiation (0 - l)

44
CLMOSW isthetime-averagedmaximumoverlapcloudfractionthat hasbeenfilledbythe:

RelaxedArakawa/SchubertConvectionschemeandwill beusedin theShortwaveRadiation

algorithm. Theseareconvectivecloudswhoseradiativecharacteristicsareassumedto be

correlatedin the vertical. For a completedescriptionof cloud/radiativeinteractions,see

Section3.2.3.

CLMOSW = C LMORAs,SW( A, ¢, a)

Note, in Version 1 of the GEOS GCM, shortwave and longwave cloud fields are computed

identically, ie.

CLMOSW = CLMOLW

47) CLROSW Random Overlap cloud fraction used in SW Radiation (0 - 1)

CLROSW is the time-averaged random overlap cloud fraction that has been filled by

the Relaxed Arakawa/Schubert and Large-scale Convection schemes and will be used in the

Shortwave Radiation algorithm. These are convective and large-scale clouds whose radiative

characteristics are not assumed to be correlated in the vertical. For a complete description

of cloud/radiative interactions, see Section 3.2.3.

CLROSW = C LRORAS, LargeScale,SW( A, ¢, a)

Note, in Version 1 of the GEOS GCM, shortwave and longwave cloud fields are computed

identically, ie.

CLROSW = CLROLW

,j /.

/] :!:?i

48) RADSWT Incident Shortwave radiation at the top of the atmosphere (Watts/m 2)

So

RADSWT = -- • cos¢z

R2a

where So, is the extra-terrestial solar contant, R_ is the earth-sun distance in Astronomical

Units, and cOS¢z is the cosine of the zenith angle. It should be noted that RADSWT,

as well as OSR and OSRCLR, are calculated at the top of the atmosphere (p=0 mb).

However, the OLR and OLRCLR diagnostics are currently calculated at p = Prop (10 mb

for the 20-level GEOS-I GCM).

49) Disabled

45
50) Disabled

51) Disabled

52) EVAP Surface Evaporation (ram/day)

The surface evaporation is a function of the gradient of moisture, the potential evapotranspiration fraction and the eddy exchange coefficient:

EVAP = p/_Kh(qsurface - qNLAY )

where p = the atmospheric density at the surface,/3 is the fraction of the potential evapotranspiration actually evaporated (fl = 1 over oceans), Kh is the turbulent eddy exchange

coefficient for heat and moisture at the surface in m/sec and qsurface and qNLAY are

the specific humidity at the surface (see diagnostic number 34) and at the bottom a-level,

respectively.

53) DPDT Total Surface Pressure Tendency (robday)

DPDT is the total time-tendency of the surface pressure due to Hydrodynamic and Analysis

forcing. There are no surface pressure changes due to physical diabatic processes, such as

changes in total precipitable water amounts.

07r c97r

DPDT Ot Dynamics Ot Analysis

where 7r ps_rf --Ptop. In forecast simulations a_ = 0. _" -D'{Analysis

54) ANALU Analysis Zonal U-Wind Tendency (._1_1_)

ANALU is the Zonal U-Wind tendency created during the assimilation cycle which is used

in the Incremental Analysis Updating (IAU) procedure.

ANALU = (uA()L¢,a,t) - uFG() , ¢, a, t)) /AtAna ysis

where U A corresponds to the u-wind after the analysis, u FG corresponds to the u-wind first

guess, and AtAnalysis is the time-step between analyses (6 hours). Note, ANALU is stored

on the staggered C-grid.

55) ANALV Analysis Meridional V-Wind Tendency (m/sec/sec)

46
ANALV is the MeridionalV-Wind tendencycreatedduring the assimilationcyclewhich

is usedin the IncrementalAnalysisUpdating(IAU) procedure.

ANALV = (vA(A,¢,a,t) - vFG()%¢,a,t)) IAtAnatysis

where v A corresponds to the v-wind after the analysis, vFa corresponds to the v-wind first

guess, and AtAnalysis is the time-step between analyses (6 hours). Note, ANALV is stored

on the staggered C-grid.

56) ANALT Analysis Temperature Tendency (mbdeg/mb'_/sec)

ANALT is the pressure weighted Potential Temperature tendency created during the assimilation cycle which is used in the Incremental Analysis Updating (IAU) procedure.

ANALT = ((Tro)A(_,¢,o',t) - (rrO)FG(A,¢,a,t))/AtAnalysis

where (r_0) A corresponds to the pressure weighted Potential Temperature after the analysis,

(_rO)FG corresponds to the pressure weighted Potential Temperature first guess, _ = T/p _,

and AtAnatysis is the time-step between analyses (6 hours).

57) ANALQ Analysis Temperature Tendency (mb. g/g/sec)

ANALQ is the pressure weighted Specific Humidity tendency created during the assimilation cycle which is used in the Incremental Analysis Updating (IAU) procedure.

ANALQ = ((Trq)A()%¢,a,t) - (Irq)Fa(_,¢,a,t)) /AtAnalysis

where (rq) A corresponds to the pressure weighted Specific Humidity after the analysis,

(Trq)FG corresponds to the pressure weighted Specific Humidity first guess, and AtAn_ty_i_

is the time-step between analyses (6 hours).

i /

58) OMEGA Pressure Coordinate Vertical Velocity (mb/day)

OMEGA, which is the change in pressure in time following the parcel, is diagnostically

obtained by integrating the GCM's continuity equation. We have

dp _ (_r& + ira) OMEGA = dt

since p = 7ra + PTop. Now,

(971"

7? - + V. V_

(gt

&r (9u&

(gt - v. (Try) - 0--7

47
Integratingthe continuityequationfrom the top of the modelto agivenmodeledge-level,

ae, and using the boundary condition & Ia=o = O, we have

_r&e = - f0a_ V. (rV)da - ae-_Or

Combining the above equations, we have

OMEGA = - V • (_rV)de + aeV •VTr

Note that in this formulation, OMEGA is computed at GCM model edge-levels. OMEGA

at the model top edge is taken to be 0, thus the first computed value of OMEGA is at

sigma edge 2, or the bottom edge of the top model sigma layer.

59) DUDT Total Zonal U-Wind Tendency (m/sec/day)

DUDT is the total time-tendency of the Zonal U-Wind due to Hydrodynamic, Diabatic,

and Analysis forcing.

Ou Ou On Ou

DUDT = O--_Dynamics -_ O---_Moist -t- O--t Turbulence _- O---_Analysis

Note, in Version 1 of the GEOS GCM there are no u-wind tendencies due to cumulus fricOn -- O. tion, (ie. 0u = 0). Also, in forecast simulations -_A_alysis -- -_Moist

60) DVDT Total Zonal V-Wind Tendency (m/sec/day)

DVDT is the total time-tendency of the Meridional V-Wind due to Hydrodynamic, Diabatic, and Analysis forcing.

Ov Ov Ov Ov

DVDT Ot Dynamics + Ot Moist -t- Ot Turbulence Ot Analysis

Note, in Version 1 of the GEOS GCM there are no v-wind tendencies due to cumulus fricOv = O. tion, (ie. 0_ -- 0). Also, in forecast simulations NAnalysis -_ Moist --

61) DTDT Total Temperature Tendency (deg/day)

DTDT is the total time-tendency of Temperature due to Hydrodynamic, Diabatic, and

Analysis forcing.

OT OT OT

DTDT - + -- + --

_t Dynamics Ot MoistProcesses Ot ShortwaveRadia_ion

OT OT OT

'_ :Or LongwaveRadiation Ot Turbulence Ot Analysis

48
• ; 'i: _

Note, this diagnostic is obtained through differentiation by parts and re-arranging the model

prognostic equations. Since the thermodynamic equation is written in the mass-weighted

flux form, we have

Ot r \-_-

where 0 = Tip _. Also, in forecast simulations, 0T --0.

Ot Analysis --

62) DQDT Total Specific Humidity Tendency (g/kg/day)

DQDT is the total time-tendency of Specific Humidity due to Hydrodynamic, Diabatic,

and Analysis forcing.

DQDT Oq

Ot Dynamics

+ Oq

Ot MoistProcesses

+ Oq

Ot Turbulence

+ Oq

Ot Analysis

Note, this diagnostic is obtained through differentiation by parts and re-arranging the model

prognostic equations. Since the moisture equation is written in the mass-weighted flux form,

we have Oq { O_q 0_ ) /r o--Y= \ ot - q-_

Also, in forecast simulations, O0-_tAnalysi s -_ O.

63) VORT Relative Vorticity (sec -1)

The relative vorticity used in the formulation of the potential energy and enstrophy conserving scheme described in Section 2 is stored as the diagnostic VORT. The relative vorticity

is given as:

VORT -

acos¢

This diagnostic is stored on the Arakawa C-grid used in the GEOS GCM at the "vorticitypoint" locations, ie., centered between the u-points and v-points, including values well

defined at both poles.

64) Disabled

65) Disabled

• )

49
66) Disabled

67) USTAR Surface-Stress Velocity (m/sec)

The surface stress velocity, or the friction velocity, is the wind speed at the surface layer

top impeded by the surface drag:

k

USTAR = CuWs where : Cu = --

Cm

C_ is the non-dimensional surface drag coefficient (see diagnostic number 10), and Ws is

the surface wind speed (see diagnostic number 28).

68) Z0 Surface Roughness Length (m)

Over the land surface, the surface roughness length is interpolated to the local time from

the monthly mean data of Dorman and Sellers (1989). Over the ocean, the roughness length

is a function of the surface-stress velocity, u,.

_tt2 C5 Z0=clu 3+c2 ,+c3u,+c4+--

where the constants are chosen to interpolate between the reciprocal relation of Kondo(1975)

for weak winds, and the piecewise linear relation of Large and Pond(1981) for moderate to

large winds.

69) FRQTRB Frequency of Turbulence (0 - 1)

The fraction of time when turbulence is present is defined as the fraction of time when the

_urbulent kinetic energy exceeds some minimum value, defined here to be 0.005 ra2/sec 2.

When this criterion is met, a counter is incremented. The fraction over the averaging interval is reported.

70) PBL Planetary Boundary Layer Depth (mb)

The depth of the PBL is defined by the turbulence parameterization to be the depth at

which the turbulent kinetic energy reduces to ten percent of its surface value.

PBL = PPBL -- Psurface

where PPBL is the pressure in mb at which the turbulent kinetic energy reaches one tenth

5O
of its surfacevalue,andPs is the surface pressure.

71) SWCLR Clear sky Heating Rate due to Shortwave Radiation (deg/day)

The net Shortwave heating rate is calculated as the verticM divergence of the net solar

radiative fluxes. The clear-sky and cloudy-sky shortwave fluxes are cMculated separately.

For the clear-sky case, the shortwave fluxes and heating rates are computed with both

CLMO (maximum overlap cloud fraction) and CLRO (random overlap cloud fraction) set to

zero (see Section 3.2.3). The shortwave routine is then called a second time, for the cloudysky case, with the true time-averaged cloud fractions CLMO and CLRO being used. In all

cases, a normalized incident shortwave flux is used as input at the top of the atmosphere.

The heating rate due to Shortwave Radiation under clear skies is defined as:

OpcpT

Ot

0 NET

- ozF(clear)sw "RADSWT,

or

SWCLR - g 0 r(clear)sNwET. RADSWT.

CpTr _o"

where g is the accelation due to gravity, Cp is the heat capacity of air at constant pressure, RADSWT is the true incident shortwave radiation at the top of the atmosphere (See

Diagnostic #48), and

r(clear)N_ = F(clear)_ W - F(clear)_sw

72) OSR Net upward Shortwave flux at the top of the model (Watts/m 2)

_NET

OSR : _'SW, top

where top indicates the top of the first model layer used in the shortwave radiation routine.

In the GEOS-1 GCM, pSWtop = 0 mb.

73) OSRCLR Net upward clearsky Shortwave flux at the top of the model

(Watts/m 2)

OSRCLR NET = F(clearsky)sw, top

where top indicates the top of the first model layer used in the shortwave radiation routine.

In the GEOS-1 GCM, PSWtop = 0 mb.

51
74) CLDMAS Convective Cloud Mass Flux (kgm/sec)

The cumulative cloud mass flux per unit time, from all convective clouds generated during

the output interval, is written:

CLDMAS = _mB

where _/is the entrainment, normalized by the cloud base mass flux, and mB is the cloud

base mass flux. ms and y are defined explicitly in Section 3.1.1, the description of the

convective parameterization.

75) UAVE Time-Averaged Zonal U-Wind (m/sec)

The diagnostic UAVE is simply the time-averaged Zonal U-Wind over the NUAVE output

frequency. This is contrasted to the instantanious Zonal U-Wind which is archived on the

Prognostic Output data stream.

UAVE = u(£, ¢, a, t)

Note, UAVE is computed and stored on the staggered C-grid.

76) VAVE Time-Averaged Meridional V-Wind (m/sec)

The diagnostic VAVE is simply the time-averaged Meridional V-Wind over the NVAVE

output frequency. This is contrasted to the instantanious Meridional V-Wind which is

archived on the Prognostic Output data stream.

VAVE = v()_, ¢, a, t)

Note, VAVE is computed and stored on the staggered C-grid.

77) TAVE Time-Averaged Temperature (Kelvin)

The diagnostic TAVE is simply the time-averaged Temperature over the NTAVE output

frequency. This is contrasted to the instantanious Temperature which is archived on the

Prognostic Output data stream.

TAVE = T()_, ¢, a, t)

78) QAVE Time-Averaged Specific Humidity (g/kg)

52
ThediagnosticQAVE issimplythetime-averagedSpecificHumidityovertheNQAVE output frequency.This is contrastedto the instantaniousSpecificHumiditywhichis archived

onthe PrognosticOutput datastream.

QAVE = q(A, ¢, a, t)

79) Disabled

80) PAVE Time-Averaged Surface Pressure - PTOP (mb)

The diagnostic PAVE is simply the time-averaged Surface Pressure - PTOP over the

NPAVE output frequency. This is contrasted to the instantanious Surface Pressure -

PTOP which is archived on the Prognostic Output data stream.

PAVE = _r(A,¢,a,t)

= ps(A,¢,a,t)--pT

81) QQAVE Time-Averaged Turbulent Kinetic Energy (m/sec) 2

The diagnostic QQAVE is simply the time-averaged prognostic Turbulent Kinetic Energy

produced by the GEOS-1 GCM Turbulence parameterization over the NQQAVE output frequency, This is contrasted to the instantanious Turbulent Kinetic Energy which is

archived on the Prognostic Output data stream.

QQAVE = qq(A, ¢, a, t)

Note, QQAVE is computed and stored at the "mass-point" locations on the staggered

C-grid.

82) SWGCLR Net downward clearsky Shortwave flux at the surface (Watts/m 2)

b: /

SWGCLR = F(clearsky)N_,ntay+l

= F(clearsky)lsW, nlay+l --F(clearsky)tsW, nlay+l

where nlay+l indicates the lowest model edge-level, or p = ps_r]. F(clearsky)SW _ is the

53
downwardclearskyShortwaveflux and F(clearsky)_ W is the upward clearsky Shortwave

flUX.

83) ANALP Analysis Surface Pressure Tendency (mb/sec)

ANALP is the surface pressure tendency created during the assimilation cycle which is

used in the Incremental Analysis Updating (IAU) procedure.

ANALP = (TrA(_, ¢, a, t) - _rF_(_,¢,_,t)) /AtAnalysis

where 7r A corresponds to the surface pressure after the analysis, 7rFG corresponds to the

surface pressure first guess, and AtAnalysis is the time-step between analyses (6 hours).

84) SDIAG1 User-Defined Surface Diagnostic-1

The GEOS GCM provides Users with a built-in mechanism for archiving user-defined diagnostics. The generic diagnostic array QDTAG located in COMMON /DIAG/, and the

associated diagnostic counters and pointers located in COMMON /DIAGP/, must be accessable in order to use the user-defined diagnostics (see Section 6). A convenient method

for incorporating all necessary COMMON files is to include the GEOS GCM vstate.com file

in the routine which employs the user-defined diagnostics.

In addition to enabling the user-defined diagnostic (ie., CALL SETDIAG(84)), the User

must fill the QDIAG array with the desired quantity within the User's application program

or within modified GEOS GCM subroutines, as well as increment the diagnostic counter at

the time when the diagnostic is updated. The QDIAG location index for SDIAG1 and its

corresponding counter is automaticMly defined as ISDIAG1 and NSDIAG1, respectively,

after the diagnostic has been enabled. The syntax for its use is given by

do j=i,jm

do i=i,im

qdiag(i,j,ISDIAGi) = qdiag(i,j,ISDIAGi) +

enddo

enddo

NSDIAGI = NSDIAGi + i

The diagnostics defined in this manner will automatically be archived by the production

output programs (see Section 11).

85) SDIAG2 User-Defined Surface Diagnostic-2

54
/

The GEOS GCM provides Users with a built-in mechanism for archiving user-defined diagnostics. For a complete description refer to Diagnostic #84. The syntax for using the

surface SDIAG2 diagnostic is given by

do j=l,jm

do i=l,im

qdiag(i,j,ISDIAG2) = qdiag(i,j,ISDIAG2) +

enddo

enddo

NSDIAG2 = NSDIAG2 + 1

The diagnostics defined in this manner will automatically be archived by the production

output programs (see Section 11).

86) UDIAG1 User-Defined Upper-Air Diagnostic-1

The GEOS GCM provides Users with a built-in mechanism for archiving user-defined diagnostics. For a complete description refer to Diagnostic #84. The syntax for using the

upper-air UDIAG1 diagnostic is given by

do L=l,nlay

do j=l,jm

do i=l,im

qdiag(i,j,IUDIAGl+L-l) = qdiag(i,j,IUDIAGl+L-l) + ...

enddo

enddo

enddo

NUDIAGI = NUDIAGI + i

The diagnostics defined in this manner will automatically be archived by the production

output programs (see Section 11).

iii

/P i / r

87) UDIAG2 User-Defined Upper-Air Diagnostic-2

The GEOS GCM provides Users with a built-in mechanism for archiving user-defined diagnostics. For a complete description refer to Diagnostic #84. The syntax for using the

upper-air UDIAG2 diagnostic is given by

55
do L=l,nlay

do j=i,jm

do i=i,im

qdiag(i,j,IUDIAG2+L-1) = qdiag(i,j,IUDIAG2+L-1) + ...

enddo

enddo

enddo

NUDIAG2 = NUDIAG2 + i

The diagnostics defined in this manner will automatically be archived by the production

output programs (see Section 11).

88) DIABU Total Diabatic Zonal U-Wind Tendency (m/sec/day)

DIABU is the total time-tendency of the Zonal U-Wind due to Diabatic processes and the

Analysis forcing.

Ou Ou Ou

+-- +-- DIABU - Ot Moist Ot Turbulence Ot Analysis

Note, in Version 1 of the GEOS GCM there are no u-wind tendencies due to cumulus fricOu --O. tion, (ie. 0u - 0). Also, in forecast simulations -O-tAnalysis -- -OTMoist --

89) DIABV Total Diabatic Meridional V-Wind Tendency (m/sec/day)

DIABV is the total time-tendency of the Meridional V-Wind due to Diabatic processes

and the Analysis forcing.

Ov Ov Ov

DIABV = O--tMoist -_ O--tTurbulence _- O--tAnalysis

Note, in Version 1 of the GEOS GCM there are no v-wind tendencies due to cumulus fricOv --O. tion, (ie. 0v = 0). Also, in forecast simulations NAnalysis -- -5-tMoist

90) DIABT Total Diabatic Temperature Tendency (deg/day)

DIABT is the total time-tendency of Temperature due to ,Diabatic processes and the

Analysis forcing.

OT OT

DIABT - + --

Ot MoistProcesses Ot ShortwaveRadiation

OT OT OT

__ +-- +--

-_ Ot LongwaveRadiation Ot Turbulence Ot Analysis

56
::::: :

!i:il'_

If we define the time-tendency of Temperature due to Diabatic processes as

OT OT OT

O---t Diabatic -- Ot MoistProcesses "_- O---t ShortwaveRadiation

OT OT

+ +-

Or LongwaveRadiation Ot Turbulence

then, since there are no surface pressure changes due to Diabatic processes, we may write

OT _ p'_ O_rO

Ot Diabatic _ Ot Diabatic

where 8 = T/p% Thus, DIABT may be written as

DIABT = -_- \ -_ Diabatic + Ot Analysis]

Note, this diagnostic is only approximate since we have neglected the affect of the analysis

on the surface pressure, ie.

OT p_ ( 0_0

Ot Analysis- \ Analysis+On

Ot Analysis L p

This term will be added in subsequent versions of the GEOS GCM. In forecast simulations

there is no error since OT = O.

-OT Analysis

i •

91) DIABQ Total Diabatic Specific Humidity Tendency (g/kg/day)

DIABQ is the total time-tendency of Specific Humidity due to Diabatic processes and the

Analysis forcing.

Oq + Oq Oq

DIABQ = -_ MoistProcesses O--_Turbuleuce -_ O---tAnalysis

If we define the time-tendency of Specific Humidity due to Diabatic processes as

0q = 0A + 0q

Ot Diabatic Ot MoistProcesses Ot Turbulence

then, since there are no surface pressure changes due to Diabatic processes, we may write

Oq 10rq

Ot Diabatic 7¢ Ot Diabatic

Thus, DIABQ may be written as

!

DIABQ = % }

7r \ Ot Diabatic Ot Analysis]

57
Note,this diagnosticis only approximatesincewehaveneglectedthe affectof the analysis

on the surfacepressure,ie.

/

Oq _ 1 ( O_q. -- q-_ Analysis

Ot Analysis 7_ \ Ot Analysis

This term will be added in subsequent versions of the GEOS GCM. In forecast simulations

there is no error since O_-_tAnalysi s -_ O.

92) Disabled

93) Disabled

94) Disabled

95) Disabled

96) Disabled

97) Disabled

98) Disabled

99) Disabled

100) Disabled

101) Disabled

102) VINTUQ Vertically Integrated Moisture Flux (m/sec. g/kg)

The vertically integrated moisture flux due to the zonal u-wind is obtained by integrating

uq over the depth of the atmosphere at each model timestep, and dividing by the total mass

58
of thecolumn.

rt°P uqpdz

VINTUQ : Jsurf

rt°P pdz Jsurf

Using p_fz = -_P = -_-_fa, we have

g g

VINTUQ = 1 uqda

103) VINTVQ Vertically Integrated Moisture Flux (m/sec. g/kg)

The vertically integrated moisture flux due to the meridional v-wind is obtained by integrating vq over the depth of the atmosphere at each model timestep, and dividing by the

total mass of the column.

rt°P vqpdz

VINTVQ- Jsurf

I't°P pdz Jsur]

Using p_z = -_ = -_-Sa, we have

g g

VINTVQ = 1vqda

104) VINTUT Vertically Integrated Heat Flux (m/sec. deg)

The vertically integrated heat flux due to the zonal u-wind is obtained by integrating uT

over the depth of the atmosphere at each model timestep, and dividing by the total mass

of the column.

rt°P uT pdz

VINTUT- J_r/

f top pdz urf

Using pSz = -@g = -_Sa, we have

1

VINTUT = uTda

105) VINTVT Vertically Integrated Heat Flux (m/sec. deg)

59
Theverticallyintegratedheatflux dueto the meridionalv-windis obtainedby integrating

vT over the depth of the atmosphere at each model timestep, and dividing by the total

mass of the column.

top Z

VINTVT = fiery vTpd

fst°p pdz urf

Using pSz = -_ = -_-Sa, we have

g g

jfO 1

VINTVT = vTda

106 CLDFRC Total 2-Dimensional Cloud Fracton (0- 1)

The 2-dimensional net cloud fraction as seen from the top of the atmosphere is given by

NLAY

CLDFRC =1- (1- MAXN_ AY [CLMO,]) 1-I

l=ll

(1 - CLROI)

For a complete description of cloud/radiative interactions, see Section 3.2.3.

107) QINT Total Precipitable Water (gm/cm 2)

The Total Precipitable Water is defined as the vertical integral of the specific humidity,

given by:

top pqdz QINT = J_r]

_01 = - qda

g

where we have used the hydrostatic relation pSz = - @ = -_-6a.

g g

108) U2M Zonal U-Wind at 2 Meter Depth (m/sec)

The u-wind at the 2-meter depth is determined from the similarity theory:

U, ttsl Cm2m

U2M = --'k Cm2m _V s -- _mst

ttsl

where Cm(2m) is the non-dimensional wind shear at two meters, and the subscript sl refers

to the height of the top of the surface layer. If the roughness height is above two meters,

6O
U2M is undefined.

109) V2M Meridional V-Wind at 2 Meter Depth (m/sec)

The v-wind at the 2-meter depth is a determined from the similarity theory:

V2M = u,/, vsl = ¢m2m VsZ

k _'_2.n W_ _p,_

where Cm(2m) is the non-dimensional wind shear at two meters, and the subscript sl refers

to the height of the top of the surface layer. If the roughness height is above two meters,

V2M is undefined.

110) T2M Temperature at 2 Meter Depth (deg K)

The temperature at the 2-meter depth is a determined from the similarity theory:

T2M = P'_(Ok(¢h2m + Cg) _- 0surf) = P'_(Os_ I +

+ ¢.

where:

where Ch(2m) is the non-dimensional temperature gradient at two meters, ¢g is the nondimensional temperature gradient in the viscous sublayer, and the subscript 8l refers to the

height of the top of the surface layer. If the roughness height is above two meters, T2M is

undefined.

....(i:i _

, i _ • /

111) Q2M Specific Humidity at 2 Meter Depth (g/kg)

The specific humidity at the 2-meter depth is :determined from the similarity theory:

k p,_ Ch2m + Cg (q_t -- qsur/)) Q2M = P_( (¢h2m + Cg) + qsurl) = (qsurf + Ch_z + Cg

where:

where Ch(2m) is the non-dimensional temperature gradient at two meters, Cg is the nondimensional temperature gradient in the viscous subtayer, and the subscript sl refers to the

61
heightof the top of the surfacelayer.If the roughnessheightis abovetwo meters,Q2M is

undefined.

112) U10M Zonal U-Wind at 10 Meter Depth (m/sec)

The u-windat the 10-meterdepthis aninterpolationbetweenthe surfacewind and the

modellowestlevelwind usingtheratioof the non-dimensionalwind shearat thetwolevels:

U, USl CmlOm

U10M = T_2mlom Ws -- _2ms I Usl

where Cm(10m) is the non-dimensional wind shear at ten meters, and the subscript sl refers

to the height of the top of the surface layer.

113) V10M Meridional V-Wind at 10 Meter Depth (m/sec)

The v-wind at the 10-meter depth is an interpolation between the surface wind and the

model lowest level wind using the ratio of the non-dimensional wind shear at the two levels:

_* Vsl CmlOrn

VIOM = TCmlOm Vsl

where Cm(10m) is the non-dimensional wind shear at ten meters, and the subscript sl refers

to the height of the top of the surface layer.

114) T10M Temperature at 10 Meter Depth (deg K)

The temperature at the 10-meter depth is an interpolation between the surface potential

temperature and the model lowest level potential temperature using the ratio of the nondimensional temperature gradient at the two levels:

T10M = P (_-(¢hl0m + Cg) + es_r/) = P'(0s_rS + CJ_lom+ ¢'(e_z - e_uri))

where:

0. - - (_'0')

where _h(10m) is the non-dimensional temperature gradient at two meters, Cg is the nondimensional temperature gradient in the viscous sublayer, and the subscript 81 refers to the

height of the top of the surface layer.

62
115) QIOM Specific Humidity at 10 Meter Depth (g/kg)

The specific humidity at the 10-meter depth is an interpolation between the surface specific

humidity and the model lowest level specific humidity using the ratio of the non-dimensional

temperature gradient at the two levels:

Q10M = _ q_ _ __ _ (qsl - qsur/)) P (_-(¢hl0m + Cg) + qsur/) = P_(qs=rf + Chl0m + Cg,

where:

q, m

where Ch(10m) is the non-dimensional temperature gradient at two meters, Cg is the nondimensional temperature gradient in the viscous sublayer, and the subscript 81 refers to the

height of the top of the surface layer.

116) DTRAIN Convective Cloud Detrainment (kg m/sec)

The cloud mass flux at the detrainment level per unit time, from all convective clouds

generated during the output interval, is written:

DTRAIN = _aDmB

where aD is the detrainment a level, ms is the cloud base mass flux, and _ is the entrainment, defined in Section 3.1.1.

117) QFILL Filling of negative Specific Humidity (g/kg/day)

Due to computational errors associated with the numerical scheme used for the advection

of moisture, negative values of specific humidity may be generated. The specific humidity

is checked for negative values after every dynamics timestep. If negative values have been

produced, a filling algorithm is invoked which redistributes moisture from below (see Section

2.3). Diagnostic QFILL is equal to the net filling needed to eliminate negative specific

humidity, scaled to a per-day rate:

QFILL = qn+l _ qn+l final initial

where

qn+l = (a.q)n+l/71.n+l

118) VINTQANA Precipitable Moisture Adjustment from Analysis (mm/day)

63
/i

For the time change in specific humidity due to the analysis increment,

vertical integral or total precipitable amount is given by:

0q the

Ot analysis '

flop Oq dz = - top Oq -- dp - 7r fo i Oq VINTQANA = Js_Tf P-_analysis Js_Tf Ot analysis g g -_analysis d(7

Defining the specific humidity analysis increment, _analysis' aS

Oq 1 (07rq 07r

O--tanalysis -- _r \ Ot analysis- q-_ analysis ]

we may write

VINTQANA= lfol(O_rq O_r ) g k, Ot analysis q-_ analysis do"

119) VINTQFIL Precipitable Moisture Adjustment from Filling (mm/day)

When advective or other processes create spurious negative specific humidity, total moisture

is conserved by filling the negative amounts by "borrowing" moisture from below. When

there is not enough moisture in a column to adequately fill the negative values, the specific

humidities are set to zero, creating an artificial source of moisture.

For a change in specific humidity due to filling, Aq]itl, the vertical integral or total precipitable amount is given by:

VINTQFIL = fop

Jsur]

pAqfilldz = lop --

J sur ]

Aq]ill dp _ _ fo 1 g g Aq]iuda

A precipitable adjustment rate is defined as the vertically integrated moisture adjustment

per Dynamics time step, scaled to mm/day.

6•4
Part III

GEOS-1 GCM: User's Guide

9 Introduction

In this section we discuss the general philosophy used in the construct of the GEOS GCM

system, the necessary components of the system, and a basic User's Guide on how to set-up

and run experiments.

10 GEOS Superstructure

The superstructure of the GEOS GCM is designed so that the entire model may be invoked

by a simple subroutine call. With each call of the main GCM driver, all model prognostic

and diagnostic fields are updated by one timestep. A variety of applications may be built

around the main GCM driver (GCMDRV), which allow for straight forecast simulations,

assimilation cyles, and dynamic initialization procedures. It should be emphasized that

what is meant by the GEOS System is simply the library of GCM subroutines which include

GCMDRV and those that are subsequently called. These subroutines will be denoted as

below GCMDRV. Those routines which are above GCMDRV include the Main Program, or

the Application, and any model Output routines.

10.1 GEOS Applications

The GEOS Application can be defined as the Main Program from which the GEOS GCM

and/or the GEOS Utility subroutines are called. All Applications are User-supplied, that

is, they do not belong to the GEOS GCM compiled library. The purpose of the GEOS

Application is a pre-determined experiment, or application, which the User wants to perform

using the GEOS GCM library.

For example, let us create an Application which performs the following funtions:

J

:/

Read a model restart from Unit1

Integrate the model forward in time for NMAX iterations

Write the resulting model restart to Unit2

65
This couldbeaccomplishedwith thefollowingFORTRANcode:

CALL RESTART(UNITI,IREAD,...)

DO N = 1,NMAX

CALL GCMDRV(Argl, Arg2, Arg3, ...)

ENDDO

CALL RESTART (UNIT2,IWRITE,...)

Here IREAD and IWRITE are read/write IO flags used by the GEOS Utility Subroutine

RESTART to read/write GEOS Model restarts, and Argl, Arg2, Arg3,... are the the

arguments to GCMDRV. Both of these argument lists are described fully in the next

section.

We can easily modify this application to run forecasts from the same initial state but using

different GCMDRV model parameters:

CALL RESTART (UNITI,IREAD,...)

DO N = 1,NMAX

CALL GCMDRV (Argl.1, Argl.2, Argl.3, ...)

ENDDO

REWIND UNIT1

CALL RESTART (UNITI,IREAD,...)

DO N = 1,NMAX

CALL GCMDRV (Arg2.1, Arg2.2, Arg2.3, ...)

ENDDO

CALL RESTART (UNIT2,IWRITE,...)

Finally, we could run the model with the same set of GCMDRV model parameters but

starting from two different initial states:

CALL RESTART (UNITI,IREAD,...)

66
DO N = 1,NMAX

CALL GCMDRV(Argl, Arg2, Arg3, ...)

ENDDO

CALL RESTART (UNIT2,IWRITE,...)

CALL RESTART (UNIT3,IREAD,...)

DO N = 1,NMAX

CALL GCMDRV (Argl, Arg2, Arg3, ...)

ENDDO

CALL RESTART (UNIT4,IWRITE,...)

These simple examples show the ease and flexibility of creating applications to perform

experiments with the GEOS GCM. In the next section we describe fully the argument

lists for subroutines GCMDRV and RESTART, and briefly describe GCMDRV's flow

structure.

10.2 Subroutines GCMDRV g¢ RESTART

Subroutine GCMDRV is the single interface between User Application routines and the

actual updating of model prognostic fields. Its flow structure can be separated into three

essential parts:

.!

Read and refresh Boundary Conditions

Compute time-tendencies of atmospheric state variables due to physical and hydrodynamicM processes

Integrate the atmospheric state variables forward in time by one

timestep

Every call to GCMDRV updates the boundary conditions to the current date and time. In

Version 1 of the GEOS GCM this is accomplished by a linear interpolation between monthly

mean data sets. Subsequent versions of the GCM will allow for boundary condition data

sets with varying fl'equency intervals. The boundary conditions which are currently updated

are:

• Sea Surface Temperature

67
• GroundWetness

• SurfaceAlbedo

• SurfaceRoughness

• Sea-IceThickness

In addition,datais provideddescribingthe modeltopography,topographyvariance,and

the Land/Water/Glacierflags. Thesedata arenot interpolatedto the currentdate and

time sincethey arefixedquantities.Finally, the incidentsolarradiationat the top of the

atmosphereis computedfor the currentdateandtime.

At prescribedtime intervals,GCMDRV calls separatemini-driversto updatethe timetendenciesof the prognosticstatevariablesdueto particularphysicalor hydrodynamical

processes.Eachmini-driverhasits owntimescMeof operation.Themini-driversassociated

with physicalandhydrodynamicalprocessesare:

Moist Processes(Updatedevery10minutes)

ShortwaveRadiation(Updatedevery3 hoursusinga normalizedincidentradiation,andadjustedeverytimestepusingthetrue incident

radiation)

LongwaveRadiation(Updatedevery3 hours)

Turbulence(Updatedevery30minutes)

EulerianHydrodynamics(Updatedeverytimestep)

The time4endenciesfortheprocessesstatedabovearealwayskeptin coreduringthe GEOS

GCM simulation.Thus,all prognosticfieldsaretime-continuouswhilethe time-tendencies

of the prognosticfieldsdue to the variousprocessesarediscontinousat the abovestated

frequencies.Thesetime-tendenciesaresummedtogetherto providea total time-tendency

to be usedin the time-integrationscheme.Currently,the GEOSGCM canberun with

eitherthe Matsunoor the Leapfrogtime scheme.The Leapfrogschemehasassociatedwith

it an Asselin(1972)time filter.

68
The callingsequenceassociatedwith GCMDRV and a description of the argument list is

provided below:

CALL GCMDRV ( KBC,N,NYMD,NHMS,SCHEME,NDT,ALPHA,QQINIT,QBUD )

PARAMETER TYPE DESCRIPTION

:i_i_ _i _i

KBC Integer

N Integer

NYMD Integer

NHMS Integer

SCHEME Character*4

NDT Integer

ALPHA Real

QQINIT Logical

QBUD Logical

KBC is the unit from which Boundary Conditions are

read. Every call to GCMDRV will initiate an update of

the boundary conditions. Boundary conditions are updated by a linear interpolation between the two months

kept in memory which strap the current date and time.

Each time a new mid-month is passed or a new KBC value

is encountered, two new months are read in. Setting KBC

equal to zero disables all boundary condition updates.

N refers to the current time index, and is equal to 1 or

2. Its value always points to the updated fields and is

modified by GCMDRV.

NYMD is the current year/month/day in YYMMDD format. Its value is modified by GCMDRV.

NHMS is the current hour/minute/second in HHMMSS

format. Its value is modified by GCMDRV.

SCHEME is the variable indicating the time scheme to be

used for the integration. SCHEME may be set to 'MATS'

for the Matsuno time scheme, or 'LEAP' for the Leapfrog

time scheme. The SCHEME variable may be changed at

any time during the forecast.

NDT is the timestep (in seconds) used for the forecast.

NDT may be positive, zero, or negative, and may be

changed at any time during the forecast.

Asselin Filter Coefficient, a, to be used with the Leapfrog

Time Scheme, where _Filt '/n = qn(1 - a) + a ( qn'[" l "_2qn--1 ).

ALPHA may be changed at any time during the forecast.

QQINIT is a Logical flag parameter used by the Turbulence parameterization scheme. QQINIT should be set to

TRUE if the GEOS GCM is starting from an initial condition which does not contain Turbulent Kinetic Energy

as a prognostic state variable. After one successful call to

GCMDRV, QQINIT should be set to FALSE.

Logical flag to compute budget diagnostics.

69
The GEOSGCMUtility subroutineRESTART maybeusedwithin User-DefinedApplications to readand/or write a GEOSGCMRestartdataset.Thecallingsequenceassociated

with RESTART anda descriptionof the argumentlist is providedbelow:

CALL RESTART( KU,IOFLAG,NYMD,NHMS,N,NYMD0,NHMS0,JOB,NTYPE)

PARAMETER TYPE DESCRIPTION

KU Integer

IOFLAG Integer

NYMD Integer

NHMS Integer

N Integer

NYMD0 Integer

NHMS0 Integer

JOB Character*8

NTYPE Integer

Unit numberfrom whichto reador to write the GEOS

Restart

IO Flagwhichis equalto: 1 to Reada GEOSRestart,2

to Write a GEOSRestart.

NYMD is the currentyear/month/dayin YYMMDD format.

NHMSis the currenthour/minute/secondin HHMMSS

format. N refersto the currenttime index, andis equalto 1 or

2. Its valuealwayspoints to the mostrecentlyupdated

fields.

NYMD0 is the beginningyear/month/dayof the current

experimentin YYMMDD format. NYMD0isusedin current ProductionApplicationsto controlOutput frequen-.

cies.

NHMS0is the beginninghour/minute/secondof the current experimentin HHMMSSformat. NYMD0 is used

in currentProductionApplicationsto controlOutput frequencies.

An ExperimentIdentifier(usuallya nameor experiment

number).

NTYPEis theGEOSRestartIdentificationindex.

NTYPE = -1 refersto a "First Guess"restart whichwas

usedfor an analysis. NTYPE = 0 refers to the "After

Analysis" or straight forecast restart. NTYPE = 1 refers

to a restart from an IAU assimilation.

It should be noted that if IOFLAG = 1, the RESTART parameters NYMD, NHMS, N,

NYMD0, NHMS0, JOB, and NTYPE are read from the GEOS Restart dataset and passed

back to the calling program. If IOFLAG = 2, the parameters NYMD, NHMS, N, NYMD0,

NHMS0, JOB, and NTYPE passed into subroutine RESTART will be written out to the

GEOS Restart.

7O
11 GEOS Production

Although the GEOS GCM system is designed to work with a User supplied "Applications"

program, several application programs have been developed by the DAO which are used

in standard production forecast simulations and analyses, and are available as templates.

These applications provide examples of convenient ways to control the input and output

streams, the frequency and selection of diagnostic and prognostic fields for output, and

the simulation start and end times. The application routine examples reside on the SGI

workstation hera.gsfc.nasa.gov, in the directory/production/geos_das/geosl.1/applications.

Some of them are described here in detMl:

gcmprod.f A simple application routine designed to read and follow namelist-supplied

directions to simulate a specific period and direct all "activated" diagnostics

to an output stream. The application reads one namelist to establish the

initial and final time and date of the simulation, the frequency of the output stream, and other run control parameters such as timestep and asselin

filter strength. Another namelist is read in which logical flags determine the

selection of diagnostics to be activated and directed to the output stream.

gcmprodn.f An application routine which is designed to handle a somewhat more complex

output stream for diagnostic and prognostic fields. The namelist controlled

interface for the diagnostics involves the specification of an HHMMSS format variable for each activated diagnostic separately. Diagnostics may be

averaged for any time period up to the length of the simulation, and the averaging period may be centered about the current time or may extend from

the previous time to the present. The application determines how many different averaging periods and types are present, and directs the output to the

appropriate number of separate output streams. Separate output routines

are called to properly direct the diagnostic output. Control over the begin

and end times of the simulation are also governed by namelist parameters.

amip.f An experiment-specific application routine that was designed to handle the

control of an extended multi-year simulation with complex diagnostic output

requirements. The application was written to direct diagnostic output to

several output data streams each with a pre-specified set of diagnostics.

Separate output routines are called for each data stream, providing them

with the required list of diagnostic quantities.

// i

? /

In addition to flexible application routines, the GEOS GCM system also incorporates user

defined output routines that can be called by the application. The output routines developed

at the DAO write a combination of prognostic and/or diagnostic quantities, on model sigma

levels or interpolated to pressure levels, in the Phoenix Output Format (see Section 11.2).

71
Justlike theapplicationroutines,severaloutputroutineshavebeendevelopedby the DAO

andareavailableasexamples.Someof themaredescribedhere:

diagprs.f An output routine whichacceptsasargumentsa specificlist of diagnostic

quantitiesandinterpolatesthemto anargumentspecifiedlist ofpressurelevels. Theroutineusesthe modelprovidedGETDIAG subprogramto obtain

the propertime-averagedquantities.The codeis written to distinguishbetweendifferentquantitiesandattemptsto useanappropriateinterpolation

algorithm.Output is written in PhoenixFormat.

diagsig.f Designedto acceptanargument-specifiedlist of diagnostics, and direct output in Phoenix Format to one output stream. The GETDIAG routine is

used to retrieve diagnostic quantities in the proper units from the model.

glaoutc.f When this output routine is cMled, M1 prognostic quantities and all activated

diagnostic quantities are interpolated to an argument-list specified set of

pressure levels. All output is directed to a Phoenix Format data stream.

sigoutc.f - All model prognostic and activated diagnostic quantities are directed to a

single sigma-level Phoenix Output Format data set.

progprs.f Designed to direct model prognostic quantities to a single output stream.

Quantities are interpolated to a list of pressure levels.

progsig.f Model prognostic quantities are simply directed to a single sigma-level

Phoenix format data set.

72
11.1 Production Namelists

Two general namelists have been designed to setup and control the Production runs of the

GEOS GCM. The namelists include information about the beginning date and time of the

forecast, the ending date and time of the forecast, the duration of each job segment, and

output characteristics.

The first namelist is called INPUT, and contains the parameters to control the beginning

and ending times of each experiment, as well as some output information. The following

INPUT namelist and subsequent description can be used as is or as a template for Users to

build upon within their own Applications:

_INPUT

JOB

XLAB

NYMDB

NHMSB

NYMDE

NHMSE

NDT

LEAP

MATS

QQINIT

ALPHA

NDSEG

NDOUT

PLEVS

_END

= ' 1234',

= 'E1234 GEOS-GCM

'RIMENT

= 850101

= 000000

= 860101

= 000000

= 225

= .TRUE.

= .FALSE.

= .FALSE.

= 0.05

= 2400000

= 060000

= 1000.0,

700.0,

200.0,

CLIMATE SIMULATION EXPE',

950.0 , 900.0, 850.0, 800.0,

600.0 , 500.0, 400.0, 300.0,

150.0 , I00.0, 70.0, 50.0,

73
PARAMETER TYPE DESCRIPTION

JOB Character*8

XLAB Character*80

NYMDB Integer

NHMSB Integer

NYMDE Integer

NHMSE Integer

NDT Integer

LEAP Logical

MATS Logical

QQINIT Logical

ALPHA ReM

NDSEG Integer

NDOUT Integer

PLEVS Real

Experiment Identifier (Usually a name or experiment

number)

User Assigned Description of Current Experiment

Beginning Date of Experiment in YYMMDD format

Beginning Time of Experiment in HHMMSS format

Ending Date of Experiment in YYMMDD format

Ending Time of Experiment in HHMMSS format

Time step (in seconds) of the hydrodynamics

Logical Flag to use the Leapfrog Time Scheme

Logical Flag to use the Matsuno Time Scheme

QQINIT is a Logical flag parameter used by the Turbulence parameterization scheme. QQINIT should only be

set to TRUE if the GEOS GCM is starting from an initial

condition which does not contMn Turbulent Kinetic Energy as a prognostic state variable. After one successful

call to GCMDRV, QQINIT should be set to FALSE.

Asselin Filter Coefficient, a, to be used with the Leapfrog

Time Scheme, where qFilt = qn(1 _ o_) + _ ( qn+l+qn--1).

Current job segment time interval in HHMMSS format. The total time of a GEOS simulation is given as

(NYMDE,NHMSE) - (NYMDB,NHMSB). This total time

is divided into a number of smaller job segments, which

run in succession on the computer. The job segment size

is determined by the User, depending on the desired size

of the Output stream and/or the CPU time restrictions

within the computer QUEUE. Note, HH may be any number of hours, and not limited to a 24 hour day (ie., 5 day

job segments would be defined as NDSEG = 1200000).

Output frequency for Prognostic fields in HHMMSS format.

Levels used for Sigma to Pressure Interpolation to create

Pressure Level Output. (Maximum number of levels =

50)

When the GEOS GCM is run during for a particular experiment, no diagnostic quantities

will be produced or saved unless the User has explicitly enabled the desired diagnostics from

the Applications program. Using the DAO supplied applications directly or as templates,

a convenient NAMELIST method has been adopted for diagnostic control. Diagnostic

quantities which are of interest to the User may be enabled by the application by including

74
i _ /

the Diagnostic Name from the Diagnostic Menu (Section 7) in the namelist INDIAG. For

those Users employing the gcmprod.fapplication described in Section 11, simply preface the

Diagnostic Name with the letter L. This variable will then be treated as a logical parameter.

Those parameters set to TRUE will be enabled for the GEOS GCM simulation. An example

of the INDIAG namelist using the Logical variable syntax is shown below:

_INDIAG

LUFLUX = .TRUE.

LVFLUX = .TRUE.

LHFLUX = .TRUE.

LMOISTT = .TRUE.

LMOISTQ = .TRUE.

LRADLW = .TRUE.

LRADSW = .TRUE.

LPREACC = .TRUE.

LPRECON = .TRUE.,

LOLR = .FALSE.

LEVAP = .FALSE.

LU2M = .FALSE.

LVRM = .FALSE.

LT2M = .FALSE.

LQ2M = .FALSE.

_END

In this example, the diagnostic quantities UFLUX, VFLUX, HFLUX, MOISTT, MOISTQ,

RADLW, RADSW, PREACC, and PRECON will all be enabled during the simulation. The remMning diagnostics OLR, EVAP, U2M, V2M, T2M, and Q2m will be

disabled during the run. It should be noted that explicitly defining the logical diagnostic

parameters to FALSE in the INDIAG namelist is not necessary since this is the default.

Using this method in conjunction with the application gcmprod.f will produce diagnostic

quantities whose output frequency is governed by the INPUT namelist parameter NDOUT.

For those Users interested in a more sophisticated diagnostic output frequency algorithm,

each diagnostic may be written out at its own unique frequency by using the DAO application gcmprodn.f. In this case, a particular diagnostic is turned on, or enabled, by specifying

a non-zero diagnostic frequency, in HHMMSS format, for that diagnostic. The diagnostic

frequency name is the same as the diagnostic name (from the Diagnostic Menu, Section

7) preceded by the letter N (eg, NUFLUX = 030000). In this example, the diagnostic

UFLUX will be accumulated and written out every 3 hours. Every diagnostics has its own

frequency, thus certain diagnostics may be written out more frequently than others. A

separate file is defined for each unique diagnostic frequency. Thus, if the user turns on a

total of 24 diagnostics, 12 with a 6 hour frequency, 10 with a 3 hour frequency, and 2 with a

half-hour frequency, a total of 3 diagnostic datasets will be created. Finally, diagnostics are

accumulated and written out with the time-stamp at the end of the accumulation period.

75
A time-stamp which occurs at the middle of the accumlating period may be obtained by

making the diagnostic frequency negative. An example of the INDIAG namelist using the

frequency variable syntax is shown below:

_INDIAG

NUFLUX = 030000

NVFLUX = 030000

NHFLUX = 030000

NMOISTT =-060000

NMOISTQ =-060000

NRADLW =-060000

NRADSW =-060000

NPREACC = 060000

NPRECON = 060000

NOLR = 060000

NEVAP = 060000

NU2M = 010000

NV2M = 010000

NT2M = 010000

NQ2M = 010000

_END

In this example, four output data streams will be produced. Diagnostics U2M, V2M,

T2M, and Q2M will be averaged over a one-hour frequency and written out with a timestamp at the end of each hour. Diagnostics UFLUX, VFLUX, and HFLUX will be

averaged over a three-hour time-period. Diagnostics PREACC, PRECON, OLR, and

EVAP will be averaged over a six-hour period and written out with a time-stamp at the

end of the averaging period. Diagnostics MOISTT, MOISTQ, RADLW, and RADSW

will also be averaged over a six-hour period but will be written out with a time-stamp

associated with the middle of the averaging period. All other diagnostics are disabled.

76
11.2 Phoenix Format

The Phoenix Format used in the archiving of the GEOS AMIP and the GEOS-DAS 5-Year

Re-analysis employs a self-documenting output stream which describes each field and its

attributes for each time period on the output dataset. At each output time period, Header

information is given which describes the horizontal and vertical dimensions of the output,

the time period of the output, and variable descriptions of the data type. The number and

type of fields are User Defined.

The Phoenix format allows for writing data in three different ways (or configurations). The

first, configuration #1, is simply a sequence of distinct 2-dimensional fields (such as sealevel pressure, topography, or selected individual levels of specific quantities such as wind or

moisture). The second method, configuration #2, is 3-dimensional data which are ordered

by quantity for each level. It is assumed that all quantities written in this manner are

defined for all levels. The third, configuration #3, is 3-dimensional data which are ordered

by level for each quantity. Here, the number of levels associated which each quantity may

vary. The specific configurations used within any dataset are determined by the User. The

following sample code may be used to read any Phoenix Format Output stream:

PARAMETER ( IDIM = 144 )

PARAMETER ( JDIM = 091 )

DIMENSION Q( IDIM,JDIM )

! IDIM=I44 for 2.8 degree Longitude

! JDIM=091 for 2.0 degree Latitude

CHARACTER*8

DIMENSION

DIMENSION

CHARACTER*8

CHARACTER*40

XLABEL(IO) , JOB, NAME

SIGE (50) , NDLEV(IO0), ZLEV(50)

IDUM(IO0), RDUM(IO0)

NAMES(IO0), NAMEU(IO0), NAMED(IO0), CDUM(IO0)

DESCS(IO0), DESCU(IO0), DESCD(IO0), DDUM(IO0)

DATA KU /10/

READ(KU,END=800) JOB, NYMD, NHMS, NYMDO, NHMSO,

XLABEL, IM, JM,

NSFLD, NUFLD, NDFLD,

PTOP, NULEV, (ZLEV(K),K=I,NULEV),

NLAY, (SIGE(K),K=I,NLAY+I),

NDUM, (DDUM(N),IDUM(N),RDUM(N),CDUM(N),N=I,NDUM)

READ(KU) (NAMES(N), DESCS(N), N=I,NSFLD ),

77
(NAMEU(N), DESCU(N), N=I,NUFLD ),

(NAMED(N), DESCD(N), NDLEV(N), N=I,NDFLD )

DO N=I,NSFLD

READ(KU) ZLEVID,NAME,Q

ENDD0

C ___________

C **** READ 3-DIMENSIONAL FIELDS (ORDERED ALL FIELDS FOR EACH LEVEL)  

DO L=I,NULEV

DO N=I,NUFLD

READ(EU) ZLEVID,NAME,Q

ENDD0

ENDD0

C ** $$ READ 3-DIMENSIONAL FIELDS (ORDERED ALL LEVELS FOR EACH FIELD) **_

DO N=I,NDFLD

DO L=I,NDLEV(N)

READ(KU) ZLEVID,NAME,Q

ENDD0

ENDD0

500 CONTINUE

STOP

END

78
PARAMETER TYPE DESCRIPTION

JOB Character*8

XLAB Character*80

NYMD Integer

NHMS Integer

NYMD0 Integer

NHMS0 Integer

IM Integer

JM Integer

NLAY Integer

SIGE Real

PTOP Real

ZLEV Real

NSFLD Integer

NUFLD Integer

NDFLD Integer

NAMES Character*8

DESCS Character*40

NAMEU Character*8

DESCU Character*40

NULEV Integer

NAMED Character*8

DESCD Character*40

NDLEV Integer

NDUM Integer

DDUM Character*40

IDUM Integer

RDUM Real

CDUM Character*8

PLEVS Real

ZLEVID Real

NAME Character*8

Q Real

Experiment Identifier (Usually a name or experiment number)

User Assigned Description of Current Experiment

Current Date of Experiment in YYMMDD format

Current Time of Experiment in HHMMSS format

Beginning Date of Experiment in YYMMDD format

Beginning Time of Experiment in HHMMSS format

Longitudinal dimension of data

Latitudinal dimension of data

Number of GCM vertical levels used in Experiment

GCM Sigma-Edge values

Model Top Pressure used in Experiment

Upper-Air Level Values archived on data (Pressure or Sigma

Level)

Number of 2-Dimensional Fields written in Config. #1

Number of 3-Dimensional Fields ordered by Level written in

Config. #2

Number of 3-Dimensional Fields ordered by Quantity written

in Config. #3

Names of 2-Dimensional Fields defined by NSFLD

Description of 2-Dimensional Fields defined by NSFLD

Names of 3-Dimensional Fields ordered by Level defined by

NUFLD

Description of 3-Dimensional Fields ordered by Level defined

by NUFLD

Number of Upper-Air Levels associated with all fields in Config.

#2

Names of 3-Dimensional Fields ordered by Quantity defined by

NDFLD

Description of 3-Dimensional Fields ordered by Quantity defined by NDFLD

Number of Upper-Air Levels associated with each field in Config. #3

Number of User-Defined Header variables. NDUM, DDUM,

IDUM, RDUM, and CDUM may be used by Users to add

REAL, INTEGER, or CHARACTER information in the Header

concerning the User's experiment.

User-Defined Description of Header variables

User-Defined Integer Header variables

User-Defined Real Header variables

User-Defined Character Header variables

Levels used for Sigma to Pressure Interpolation to create Pressure Level Output. (Maximum number of levels = 50)

Level Indentifier (Pressure or Sigma Level) of Quantity being

read

Character Name of Quantity being read

Two-Dimensional Quantity being read

79
It should be noted that the Datasets produced for the GEOS-1 DAS 5-year re-analysis and

the GEOS-1 GCM AMIP simulation have the following characteristics and User-Defined

Header information:

The horizontal resolution for all data is constant for all time periods.

All horizontal data, q(i, j), is written with i = I corresponding to Longitude =

-180 (International Dateline), and j = 1 corresponding to Latitude = -90

(South Pole).

The data written in Configuration # 1 are generally surface data (such as

Sea-Level Pressure, Surface Topography, Land/Water Masks, etc.).

The data written in Configuration # 2 are generally Upper-Air prognostic

quantities (such as wind, temperature, moisture, etc.).

The data written in Configuration :#: 3 are generally GCM Model Diagnostic

quantities (such as Diabatic Heating). Some Diagnostics are only defined at

one level (such as Precipitation, Outgoing Longwave Radiation, etc.).

The User-Defined Header information has been defined as:

NDUM =

DDUM(1) =

IDUM(1) =

lZDUM(1) =

CDUM(1) =

DDUM(2) =

IDUM(2) =

RDUM(2) =

CDUM(2) =

'OUTPUT RECORD TYPE (INTEGER)'

-1 FIRST GUESS

0 AFTER ANALYSIS (or STRAIGHT FORECAST)

1 INITIALIZED ANALYSIS

0.100000E+16

'LOGSR'

'UNDEFINED FILL VALUE (REAL)'

0

0.100000E+16

'UNDEF'

80
11.3 Geopotential Heights on Sigma Surfaces

In the GEOS-1 GCM, the pressure gradient computed on sigma surfaces is derived by

removing the large-scale component due to a mean potential temperature and surface pressure. Therefore, the perturbation geopotential heights have been archived on the GEOS-1

GCM sigma level datasets.

The pressure gradient on a constant pressure surface is related to the gradient on a constant

sigma surface by

VvO = V¢O + cpOV_p _ , (25)

while the hydrostatic equation is given as

- %0 •

Op_

Defining _ = • + _' and O = 0 + 0', we may write

(26)

(27)

or

- cp-O _ -_=-Cp-O(p'_-p'_) (28)

Op ,¢

- CpO'. (29)

Op_

Thus, we may express the pressure gradient on sigma surfaces as

Va¢ + cpOV_p _ :

= V-_ + V_' + %-OVp '_+ %0'Vp '_

= -Cp-_Vp '_ + V¢t + cp-OVp_ + cpS'Vp _

= V_' + %0'Vp '_ (30)

81
where

m

(I)_= (I)e - (I)e = v_ - _fi(p_ - p[)

; O= 280/1000 _

Po = 1000. (31)

82
11.4 Interpolation from Sigma to Pressure Surfaces

While the GEOS GCM uses the sigma coordinate in the vertical for the discritization of

the primitive equations, the current analysis system requires "first guess" fields to be on

prescribed pressure surfaces. In addition, it is often desirable to process model output on

pressure surfaces for scientific studies and comparisons to observations. The method used

for the sigma coordinate to pressure coordinate interpolation is consistant with the form of

the hydrostatic equation used in the GEOS GCM.

The vertical grid and index notation for sigma levels used in the GEOS GCM is given by

.................. _¢__PL..O..q. ............ e- 1

q,

ge

.................. ¢....PZ.q..q_............ e

Figure 6: Vertical placement and index notation for sigma levels in the GEOS GCM

The form of the hydrostatic equation used in the GEOS GCM is given by

Op,_ - -cpe (32)

where _ is the geopotential height, and _ = T/P '_. Similarly, for an arbitrary quantity, q,

we may define

Oq

-:C

where c is a constant. Thus, we may write

(33)

or

Oq

--C

qe- 1 - qg qg.- 1 - q,

P_-I - P_ P_-I - P,_

(34)

83
q* : qg--1 [P;-'- P:] (35)

+ PFj

This method is used for the winds, temperature, and relative humidity. Specific humidity

at pressure levels is obtained from the previously interpolated temperature and relative

humidity. For target pressures above the first model level, a linear extrapolation in P_

from levels 1 and 2 is used. For target pressures between the surface pressure and the

lowest model level, the quantity is assumed constant and is equal to its value in the lowest

model level. For target pressures greater than the surface pressure (below topography), the

quantity is undefined. It should also be noted that all GEOS GCM diagnostic quantities

on pressure surfaces have been similarly linearly interpolated in P_.

The method used in the interpolation of geopotential heights, in particular, must be intimately linked to the finite-differencing scheme used in the background model in order to

minimize computational error and prevent spurious height gradients (or geostrophic imbalances). The technique described here, formulated at UCLA by Arakawa and Suarez (personal communication), uses two separate interpolation algorithms depending on whether

you are above or below the nearest sigma-level interface, ge, and has the following properties:

@ P'_ O

Q

_-1

_e

Figure 7: Index notation for vertical interpolation of Geopotential Heights

1) _p_ = _ for P1 = P_

2) _P2 : _g-1 for P2 -- Pe-1

3) OP1 = OP2 for PI = P2 =/5_e

4) Vp_p1 = Va g2i q- CpOiVaP_ for P1 = Pl

5) VpCbp2 = Va¢_-_ + CpOe-l V,TP__l for P2 = P/-1

84
Thus, the interpolation not only preserves the heights on sigma surfaces whose target pressure coincides with the sigma level pressure, but also the horizontal gradient of the height

on pressure surfaces is defined accurately as the horizontal gradient of the sigma level height

plus the term corresponding to the horizontal gradient of the surface pressure.

We begin by writing the hydrostatic equation in the form

0¢

- %6. (36) OP"

Following Arakawa and Suarez (1983), we may write

(_g-1 = (_g2V£pO_e(P_-/P_-I)' (37")

tge_ = (p_ - Pt_-l)0e-1+ (P_ - P_)Oe (38)

(P?- PL1)

Thus, for the two regions in Figure 7, we have:

1) ForP, > P1 > Pee

2) For/Se e > P2 _ Pg-I

OP1 = _e + Cp 01dP _ (39)

where

OP2 = Og-1 + Cp 02dP '_

--1

(40)

/,?i

?'j,

OI(P_) = (P_-P_)Oee-I-(P'_-D_)Oe

(p__ p_) (40

(P_ - P')O,.l + (P" - PL1)&_ O2(P_) = (42)

Here OI(P a) and O2(P a) are the interpolated 0 values from levels g and e- 1 to the pressure

P. Using (42)in (39), we may define

85
or

I1 = 01dP '_

fp_(P_ P_)G + (P_ ^_

- -P_)Oedp,_

(Pi - G)

(43)

t¢ Ate (pl -Pio)

(44)

(45)

_ 1

-- 2

(46)

Similary, we may define

or

--1

--1

02dP _

(P[_- P_) o___+ (P_ - P__x)G dp ,_ (47)

A

(P_e-P?-I) [0g__ 1 (p_t_[p_ _ P?-l]-2[P_ 2 - P?-12])

--_O_e (l[p_2 - P__12]- IP?_I[P_- P?_I])]

(48)

86
= \ _-_-1 _[P_ +

(49)

(5o)

Therefore, for the interpolation of heights within regions (1) and (2) from Figure 7, we may

write:

1) For P_ > P1 _> /_ee

_P1 = _ + -_ \ p_ p_]

(51)

2) For/Stc _> P2 _> P_-I

Cp ( P_ - r_-l ) log_ 1 ( 2 _e _ g_ _ p__ l ) .__ _ge ( g_ _ g__ l ) ]

(52)

;, ::L _ ,.. ................................................................

87
11.5 Sea-Level Pressure

In addition to the archiving of surface pressure, which is a prognostic variable in the GEOS1 GCM, on all GEOS Production output streams, the sea-level pressure is also computed

and archived. Also, during the GEOS data assimilation cycle, sea-level pressure rather than

surface pressure is analyzed. The method used to compute sea-level pressure from surface

pressure involves integrating the hydrostatic equation from sea-level to the surface, and an

extrapolation of the temperature profile below topography.

The hydrostatic equation may be written as:

0p _ pg (53)

Oz Pg - - RT

Integrating (53) from sea-level, sl, to the surface, s, we may write

f s dp = _ fs s dO (54) t P tRT

where we have used dO = g dz. Assuming a mean temperature, T, between the sea-level

and surface, (54) may be written as

or

Ps ¢_ (55)

P_l = P_e¢_/nT • (56)

The mean temperature T is defined as the average of the surface and sea-level temperatures:

-- 1 T

T = _ ( _+ Ts_) (57)

where T_l is obtained by assuming the temperature follows the moist adiabatic lapse rate

dT /3_ _

__ = -/3 _ T_l= T_+ -- (5S)

dz g

where/3 to K The surface temperature used for the sea-level pressure calculation is.... _-_.

obtained by averaging the potential temperature, 0, in the lowest 100 mb of the model.

88
This allowsvariableresolutionchangesin the PBL whilemaintainingconsistentsea-level

pressurevaluesbetweenresolutions.Wehave

where

suchthat

Ts =-Op'_ ; 0 = T/p '_ (59)

_ E_ 0e,5_r_ (60)

E_ A,rl

100 mb

Aae- EeApe _ -- (61)

Here zr = Ps - Prop.

)(i,!i_

89
12 GEOS Unix Script

The GEOS GCM is currently being run on NASA/GSFC's Cray C-90 computer system

under the Cray UNICOS environment. There are five essential components to the basic

script which are needed in order to perform experiments with the GEOS GCM:

Create GEOS GCM NAMELIST for model parameters and diagnostics

Create DIAGSIZE Subroutine for enabled diagaostics

Remote copy frozen GEOS Libray, Application and Output Routines

from geos_das@dao

Compile User Application and Output Routines

Load and Run GEOS simulation

Each of these elements are briefly discussed in the following subsections.

12.1 Creating the GEOS GCM Namelist

The following example shows the creation of a typical GCM namelist to be used during a

GEOS simulation:

##################################################################

######

###### Setup GCM Namelist

######

##################################################################

/bin/rm gcmrun.namelist

cat > gcmrun.namelist <<'EOF'

kINPUT

JOB 1234',

XLAB = 'E1234 GEOS-GCM

'RIMENT

NYMDB = 850101

NHMSB = 000000

NYMDE = 860101

NHMSE = 000000

NDT = 225

LEAP = . TRUE.

MATS = . FALSE.

CLIMATE SIMULATION EXPE',

9O

2  ̧
ALPHA

NDSEG

NDOUT

PLEVS

= 0.05 ,

= 2400000 ,

= 060000 ,

= I000.0,

700.0,

200.0,

_END

_INDIAG

NUFLUX = 030000

NVFLUX = 030000

NHFLUX = 030000

NMOISTT =-060000

NMOISTQ =-060000

NRADLW =-060000

NRADSW =-060000

NPREACC = 060000

NPRECON : 060000

NOLR = 060000

NEVAP = 060000

NU2H = 010000

NV2M = 010000

NT2M = 010000

NQ2M = 010000

_END

'EOF'

950.0 , 900.0, 850.0, 800.0,

600.0 , 500.0, 400.0, 300.0, 250.0,

150.0 , 100.0, 70.0, 50.0, 30.0, 20.0

For a complete description of GCM Namelist parameters, see Section 11.1.

12.2 Creating th_ DIAGSIZE Subroutine

Using the DIAGSIZE subroutine is a convenient method to allocate just enough permanent memory during a run to be used for diagnostic storage. All GEOS GCM diagnostic

quantities are stored in the single diagnostic array QDIAG which is located in the DIAG

COMMON. In order to minimize the maximum memory requirement used by the model,

the GEOS GCM libraries are compiled with room for only one horizontal diagnostic array.

All production DAO applications contain a call to subroutine DIAGSIZE, which expands

the size of the DIAG COMMON to include the total number of surface and upper-air diagnostics the User has enabled. For the diagnostics enabled in the preceeding section, the

appropriate DIAGSIZE subroutine should be defined as:

##################################################################

######

###### Create DIAGSIZE Routine to control Total Diagnostics

91
######

##################################################################

/bin/rm diagsize.f

cat > diagsize.f <<'EOF'

subroutine diagsize

parameter ( imr = 144 )

parameter ( jmr = 90 )

parameter ( nlayr = 20 )

parameter (ndiags = II )

parameter ( ndiagu = 4 )

common /diag/ ndiagmx,qdiag(imr,jmr+l,ndiags+ndiagu*nlayr)

ndiagmx = ndiags + ndiagu*nlayr

return

end

'EOF'

We see that the User needs to simply assign the appropriate values for the number of

surface diagnostics (NDIAGS) and the number of upper-air diagnostics (NDIAGU) which

are enabled in the User's n_melist.

12.3 Remote copy frozen GEOS Library

In order to gain access to the GEOS GCM, the User must remote copy the frozen production

system from the DAO production workstation, geos_das@dao. In addition to the GEOS

GCM library, the User m_y _lso remote copy production application programs and output

routines. The COMMON files associated with these programs must also be copied to the

User's local directory. Finally, the GEOS GCM boundary conditions and other GEOS

input datasets must be obtained from the production workstation. An example of the

syntax necessary for these remote copies is as follows:

##################################################################

######

###### Get Frozen GCM, Applications, and Output Koutines

######

###############$##################################################

rcp geos_das_dao:/production/geos_das/geosl.1/gcm/vc5.2b20.macro vcS.2620

rcp geos_das@dao:/production/geos_das/geosl. I/applications/gcmprodn.f.

rcp geos_das@dao:/production/geos_das/geosl. I/applications/progprs.f •

rcp geosdas©dao:/production/geos_das/geosl.1/applications/progsig.f •

92
rcp geos_das@dao:/production/geos_das/geosl. I/applications/diagprs.f .

rcp geos_das©dao:/production/geos_das/geosi. I/applications/vrslvb20.com vrslv.com

rcp geos_das@dao:/production/geos_das/geos1.1/applications/vstate.com .

rcp geos_das@dao:/production/geos_das/geosl. I/applications/vreal.com .

rcp geos_das@dao:/production/geos_das/geos1.1/applications/vphyscom.com .

rcp geos_das@dao:/production/geos_das/geosl.i/applications/vmetric.com .

rcp geos_das%dao:/production/geos_das/geos1.1/applications/metricc.com .

rcp geos_das_dao:/production/geos_das/geos1.1/applications/vbudget.com .

rcp geos_das@dao:/production/geos_das/geosl. I/applications/pblparm.com .

rcp geos_das@dao:/production/geos_das/geosl. I/applications/iau.com .

rcp geos_das_dao:/production/geos_das/geosl.2/data/bc14409I.y1985 gcmbcs

rcp geos_das@dao:/production/geos_das/geosl.2/data/gcmo3.data .

rcp geos_das@dao:/production/geos_das/geosl.2/data/zonalT.data .

12.4 Compile User Application and Output Routines

Since the GEOS GCM library only contains subroutines from GCMDRV and below, the User

must provide compiled Main Programs, or Applications, to perform specific experiments.

For complete consistancy with the GEOS library, the Application and Output routines

should be compiled in the same manner as the background model. An inline compile routine

to accomplish this is shown below:

##################################################################

######

###### Create COMPILE Routine and Compile Userappl _ Output

######

##################################################################

/bin/rm compile

cat > compile <<'EOF'

/lib/cpp -P -I . $I.f $1.cpp

fpp -dec $1.cpp > $1.fpp

cft77 -i 64 -a static -V -exm $1.fpp

my $1.fpp.o $1.o

'EOF'

chmod +x compile

./compile gcmprodn

./compile progprs

./compile progsig

./compile diagprs

./compile diagsize

93
Here we have compiled those routines which were remote copied in Section 12.3.

12.5 Load and Run GEOS Simulation

Following successful compilations, the User must load the Application and Output routines

together with the GEOS GCM library, creating an executable a.out. Then, assuming that

an appropriate GEOS GCM Restart has been brought to the current local directory as

gcmrst, the GEOS simulation may be performed:

##################################################################

######

###### Load GCM, Applications and Output Routines, and Run

######

##################################################################

segldr -M,s -i ./vcS.2620 -- gcmprodn.o progsig.o progprs.o diagprs.o diagsize.o

assign -a gcmrun.namelist fort.ll

assign -a gcmrst fort.12

assign -a prog.sig

assign -a prog.prs

assign -a diag.prs.a

assign -a diag.prs.b

assign -a diag.prs.c

-N ieee -F f77 fort.34

-N ieee -F f77 fort.35

-N ieee -F f77 fort.41

-N ieee -F f77 fort.42

-N ieee -F f77 fort.43

assign -a gcmbcs -N ibm -F COS fort.16

assign -a gcmo3.data fort.20

assign -a zonalT.data fort.52

setenv NCPUS 4

a.out

Note that NCPUS has been set as an environmental variable equal to 4 in order to take

advantage of the macro-tasking properties of the GEOS GCM. In this example, file units 34,

35, 41, 42, 43 are associated with Output data streams, and file units 16, 20, 52 are input

datasets and boundary conditions for the GEOS GCM. Once the GEOS GCM simulation

has been completed, the resulting output data streams should be moved to permanent

storage.

94
Acknowledgements

The model development and research documented in this report were supported by the EOS

Interdisciplinary Science Program and by the NASA Research and Applications Program.

Computer resources and funding were provided by the EOS Project through the Scientific

Computational Facility of the Data Assimilation Office.

95

References

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. Wen. Rev., ii1, 34-45.

Asselin, R., 1972: Frequency fiIter for time integrations. Mon. Wea. Rev., i00, 487-490.

Blackadar, A.K.,1977:: High Resolution Models of the Planetary Boundary Layer. Advances

in Environmental Science and Engineering, Vol i Editors Pfaffiin and Zeigler, Gordon and Breach, Scientific Publishers.

Bloom, S.C., L.L. Takacs, and E. Brin, 1991: A scheme to incorporate analysis increments

graduMly in the GLA assimilation system. Ninth Conf. on Numerical Weather

Prediction, Denver, CO, Amer. Meteor. Soc., I10_I12.

Burridge, D.M. and J. Haseler, 1977: A model for medium range weather forecasting--

adiabatic formulation, Tech. Report. No, 4, European Center for Medium Range

Weather Forecasts, Brachnell, Berkshire, UK.

Chou, M.-D., 1984: Broadband water vapor transmission functions for atmospheric IR flux

computation. J. Atmos. Sci., 41, 1775-1778.

Chou, M. -D., and L. Peng, 1983: A parameterization of the absorption in 15-micron C02

spectral region with application to climate sensitivity studies. J. Atmos. Sci:, 40,

2183-2192.

Clarke, R. H., 1970: Observational studies in: the atmospheric boundary layer. Quart. J.

Roy. Meteor. Soc., 96, 91-114.

Davies, R., 1982: Documentation of the solar radiation parameterization in the GLAS

climate model. NASA Tech. Memo. 83961, 57 pp. [Available from the U.S. Department of Commerce,• NationM Technical Information Service, 5285 Port Royal Road,

Springfield, VA 22161]

Dorman, J. L., and P. J. Sellers, i989: A global climatology of albedo, roughness length

and stomatal resistance for atmospheric general circutation models as represented by

the Simple Biosphere model (SIB). J. Appl, Meteor., 28, 833-855.

Fox-Rabinovitz, M., It. M. Helfand, A. Hou, L. L. Takacs, and A. MoIod, 1991: Numerical

experiments on forecasting, climate simulation and: data assimilation with the new

17 layer GLA GCM. Ninth Con:ference on Numerical Weather Prediction. 21-25

\

PAGE BLAb_ NOT FR.MED

97

.__ _!
October1991,Denver,CO,506-509.

Gates,W. Lawrence,1992:AMIP: TheAtmosphericModelIntercomparisonProject. Bull.

Am Met. Soc., 73, 1962-1970

Harshvardhan, R. Davies, D. A. Randall, and T. G. Corsetti, 1987: A fast radiation parameterization for atmospheric circulation models. J. Geophys. Res., 92, 1009-1016.

Helfand, H. M., and J. C. Labraga, 1988: Design of a non-singular level 2.5 second-order

closure model for the prediction of atmospheric turbulence. J. Atmos. Sci., 45,

113-132.

Helfand, H. M., M. Fox-Rabinovitz, L. Takacs, and A. Molod, 1991: Simulation of the

planetary boundary layer and turbulence in the GLA GCM. Proceedings of the AMS

Ninth Conference on Numerical Weather Predicition, 21-25 October 1991, Denver,

CO, 514-517.

ttelfand, H. M., M. and S. D. Schubert, 1994: Climatology of the simulated Great Plains

low-level jet and its contribution to the continental moisture budget of the United

States. J. Climate, submitted.

Higgins, R. W., and S. D. Schubert, 1993: Low-frequency synoptic-eddy activity in the

Pacific storm track. J. Atmos. Sci., 50, 1672-1690.

Kalnay, E., M. Kanamitsu, J. Pfaendtner, J. Sela, M. Suarez, J. Stackpole, J. Tuccillo, L.

Umscheid, and D. Williamson, 1989: Rules for the interchange of physical parameterizations, Bull. Am Met. Soc., 70, 620-622.

Kondo, J., 1975: Air-sea bulk transfer coefficients in diabatic conditions. Boundary Layer

Meteor., 9, 91-112.

Lacis, A. A., and J. E. Hansen, 1974: A parameterization for the absorption of solar

radiation in the Earth's atmosphere. J. Atmos. Sci., 31, 118-133.

Lanczos, C., 1966: Discourse on Fourier Series. Hafner Publishing, 255 pp.

Large, W. G., and S. Pond, 1981: Open ocean momentum flux measurements in moderate

to strong winds. J. Phys. Oceanogr., 11,324-336.

Meador, W. E., and W. R. Weaver, 1980: Two-stream approximations to radiative transfer

in planetary atmospheres: A unified Descriptoin of existing methods and a new

improvement. J. Atmos. Sci., 37, 630.643.

Mintz, Y. and Y. Serafini, 1984: Global fields of monthly normal soil moisture as derived

from observed precipitation and an estimated potential evapotranspiration. Final

98
scientificreport underNASA grantNAS5-26,Part V, Dept. of Meteorology,Universityof Maryland,CollegePark,MD

Moorthi, S.,and M. J. Suarez,1992:RelaxedArakawaSchubert:A parameterizationof

moistconvectionfor generalcirculationmodels.Mon. Wea. Rev., 120, 978-1002.

Panofsky, H. A., 1973: Tower Micrometeorology. In Workshop on Micrometeorology, D. A.

Haugen (ed.), American Meteorological Society, Bosston, 392 pp.

Phillips, N. A., 1957: A coordintate system having some special advantages for numerical

forecasting. J. Meteor., 14, 184-185.

Platt, H. R., and Harshvardahn, 1988: Temperature Dependance of Cirrus Extinction:

Implications for Climate Feedback. J. Geoph. Res., 93, 11051-11058.

Posey, J. W., and P. F. Clapp, 1964: Global distribution of normal surface albedo. Geofis.

Int., 4, 33-48.

Reynolds, R. W., 1988: A real time global sea surface temperature analysis. J. Climate, 1,

75-86.

Rodgers, C. D., 1968: Some extension and applications of the new random model for

molecular band transmission. Quart. J. Roy. Meteor. Soc., 94, 99-102.

Rosenfield, J. E., M. R. Schoeberl, and M. A. Geller, 1987: A computation of the stratospheric diabatic circulation using an accurate radiative transfer model. J. Atmos.

Sei., 44, 859-876.

Sadourney, R., 1975: The dynamics of finite difference models of the shallow water equations, J. Atmos. Sci., 32,680-689.

Schemm, J., S. Schubert, J. Terry, and S. Bloom, 1992: Estimates of monthly mean soil

moisture for 1979-1989. NASA Tech. Memo. No. 104571, Goddard Space Flight

Center, Greenbelt, MD 20771.

Schubert, S. D., J. Pfaendtner and R. Rood, 1993: An assimilated data set for Earth Science

applications, Bull. Am Met. Soc., 74, 2331-2342

Schubert, S., M. Suarez, C. K. Park, and S. Moorthi, 1993: GCM Simulations of intraseasonal variability in the Pacific/North American region. J. Atmos. Sci., 50, 1991-

2007.

Shapiro, R., 1970: Smoothing, filtering and boundary effects. Rev. Geophys. Space Phys.,

8, 359-387.

99
Suarez,M. J., andL. L. Takacs,1994:Documentationof the Aries/GEOSDynamicalCore

Version2, NASA Tech.Memorandum104606-Vol(inpress),NASA,GoddardSpace

Flight Center,Greenbelt,MD.

Sud,Y. C.,andA. Molod,1988:Therolesof dry convection,cloud-radiationfeedbackprocessesandtheinfluenceofrecentimprovementsin theparameterizationof convection

in the GLA GCM. Mon. Wea. Rev., 116, 2366-2387.

Sud, Y. C., P. J. Sellers, Y. Mintz, M. D. Chou, G. K. Walker, and W. E. Smith, 1990:

Influence of the Biosphere on the Global Circulation and Hydrological Cycle - A

GCM Simulation Experiment. Agric. and Forest Meteo., 52, 113-180.

Yaglom, A. M., and B. A. Kader, 1974: Heat and mass transfer between a rough wall

and turbulent fluid flow at high Reynolds and Peclet numbers. J. Fluid Mech., 62,

601-623.

100

•

__ -_I. i____, I_ _i_ _I_  ̧ ̧,_

ri_-I_i__  ̧ ̧I__•

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for 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

September 1994 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

Technical Report Series on Global Modeling and Data Assimilation

Volume 1 - Documentation of the Goddard Earth Observing System

(GEOS) General Circulation Model - Version 1

6. AUTHOR(S)

Lawrence L. Takacs, Andrea Molod, and Tina Wang

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

94B00128

10. SPONSORING/MONITORING

ADGENCY REPORT NUMBER

TM- 104606, Vol. 1

11. SUPPLEMENTARYNOTES

L. Takacs, A. Molod, and T. Wang: General Sciences Corporation, Laurel, Maryland (at GSFC Data

Assimilation Office, Greenbelt, Maryland)

M. Suarez: Goddard Space Flight Center_ Greenbelt_ Maryland

12a. DISTRIBUTION/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

This technical report documents Version 1 of the Goddard Earth Observing System (GEOS) General

Circulation Model (GCM). The GEOS-1 GCM is being used by NASA's Data Assimilation Office

(DAO) to produce multiyear data sets for climate research. This report provides a documentation of

the model components used in the GEOS- 1 GCM, a complete description of model diagnostics

available and a User's Guide to facilitate GEOS-1 GCM experiments.

14. SUBJECT TERMS

GEOS, EOS, GCM, Documentation, Diagnostics, User's Guide, DAO

17. SECURITY CLASSIFICATION

OF REPORT

Unclassified

NSN 7540-01-280-5500

18. SECURITY CLASSIFICATION

OF THIS PAGE

Unclassified

19, SECURITY CLASSIFICATION

OF ABSTRACT

Unclassified

15. NUMBER OF PAGES

114

16. PRICE CODE

20. LIMITATION OF AB_ I RAC'f

UL

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39.18

298-102

. : i. .,,.