NASA Technical Memorandum 104606, Vol. 11

Technical Report Series on

Global Modeling and Data Assimilation

Max J. Suarez, Editor

Volume 11

Documentation of the Tangent Linear and

Adjoint Models of the Relaxed

Arakawa-Schubert Moisture

Parameterization Package of the NASA

GEOS-1 GCM (Version 5.2)

Weiyu Yang, I. Michael Navon, and Ricardo Todling

May 1997

DATA ASSIMILATION OFFICE
NASA Technical Memorandum 104606, Vol. 11

Technical Report Series on

Global Modeling and Data Assimilation

Max J. Suarez, Editor

Goddard Space Flight Center

Greenbelt, Maryland

Volume 11

Documentation of the Tangent Linear and

Adjoint Models of the Relaxed

Arakawa-Schubert Moisture

Parameterization Package of the NASA

GEOS-1 GCM (Version 5.2)

Weiyu Yang

I. Michael Navon

Supercomputer Computations Research Institute

Florida State University

Tallahassee, Florida

Ricardo Todling

General Sciences Corporation

Laurel, Maryland

National Aeronautics and

Space Administration

Goddard Space Flight Center

Greenbelt, Maryland

1997
Abstract

A detailed description of the development of the tangent linear model (TLM) and its

adjoint model of the Relaxed Arakawa-Schubert moisture parameterization package used

in NASA GEOS-1 C-Grid GCM (Version 5.2) is presented. The notational conventions

used in the TLM and its adjoint codes are described in detail.

iii
Contents

List of Figures vii

1 Introduction 1

Description of the Moisture Parameterization Scheme 1

2.1 Cumulus Parameterization Package of the NASA GEOS-1 GCM ...... 5

2.2 Description of the Discretization of the Moisture Physics Parameterization . 12

Tangent Linear Model of the Moist Process Physics Package 17

3.1 Linearized Discrete Dynamical Equations ................... 23

3.2 Coding of the Tangent Linear Model ...................... 23

3.3 Notational Convention for Variables and Subroutines Used in the Tangent

Linear Model Code ................................ 24

4 Adjoint Model of the Moist Process Physics Package

4.1 Using the Adjoint Method to Calculate the Gradient of a Cost Function . .

4.2 Coding of the Adjoint Model ..........................

4.3 Notational Convention for Variables and Subroutines Used in the Adjoint

Model Code ....................................

4.4 Verification of the Correctness of the Adjoint Model .............

Acknowledgments

References

24

24

25

26

26

31

33
List of Figures

Variation of the _(_) with respect to log g (gradient check of correctness of

adjoint model). Integration period is 6 hours and t = 6 hours model generated

observations were used. January 1, 1985 00Z DAO's data was used as t = 0

observations. The first guess is the shifted 6-hour initial condition. The time

integration scheme employed is the leapfrog scheme ..............

As in Figure 1, but for the Matsuno time integration scheme .........

Variation of the log I_(g) - 11 with respect to log g. Integration period is 6

hours and t = 6 hours model generated observations were used. January 1,

1985 00Z DAO's data is used as t = 0 observations. The first guess is the

shifted 6-hour initial condition. The time integration scheme employed is the

leapfrog scheme ..................................

As in Figure 3, but for the Matsuno time integration scheme .........

28

28

29

29

vii
1 Introduction

The GEOS-1 C-Grid GCM was developed by the Data Assimilation Office (DAO) at Goddard Laboratory for Atmosphere (GLA), NASA,/GSFC (Takacs, et al., 1994; Suarez and

Takacs 1994) to be used in conjunction with an analysis scheme to produce a multi-year

global atmospheric data set for climate research (Schubert et al., 1993). It has an advanced

structure, i.e., a "plug-compatible" structure. It means that if "plug-compatible" rules are

followed in coding different GCMs and parameterizations, codes can be "unplugged" from

one model and "plugged" into another with little coding effort. Thus each part of the

GEOS-1 C-grid GCM can be used independently at another GCM. For instance, the full

physics package of GEOS-1 C-grid GCM has been used into NASA/GLA Semi-Lagrangian

Semi-Implicit (SLSI) GCM. The DAO and the Climate and Radiation Branch at GLA,

NASA/GSFC have produced a library of physical parameterizations and dynamical algorithms which may be utilized for various GCM applications.

There are four physics parameterization packages in the GEOS-1 GCM, i.e., the Relaxed

Arakawa-Schubert (RAS) subgrid moisture parameterization and large-scale convection

schemes, the shortwave radiation scheme, the longwave radiation scheme, and the turbulence parameterization scheme. Recently, the land surface model (LSM) described by

Koster and Suarez (1992) has been fully coupled to the Aries GCM at NASA/GSFC, and

the coupled models have been used to address a number of climate-related problems (Koster

and Suarez 1996). In the physics packages, the moisture process plays an essential role towards improving the quality of the GCM productions.

In order to obtain a 4-D variational data assimilation system based on the NASA GEOS-1 Cgrid GCM including the moisture parameterization process, we need to develop the tangent

linear and corresponding adjoint model of this physics package. This document describes the

development of the tangent linear and adjoint models of the RAS and large-scale convection

schemes as well as the formation of clouds which provides the cloud information used for

cloud-radiative interactions.

In Section 2 we provide a condensed description of the moisture parameterization package,

including the basic original equations and their discrete forms. Section 3 describes and

documents in detail the derivation and coding at the tangent linear model (TLM), and the

notational conventions used. Section 4 describes the derivation of the adjoint model code

and discusses its validity.

2 Description of the Moisture Parameterization Scheme

In a large-scale disturbance, there are many individual cumulus clouds whose spatial and

temporal scales are much smaller than the disturbance itself. Because of this scale separa-
tion, it may be possible to predict the time evolution with the collective influence of all the

smaller scale cumulus clouds. This is the goal of cumulus parameterization.

Subgrid scale cumulus convection parameterization is a crucial component of any GCMs.

Without including it in a GCM, one cannot hope to simulate the right general circulation patterns especially over the tropical region, and some other regions, where cumulus

convective activities are frequent and strong.

The necessity of parameterization of cumulus convection was understood and became clear

after early research efforts failed to explain theoretically the size and growth rate of tropical

cyclones (Lilly 1960; Yanai 1964). Among of the early parameterization researches, Charney

and Elisaaen (1964) and Ooyama (1964) presented their classical papers in which the concept

of %onditional instability of the second kind _ (CISK) was first introduced. CISK describes

the cooperative interaction between the cumulus scale and large scale, i.e., the large-scale

circulation is responsible for organizing and maintaining cumulus convection by providing

the necessary horizontal transport of water vapor, while the cumulus-scale drives the largescale circulation through the release of latent heat in deep convective elements. CISK

mechanism treats cumulus activities to be a function of the large-scale fields and it is now

commonly referred to as cumulus parameterization. The simple parameterizations used in

these papers led to considerable success in the numerical simulation of tropical cyclones

(e.g., Ooyama 1969). Due to the high degree of empiricism and intuition, and the lack of

theoretical framework for describing the mutual interaction between a cumulus ensemble

and the large-scale environment, these early parameterizations were too crude to be used

in GCM.

Ooyama (1971) first proposed a cumulus parameterization theory which took into account

the coexistence of a spectrum of clouds. He assumed that cumulus clouds can be represented as non-interacting spherical bubbles, dispatched from the mixed layer. He concluded

that the problem of parameterization of cumulus convection reduces to a determination of

the dispatch function thus his parameterization scheme is not closed due to the unknown

dispatch function.

Before 1974, the basic physical image related to the cumulus parameterization, was that

an existing cumulus cloud ensemble produces time changes in the large-scale temperature

and moisture fields (Arakawa 1969, 1971, 1972; Betts 1973a, b; Gray 1972; L6pez 1972a, b;

Ooyama 1971; Yanai 1971a, b; Yanai et al. 1973). Cumulus convection modifies the largescale temperature and moisture fields through detrainment and cumulus-induced subsidence

in the environment. The detrainment causes large-scale cooling and moistening, and the

cumulus-induced subsidence causes large-scale warming and drying. Based on 1956 Marshall

Island observation data, Yanai et al. (1973) quantitatively derived these effects from a

combination of observed large-scale heat and moisture budget over an area covered by

cloud cluster and a cumulus ensemble model which is similar to the one used in Arakawa

and Schubert (1974). Their results show the importance of coexistence of shallow clouds
with deepcloudsin maintainingthe large-scaleheatand moisturebudgets.

Basedonthe researchresultsof the cumulusparameterizationmentionedabove,Arakawa

andSchubert(1974,hereafterAS) presenteda remarkableachievement,namelythey proposeda closedcumulusparameterizationtheory that describesthe mutualinteractionof

anensembleof cumuluscloudswith the large-scaleenvironment.This is perhapsthe most

physicallycompleteapproachto the issueof cumulusparameterizationup to now. The

cloudensembleis representedby a spectrumof idealizedmodelcloudsub-ensembles. Each

of the sub-ensembleshasits ownmass,heatand moisturebudget.Theverticaltransports

accomplishedby this ensembleof modelcloudsis ultimatelydeterminedby the cloudbase

massflux for eachmemberof ensemble.In orderto determinethis massflux distribution, ArakawaandSchubertintroducedtheconceptof quasi-equilibriumof the cloudwork

function,whichleadsto anintegralequationmakingthe cloudbasemassflux distribution

relatingto large-scalethermodynamicprocesses.

Kuo (1965,1974)introducedcumulusparameterizationfor usein tropicalcyclonemodeling

whichweresubsequentlyemployedin large-scalenumericalpredictionmodels(e.g.,Krishnamurti 1969;Krishnamurtiet al. 1979;the NMC spectralmodeland the limited-area

nonhydrostaticMM5 model).TheKuoparameterizationschemeusesdifferentmechanisms

to describephysicallythe scaleinteractionbetweenthe cumulusensembleand the environmentbasedon similar assumptionscomparableto the AS parameterizationscheme

(Fraedrich1973). The mechanismappliedin the AS schemeis the verticalmassflux relating the fluxesinsideandoutsidea convectiveelement,i.e.,the vertical masstransport

as the representativevariable. Yet the mechanismappliedin the Kuo schemeis based

on a non-steadydeepcumulusmodel,usingthe temperaturedifferencebetweenthe cumuluscloudandthe undisturbedenvironmentandthe large-scaleconvergenceof moisture

asindicators. The advantageof Kuo'sschemeis that, as a parameterizationprocedure,

it providesimmediatemeasuresof the cumulus-scaleheat and moisturefluxesin terms

of measurablelarge-scalevariables,without havingto computeclouddynamicalprocesses

(suchasentrainment,detrainmentanddowndrafts)andcloudmicrophysicsprocesses.

Manabeet al. (1965)introducedanadjustmentschemeto simulatesubgrid-scalemoistconvection.Theiradjustmentincludesthedry convectiveadjustmentandthe moistconvective

adjustment.The moistadjustmentis performedonly whenthe relativehumidity reaches

100%andthe lapserate exceedsthe moistadiabaticlapserate. With theseartificial adjustments,theyintroducedasimplifiedhydrologiccycle,whichconsistsof the advectionof

watervaporbylarge-scalemotion,evaporationfromthe surface,andprecipitationprocess,

to simulatethe processof moistconvection.

The basicclosureassumptionfor the AS parameterization,i.e., the cloud-workfunction

quasi-equilibriumassumption,wasexaminedby Lord and Arakawa (1980). They justified

this assumption by considering the kinetic energy budget of a cumulus subensemble. That

is, the generation and dissipation of kinetic energy per unit cloud-base mass flux should
approximatelybalanceovertime scalesof thelarge-scaleprocesses,andsuchkineticenergy

generation(the cloud-workfunction) and dissipationper unit cloud-basemassflux for a

givensubensembleshouldnot dependsubstantiallyon the large-scaleconditions. They

calculatedcloud-workfunctionsfrom a varietyof data setsin the tropicsand subtropics

includingthe GATE,AMTEX, VIMHEX aswellascompositedtyphoondata. Theirresults

confirmedthe correctnessof the the cloud-workfunction quasi-equilibriumassumption.

Lord (1982)continuedto verify the AScumulusparameterizationusinga semi-prognostic

approachappliedGATE PhaseIII data. His resultsshowthat the calculatedprecipitation

agreesvery well with estimatesfrom the observedlarge-scalemoisturebudgetandfrom

radar observations.The calculatedvertical profilesof cumuluswarmingand drying are

alsoquite similar to the observation.His resultsalsoshowthat the error causedby the

cloud-workfunctionquasi-equilibriumassumptionis generallylessthan 10%.In hispaper,

someadditionalexperimentswerealsocarriedout to investigatethe sensitivityof the AS

schemeto someof the arbitrary parametersandassumptions.After theseverificationsand

investigations,Lord et al. (1982)incorporatedthe discretizedform of the ASschemeinto

the UCLA GCM.

TheAScumulusparameterizationhasbeenwidelytestedandappliedfor variouspurposes.

Ramanathan(1980)appliedthe ASschemein a semi-prognosticasestudyof a monsoon

depressionwhenit wasforming.Hefoundthat quasi-equilibriumheldevenin thisdisturbed

situation,and precipitationand cumulusheatingrateswerereasonable.Krishnamurtiet

al. (1980)comparedfivecumulusparameterizationschemesusingthe semi-prognosticapproach.Thecalculatedrainfallratewerecomparedwith theobservedestimatesfrom GATE

A/B data. Moorthi andArakawa(1985)usedthe ASparameterizationto carryout a systematicinvestigationon howcumulusheatingaffectsbaroclinicinstability. KaoandOgura

(1987)testedthe ASschemethrougha semi-prognosticapproachusingthe tropicalcloud

banddata and the tropical compositeeasterlywavedisturbancedata. They found that

cloudheatingand drying effectsaswell as the predictedcloud populationagreedrather

wellwith the observations.Sudet al. (1991)testedthe ASschemeusingthe GLA/GCM.

They modifiedsomeparametersusedin the ASschemeto improveparameterization,and

foundthat theroleof cumulusconvectionin maintainingthe observedtropicalrainfalland

850mbeasterlywindswassatisfactory.

Dueto the complexityof the originalAScumulusparameterization,someresearchefforts

werecarriedout aimedat simplifiedparameterizationschemeswhich wouldstill provide

realisticvaluesof thethermalforcingby convectionundervarioussynopticconditions(e.g.,

Ceselski1974;Hacket al. 1984;Tiedtke 1989;Moorthi and Suarez1992). Hacket al.

(1984)useda convectiveflux form of the AS schemeinsteadof the originaldetrainment

form. This flux form is moreconveniento usein numericalweatherpredictionmodels.

Tiedtke(1989)introducedamuchsimplerparameterizationschemeandcomparedits results

with the conventionalconvectionschemeusedin NWP at ECMWF. In 1992,Moorthi

andSuarezintroducedthe "RelaxedArakawa-Schubert _ (RAS)cumulusparameterization

schemewhichis asimplifiedASparameterizationscheme.It is veryefficient,andproduces
results very close to those of the standard AS implementation. The RAS scheme is used

in the NASA GEOS-1 GCM and is described in next subsection, following closely Moorthi

and Suarez (1992).

2.1 Cumulus Parameterization Package of the NASA GEOS-1 GCM

The cumulus parameterization package of the NASA GEOS-1 GCM is the RAS scheme

proposed by Moorthi and Suarez (1992). RAS makes two major simplifications in the

standard AS implementation:

1) It modifies the entrainment relation to avoid the costly calculation that is required to

find the entrainment parameter of clouds detraining at each levels.

2) "Relaxes" the state toward equilibrium instead of requiring "quasi-equilibrium '_ of the

cloud ensemble to be achieved each time the parameterization is invoked.

The iteration method chosen in the RAS scheme considers one cloud type at a time and

computes the cumulus mass flux that would be required to maintain the invariance of the

work function if there were no other clouds present. It then uses the heating and drying

effects of the given type of cloud to change the large-scale environment fields; the same

thing is done for another cloud type based on the new modified environment fields. With

this procedure, each step is in the direction of a single-cloud equilibrium, but in the course

of iteration, all cloud types affect each other by modifying the environment.

a. Cloud modal

The cloud model used by RAS is a simplified form of that in the AS scheme. The first

major simplification introduced is to assume that the normalized mass flux for each cloud

type is a linear function of height instead of the exponential function used in AS. Thus

_ :, (1)

Oz

where _lx(z) is the normalized mass flux for cloud type A at height z, with boundary condition

= 1 (2)

The hydrostatic equation is used in the form

Oz _ cp0 (3)

OP g

where P = (p/po) R/cp, p is the pressure, R is the gas constant, cp is the specific heat at

constant pressure, g is the gravitational acceleration, 0 is the potential temperature, and
po = looo ,,_b.From Eqs. (3) and (1), we have

&>(P) _ cp e,_ (4)

OP g

and integration gives

,x(P) = 1+ _Ajp _ e dP (s)

g

where PB = P(zB), is the pressure at the cloud base height.

As in the AS scheme, the large-scale budget of moist static energy and total water substance

for each cloud type are

a Orlx(P)h(p ) (6) aP [r>(P)t_(P)]- aP

and

0

aP {r/x(P)[q_,(P) + I_,(P)]} - &Ix(P) aP q(P) (7)

where h_(P), q_,(P), and I_(P) are the cloud moist static energy, specific humidity, and

liquid water mixing ratio for cloud type ,_ at level P, respectively, h(P) and q(P) are the

moist static energy and specific humidity in the environment, respectively.

In Eq. (7), the precipitation term has been neglected for simplicity. It is assumed that all

liquid water is carried to the cloud top where part is precipitated and part is evaporated,

depending on the cloud type. Since the liquid water loading and the precipitation effects

are excluded from Eq. (7), there is no need to specify the vertical distribution of the

precipitation and of the liquid water.

With these assumptions, the detrainment level of the cloud type A is the level at which the

moist static energy within clouds equals the saturation moist static energy of the environment. That is

where PD = PD (A) is the detrainment level.

Integrating Eq. (6) from PB to PD, we obtain

,x(P>)t4(P>) - t_B= -_;_g Pjp:et_(P)dP (9)

Combining Eqs. (5) and (9), Eq. (8) can be solved directly for the value A corresponding

to clouds that detrain at a given level PD:

(lo)
Similarly, assuming the cloud air is saturated at the level of non-buoyancy,

q_(r_) = q*(r_) (11)

and integrating Eq. (7) from PD to PB, the liquid water mixing ratio at the detrainment

level I(PD) is calculated from

1

l(P_) -- l_(P_)= ,_x(P_)

b. Cloud work" function

Following AS and neglecting the effects of water vapor and liquid water on the buoyancy,

the cloud work function Ax is

_ZD g Ax = B Cp_(Z)'lX(Z)[_I(Z)- _(_)] dz (13)

for cloud type A, where T(z) is the temperature in the environment at height z, and s_(z)

and s(z) are the cloud's and the environment's dry static energies, respectively.

Using the hydrostatic equation (3), we have

Ax = j;% ,Ix(P)[_i(P)p- _(P)] dP (14)

D

To obtain the cloud work function in terms of the moist static energy, let us approximate

the static energy difference as

1

st(P) - s(P) _ 1 + _(P)[hi(P) - h*(P)] (15)

where G(P) = (L/cp)[dq*(P)/dT], L is the latent heat of condensation of water vapor and

h* and q* are the saturation moist static energy energy and the saturation specific humidity

of the environment. Finally the work function form used in the GCM is

Ax= j;%D 1T-EP) 'Ix(P) tffx(P)-h*(P) P dP (16)

c. Cumulus efJ)cts on the lar'ge-scale budgets

The rate of change of dry and moist static energies due to cumulus convection can be written

in the form

= g_ gLD(P)I(P)[1- _.(P)] (17)

op
and

""(Oh) 57 V Oh _ gD(P)(,_ is) (18) c= v cot

where M_(P) is the total cumulus mass flux per unit horizontal area at level P, D(P) is

the detrained mass per unit area and unit pressure depth, and r(P) is the fraction of the

detrained liquid water which is precipitated. In the NASA GEOS-1 GCM (version 5.2),

r(P) is calculated by

,.(P)= VB(A)t(P.) (19)

Mo(P) involves contributions from all cloud types penetrating level P, i.e.,

_(P) = x(P) _>(P)-_B(A)dA (20)

where A(P) is given by gq. (10), and roB(A) is the cloud-base mass flux per unit k. The

mass detrainment rate

D(P) = _>(P)._B(a) da(P) (21)

dp

In the NASA GEOS-1 GCM, the continuous spectrum of clouds is divided into subensembles

of finite AA and the effects of each subensemble are considered independently. Thus for the

ith spectral band (from Ai - AAi to Ai), we have

_;(P) = r_(_) (22)

0,

and

P > P. (hi),

P. (hi) >_P > P. (hi- _xh_),

otherwise

Di(P) J @' - - (23) 0, otherwise

where PD (hi) is the detrainment level of clouds with h = hi. It may be obtained by solving

Eq. (10).

When hi is small, we may assume that roB(h) --_roB(hi), thus neglecting terms in (Ahi) 2,

Eq. (22) becomes

M_(P) = rlxi(P)mB(hO[h(P)- hi + Ahi],

O,

P > P. (hi),

P. (hi) _>P _>P. (hi - :xh_),

otherwise

(24)

Similarly, approximating the dh/dp in Eq. (23) by Ahi[PD(hi) -pD(hi- Ahi)] -_, we have

Di(P) = { O,rlxi(P)mB('XOA'Xi[PD('XO -- pD(,Xi -- Ahi)]-',

P. (hi) >_P > P. (hi - _xh_),

otherwise

(25)
Then Eqs.

and

where

r_(P) = {

and

rt_(P) = {

(17) and (18) can be rewritten as

(0s) =F_(P)rnB(Ai)AAi _c

(26)

(27)

g,x,(P)[A(P) - A_+ AA_](AA_)-1 o_Op

+g.x, (P) [p. (A_)- p. (A_- AA_)]-I

xlx,(P_)L[1 - ,.(P_)],

0,

P > PD (hi),

p.(A_) > p > p.(A_ - exam),

otherwise

(28)

g_]_,(P)[,_(P)- ,_i + AAi](AAi)-I_

+g._, (P) [p. (A_)- p. (A_- LxA_)]-'

x F_*(P.) - t_(P.)],

O,

P > PD (,_i),

p. (A_)> p > p. (A_- exam),

otherwise

(29)

Finally, the rate of change of the large-scale prognostic variables, the potential temperature

and specific humidity due to the ith cloud subensemble can be written as

oo) _ ._(A_)_x_r_(p ) (a0) -5-{ _ %P

(_-t)c= LrnB(Ai)AAi[Ph(P)- Ps(P)]

and

(31)

d. Mass-flux kernel and cloud-base mass flux

In AS, the change rate of cloud work function is expressed as

27 \--gi-,, _+ \--gi-,, _

where the subscripts c and ls denote the contributions from the cloud-scale and the largescale processes, respectively.

From Eq. (16) and together with some assumptions, i.e., the time variations of pressures

at the cloud base and top equal to zero, ignoring the time dependence of 0 in Eq. (5), and

using

Oh*(P) as(P) (33) O_ _ [I+v(P)] at
the change rate of cloud work function in terms of the rate of change of the large-scale

variables, h and s, can be approximately obtained as

dA_ _ [PB dP

JPD P[1 + v(P)]

× at [1+ v(P)] _,_ + -Ag 0 \ _ [1+ v(P)]--b-_- 34)

In AS, the cloud-scale contributions on dAx/dt is

--J-Jc Jo I(_,_,-_B(A')dA' (35)

where the kernel Kxy represents the rate of change of the cloud work function of cloud

type A per unit cloud-base mass flux of cloud type A_. These are not direct cloud-cloud

interactions, but indirect effects of the various cloud types on each other through their

environment. For quasi equilibrium to hold for the cumulus ensemble as a whole, these

interactions must occur quickly compared to changes in the large-scale forcing. The standard

implementation assumes that they occur instantaneously, resulting in a quasi-static balance

between the cloud ensemble and the large-scale forcing. It is this assumption that results in

an ill-posed problem, with the possibility that either no mass-flux distribution can produce

an exact balance for all clouds with positive buoyancy or that (most frequently) multiple

distributions can satisfy an "overadjustment" problem (Silva-Dias and Schubert 1977) in

which some of the possible cloud types are overstabilized by the effects of other cloud types.

The main assumption of RAS is that the interaction between clouds, represented by the

off-diagonal terms of K in Eq. (35), occurs over a short but finite time and that at

any instant the computations for each cloud and each cloud type are just based on the

"current" environment which already has been affected by influences of some other type

of cumulus convections. In this way, the cloud interactions are taken into account. The

ill-posedness of the original AS implementation is thus removed by solving an initial value

problem that selects an equilibrium distribution that depends on the time scales specified

for the adjustment of the individual cloud types. In RAS, considering the effects of a single

subensemble on the cloud work function, Eq. (35) reduces to

I_,,_,- ._B(A_)AA_ c (36)

Finally, fl'om Eqs. (34), (26) and (27), the approximate I(_,_ is given by

I(_,,_, = D P[1+ v(P)] ×

{F/_(PB)- [1 + "/(P)]Fs(P) + CP _ig JP/'PB o (Ph(P') - [1 q- _/(P)]Ps(P)) dP'} (37)

10
The subensemble cloud-base mass flux fnB (Ai)AAi is obtained by equating the large-scale

and cloud-scale changes of A, i.e., making dAxi/dt = 0, we have

> o, mB(A_)AA_ = _ d_ J_s (38)

O, otherwise.

In the GEOS-1 GCM, the large-scale forcing (dA_dt)ls is calculated directly from time

difference:

(_t_) = A_i(t + At) - A_i(t) (39)

l_ At

where A_i(t 4- At) is the cloud work function calculated from the profiles of 0 and q after

they are modified by the large-scale processes over a time interval At. This is a good

approximation as long as At is small. In the GCM, A_(t + At) is replaced by the value

of the cloud work function after modification of the environment by large-scale effects and

using a cloud-type dependent critical value of the work function instead of Axi(t). Due to

the convective parameterization is designed to nearly neutralize the instability, the critical

value of the work function is near zero.

In RAS, the cloud type is determined by the height of cloud-top level and the kernel is

computed explicitly since only the diagonal elements are required.

In addition to RAS scheme, GEOS-1 GCM employs a Kessler-type scheme for the reevaporation of falling rain (Sud and Molod 1988). The scheme accounts for the rainfall

intensity, the drop size distribution, and the temperature, pressure and relative humidity

of the surrounding air.

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 (about 50 mb thick).

Each time RAS is invoked (every ten simulated minutes), the possibility for shallow convection is checked for the two layers just above the cloud base. RAS also randomly chooses

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

to the model top layer.

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

specific humidity in any grid point 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.

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

the shortwave and longwave radiation packages.

Supersaturation or large-scale cloudiness is defined whenever the large-scale precipitation

scheme determines that the grid box at any level becomes supersaturated. In order to ensure

11
that at any instant the total cloud fraction is less than or equal to one, supersaturation

clouds are only prescribed when there are no deep convective clouds.

Since in the current moist version of GEOS-1 GCM which includes just the moist process,

in the tangent linear model and its adjoint code the cloud formation part is not employed

(i.e., all the related lines are commended out) although this part of code has already been

included.

2.2 Description of the Discretization of the Moisture Physics Parameterization

a. Cloud model

All clouds are assumed to have the same base. We will refer to that cloud type with its

detrainment level in layer i, as the ith cloud type.

The normalized mass flux for each cloud type is a linear function of height,

rib,k-i/2 - rl_,k+l/2 = )_(zk_l/2 - zk+l/2), /;.= i + 1, i + 2, ..., K - 1 (40)

where _li,k+l/2 is the cloud mass flux of the ith cloud type at level k + 1/2 normalized by its

value at the cloud base, Ai is its entrainment rate, and Z1_+1/2 is the height of level k + 1/2.

Equation (40) applies from the layer immediately below the detrainment layer to the layer

immediately above the cloud base which is at K - 1/2. We assume the detrainment occurs

at the middle of the detrainment layer and therefore that there is an additional half-layer

at the top over which the cloud entrains,

,_,_ - ,_,_+1/2 = ),_(z_ - z_+1/2) (41)

The vertical coordinate is specified by the pressure at the half-integer levels (Pk+l/2, /; =

1, 2, ..., K). The discrete form of hydrostatic equation over the full layers is,

Zk-1/2 - Zk+l/2 = % G(G+l/2 - G-1/2),

g

For the hydrostatic equation over the lower half of each layer,

Zk - Zk+l = CPok+l/2(Pk+l -- Pk).

g

0 is the potential temperature, and P is the form as

5+1/2 = (pk+l/2/po) _

and

G-

/_'= 1,2,...,K (42)

/_.= 1, 2,...,K (43)

1 q- t; \ Pk+l/2 Pk-1/2

(44)

(45)

12
which is the form suggested by Phillips (1974) and used in Arakawa and Suarez (1983).

Here _; = t_/cp.

From Eqs.

where

(4o)to (42),wehave

Vi,k-1/2 - Vi,k+i/2 = flkOkAi, k = i+ 1, i+ 2,...,K- 1

flk = _(G+1/2 - G-1/2)

For the last half-layer up to the detrainment level,

_li,i = _li,i+a/2 + fliOiAi

cp

9

(46) and combining with Eq.

_li,i = 1 + Ai

where

Summing Eq.

the detrainment level,

(46)

(47)

(48)

(49)

(48), we obtain the normalized mass flux at

i

Z _0_ (50)

h=K-1

where _1i,I<-1/2 = 1 is used.

The discrete form of the moist static energy budget of each cloud type is

k = i+ 1, i+ 2,...,K- 1 (51)

where hk is the environment moist static energy of layer k, h_,k+l/2 is the cloud moist static

energy of the ith cloud type at level k + 1/2. For the half-layer at the cloud top we have

l]i,ihc, i -- l]i,i+l/2hc, i+l/2 = (l]i,i -- l]i,i+l/2)hi (52)

Summing Eq. (51) and combining Eq. (52), we obtain the cloud-top moist static energy

expression as

i+1

,_,_t¢,_= t_,,_+ E [(,_,J-,,- ,_,j+,,)t_j] + (,_,_- ,_,_+1/_)t_ (sa)

j=K--1

Ignoring precipitation, the cloud total water is

qC qC ,

l]i,h-1/2( i,h-1/2 q- li,h-1/2) -- l]i,h+l/2( i,h+l/2 q- li,h+l/2) = (l]i,h-1/2 - l]i,h+l/2)qh

k.=/+ 1,i+ 2,..., I<- 1 (54)

where q[,k+l/2 and lijc+U2 are the specific humidity and the liquid water mixing ratio of

the ith cloud type at level k d- 1/2. The rhs of Eq. (54) assumes no liquid water in the

13
environment.Assumingqi_,i= q*, where q* is the saturation specific humidity of layer i, we

obtain,

-- q,,+ E - + - q: (551

_]i,i j=K-1

where li,i is the liquid water mixing ratio of the detraining air for the ith cloud type. Since

0 and q are known and the saturation specific humidity can be calculated and li,i can be

calculated if Ai is known.

Assuming at the detrainment layer

h_,i = h_ (56)

where h* is the saturation moist static energy of layer i. In the AS, Eq. (56) is a polynomial

in A_ whose degree depends on the height of the detrainment level. From Eqs. (40), (41),

(53) and (56), the expression for Ai is obtained as

Ai = hi,_- h_ (57) i

b. Cloud work function

The discrete form of the cloud work function is obtained by discretizing Eq. (16) as

i+1

A i z

j=K-1

hc c + _j,_,j_,/_( _,j_,/2- h;)] + _,_,_+,/_(h_,_+,/_- hT)

where Ai is the cloud work function for the ith cloud type and _ and ¢ are defined as

(ss)

and

Using Eq.

environmental quantities as:

Pj+_/2 - Pj

_J- Pj(l+_j) ' j = 1,2, ...,I_ - 1 (59)

Pj - Pj-1/2

_J - Pj(1+ _j) ' j = 1,2, ...,I_ - 1 (60)

(51) to eliminate h _, the cloud work function can finally be written in terms of

A i z

+

_K-1 hK -- _i_]i,i+l/2 h*

h=K-1 j=K-1

i+1

h=K-1

(61)

14
c. Cumulus cfJ)cts on the lar'ge-scale budgets

The cumulus effects on the budgets of dry and moist static energies of the environment are

discretized as

g F 3

Ot )c Apk

(62)

and

-0-[-) c Apk (63)

where (Osk/Ot)_ and (Ohk/Ot)_ are the rate of change of dry and moist static energies of

layer k, Mk+l/2 is the cumulus mass flux at level k + 1/2, Dk is the detrained mass at level

k and Apk = Pk+l/2 - Pk-1/2. The last term on the rhs of Eq. (62) represents cooling from

the re-evaporation of liquid water detrained to the environment. Here lk is the liquid water

mixing ratio, r'k is a cloud-type-dependent precipitation fraction.

In AS implementation, Mk+ll 2 would be the total mass flux of all cloud types penetrating

the level k + 1/2. However in RAS, only one cloud type is considered at a time, then Mk+l/2

has the form

Mh+l/2 = { MB(i)_li,k+l/2 O, , ii <> k, ]% (64)

and

¢ MB(i)_li,i, Dk =

k O,

Also lk = li,i, for i = k, and lk = 0 otherwise. Eqs. (62) and (63) can be rewritten as

(65)

Ot ]o = MB(i)F_(k) (66)

and

where

and

(_) = MB(i)Ft_(k) (67)

c

9 [7]i,k_l12(Sk_l/2 8k ) @ 7]i,k+ll 2

×(_ - _+1/_)- ,_,it_,_L(1-_._)_)],

O,

9 [1]i,k_l/2(hk_l/2 hk) _- 1]i,k+l/2

×(h_- h_+_/_)+ ,_,i(h_- hd_)],

O,

k=i,i+ 1,...,K, (68)

k=1,2,...,i-1

k=i,i+ 1,...,K, (69)

k=1,2,...,i-1

15
where6/k is the Kronecker delta function and it is assumed that _li,i-U2 = _li,K+U2 = O.

Finally the rate of change of potential temperature and specific humidity due to the ith

cloud type can be obtained as

- k= 1,2,. (70)

c CpPk ""

and

c = 7 B(/)[G(k') - G(k')], k.= 1,2,..., I< (71)

d. Mass-flux kernel and cloud-base mass flux

In RAS, for a single cloud type, only the diagonal element Ki,i is required, which is given

by

1 (72)

Ki,i- MB(i) \ dt /_

Using the approximate relation

oh; _ (1+ w) o_ (73)

Ot Ot

andEqs. (57), (61), (66), (67),

[£i,i : (_i,__,+ _)ct_(I0 - (_,_,_+,/2+ _,_,_)(1+ w)c_(/)

+ _ (_k-, + #:_) r:_(IC) + _ (']i,j-,/2 - ']i,j+,/2)r:_(j)

k=K-1 j=K-1

+ 0(_,k__/2 - _<k+_/2)l%_(k)} + 0(_,_ - _<_+_/2)G(i)

i+1

- Z (_'<_+,/2 + _'<_-,/_)(1 + v_)c_(k.),

k=K-I

(74)

where

+

+

(_]i,j-1/2- _]i,j+l/2)hj

_j=K--1

-Gh*_]i,i+l/2

i+1

Z (_-' + _)

k:K-I

i+I

k=K-1

i+1

k:K-I

(75)

16
For simplicity, the rate of change of ¢ and # have been ignored. The rate of change of Ai is

included through the terms involving 0.

With the quasi-equilibrium of cloud work function we obtain,

MB(i) = - k dt /_,

when the rhs of Eq. (76) is positive, otherwise MB(i) = O.

(dAi_ can be calculated by Eq. (39). the ith cloud type \ dt ]l_

cumulus-induced changes in 0 and q can be calculated.

(76)

The large-scale forcing of

Once MB(i) is known, the

In RAS, we assume that all liquid water formed inside a cloud is carried to the top detrainment level (i.e., 1_,_in Eq. (55)), and a fraction of this detrained liquid water is precipitated

and the rest is evaporated within the detrainment layer. Another assumption is the precipitation simply falls to the ground without evaporation. Thus the precipitation Ri for the

ith cloud type can be written as

= (77)

The cloud-type-dependent parameter r'i is set to one for every cloud type in the NASA

GEOS-1 GCM Version 5.2.

3 Tangent Linear Model of the Moist Process Physics Package

The linearized discrete RAS parameterization equations (40) - (77) are derived as follows.

We use {} to describe the basic state trajectory terms and O' to denote the perturbation

variables terms.

For Eq. (44), the corresponding tangent linear formula is

(G+,/_)'= k77_ j (w+,/_)'

{pk+l/2} (78)

For Eq. (45),

(P_)' =

+

-

1+_

_ /9 t {w-,.}(

l+n

- (79)

17
For Eq.

where

For Eqs.

and

For Eq.

(46),

(_,k-1/2)' - (,_,k+_/2)'=

+ {/_}{),d(o_)' + {O_}{),d(/_k)',

k = i+ 1, i+2,...,N- 1

(48) and

= _ [(p_+,,),_ (P__,,)'] (_)' 7

(8o)

(50),

(_]i,i) t =

(81)

+ {/_d{o_}(),d'

+ {;_d{_,d(oo' + {o_}{_,_}(_)', (82)

(9_)' 7

(83)

i

k=K-1

i

+ (a_)' z {_}{o_} (84)

k=K-1

The expression Eq.

cloud type li,i can be linearized as

( _,_)t ' = -(q_)*' + {,_,d-_/ (q_;)'

t.

i+l

j=K--1

+ (1,_,_}- 1,_,_+i/2}) (q_)'

-{q_}(,,_,_+i/2)'}+{,,_,_}-_{{,,_,_+,/_}{q_}-{q,, _}

- i+1 _ ({,_,j-,/2}- {,_,j+,,}){,lj}} (,_,_)'

j=K-1

(55) of liquid water mixing ratio of the detraining air from the ith

(85)

18
Thecloudwork functionexpressionEq.

where

and

(Ad' =

(61) can be linearized as

{_d{hT}0_,_+,/2)'- {,_,_+,/2}{h7}(_)'

j=K--1

j=K--1

, /1

(_j) Iz {Pj+,/2} (pj), [{Pj}(1+ {-_;})]-'(P_+,/2)'- {_}2(1 + {._;})

{_+,/2} - {P_}

{-vT}V ; 7_7__(_)'

(#j)_ =

+

, {r__,/2} (r_)'

- [{r_} (i + {v;})]-' (;-,/2) + {r_}2(i + {v;})

{r__,/2}- {r_}

The expressions of cumulus effects on the large-scale budgets (Eqs.

linearized as:

(o(o_)"_ = (r_(_.))'

{_/_(_)}

(70) and

(86)

(87)

(88)

(71) are

19
and

where

+

{_B(i)}

_p{p_}2{Y_(A.)}(P_)' a.= 1,2,...,1; (s9)

+

{{M,(i)} [(r_(k.))'-(r_(k.))'] L

(:_. (i))'[{r_(k.)}- {r_(k.)}]},

h' = 1_2, .., K

_{,_,_}{_,_}L (1- {_._})a_ /

(Ap_),

+ A--A_k { }k}

O,

k = i,i+ 1,...,K,

k=1,2,...,i-1

(90)

(91)

2O
and

t.

1

(rli,k+i/2)

- {,_,_}((t_t)' - ,n_a¢_,j _')

q

- ' h* - ]

d

O,

(92)

Finally, the cloud work function kernel and cloud-base mass flux expressions t3qs. (74) and

(76), respectively, can be linearized as

/,,_/'-- /_,,,_+_</_,_/,_//'+ (/_,-,/'+/_/')_,_/,,/_

i+1 /

k

+ z [({<,,_}-{<+,,_})/_,_/;//'

j=K--1

k

21
and

where

(0/=

+ {0} ((,_,_-1/2)'-(,_,_+1/2)'){v,,(a.)}

(0/({_,___/_} - {_,_+_/_}){v,_(a.)}} +

+ {0} ({_,_}- {_,_+_/_})(v,_(0)'

+ (0/({_,_} - {_,_+_/_}){v,_(i)}

+ {0} ((_,_)'-(_,_+,/2)') {v,_(O}

- Z _) {_,_+_/_}+ {_}(_,_+_/_)'

_=_,__,

+ {_}(_,___/_)'+ (_)'{_,__1/_}) (1+ {-_}){F_(k.)}

+ ({_}{_,_+_/2}+ {_}{_,___/_})(-_)'{V_(k.)}

+ ({_}{_,_+,/_}+ {_}{_,__,/_}) (_+ {-_})(V_(k.))'_

N

J

(d{A_}'_ {i<,d__(i<,y (d(A,_)"_ {i<,d_1 (MB(i))'=\ gt /_ -\ gt /_

+

+

+

+

+

+

h* i _ i h* • I

-{ _}{ _}(_,_+_/2)- (_) { _}{_,_+_/_}- {_d(h_){_,_+_/_}

i+1 /

•(_,j-i/_)'- (_,_+_/_)'{h _}]

..I

i+1 I( {_}(_,_+,/2)'+ t

(_) {_,_+,/_}

k=K-1

{_}(_,__,/_)'+ (_) '{_,__,/_}) {h_}*

(_ .V.fh*-_

(93)

(94)

22
h=K-1

(95)

3.1 Linearized Discrete Dynamical Equations

3.2 Coding of the Tangent Linear Model

For coding the tangent linear model, we linearize the original nonlinear forward model code

line by line, do loop by do loop and subroutine by subroutine. This amounts to obtain the

exact same tangent linear model as by coding directly from the original linearized model

dynamical equations.

The tangent linear model is the linearized nonlinear forward model in the vicinity of a basic

state which is a model trajectory. Any original code line, we may write it as

V = f(X) (96)

where

X= (z,, z2, "", z_) T (97)

and U is a new derived variable related to the original control variables of the nonlinear

forward model, i.e., it may be one of the original control variables or an intermediate variable

which is a function of the original control variables. Here zl, z2, ..., _'_ (the components

of the vector X) are the required variables to derive U, which may consist of either the

original model control variables or of the intermediate variables derived from the original

control variables; fn is the number of the required variables.

The corresponding tangent linear code assumes the form:

X:Xbasic state X:Xbasic state

_-(_'_ _ X:Xbasic state

of

where X = Xb,sic state means that in the expression _, i = 1, 2, ..., rn, all the values of

the required variables xl, x2, ..., x,_ are chosen to have the exact same values as those of

the basic state trajectory as in the nonlinear forward model to ensure that the basic state

of the integration of tangent linear model is exactly the basic state of the nonlinear model

integrating trajectory. Here 5U and 5xl, 5x2, • •., 5x,_ are the corresponding perturbation

variables of U and xl, x2, ..., x,_, respectively.

In order to obtain the necessary values of Xbasic state, the nonlinear model integrating

trajectory, for the tangent linear model, we must apply the parallel method. This method

consists of calculating in parallel the nonlinear model trajectory as the basic state Xbasic state

and the integration of perturbation variables as well, such as _U, in the tangent linear model.

23
3.3 Notational Convention for Variables and Subroutines Used in the

Tangent Linear Model Code

Forconvenience, thesameoriginalnamesusedin the nonlinearforwardmodelareemployed

forthecorrespondingperturbationvariablesin thetangentlinearmodelcode.Forinstance,

weuse"U" for "5U", "PKHT" for "5(PKHT)", "AKM" for "5(AKM)", etc.. This also

meansthat theperturbationcontrolvariablesin theTLM sharethesamecommonstructure

andsamecommonblocknamesasthe full variablesof the GCM itself. Thus,oneneedsto

payattentionto this issuewhenrunningthe TLM in conjunctionwith the originalGCM.

Weappenda "0" at theendof a variablenamein the originalnonlinearforwardmodelto

representhecorrespondingbasicstatevariable,suchas"U0" for "Ub_sicst_t¢","PKHTO"

for "(PKHT)b_c _" , "AKMO" for "(AKM)b_ _", etc..

For namingsubroutinesin the tangentlinearmodel,wesimply appenda "L" at the beginningof the original namesof subroutinesof the nonlinearforwardmodel. To conform

with ANSI FORTRAN77language,if thenewnameof atangentlinearsubroutineexceeds

six letters,wejust retain its first six letters. For instance,the originalsubroutinesof the

nonlinearmodel "RANG", "MOISTIO" and "RNEVP", havecorrespondingnamesin

thetangentlinearmodelas "LRASG", "LMOIST" and "LRNEVP", respectively.

4 Adjoint Model of the Moist Process Physics Package

4.1 Using the Adjoint Method to Calculate the Gradient of a Cost Function

The practical determination of the adjoint model of the moist physics package used in NASA

GEOS-1 C-Grid GCM is the key computational method enabling us to calculate the gradient

of a cost function with respect to initial conditions (or other control variables) for carrying

out a 4-D variational assimilation. In 4-D variational assimilation, the cost function, which

measures the weighted difference between observations and forecasts in an adequate norm,

is minimized by using a large-scale unconstrained minimization method iteratively which

requires for its implementation the gradient of the cost function with respect to the control

variables. Finally, the optimal state defines a trajectory which passes as close as possible

in a least-squares sense to the observations while satisfying the system of coupled partial

differential equations of the numerical weather prediction model as strong constrains.

Assuming that the cost function consists of a weighted least square fit of the model forecast

24
to the observations, it has the form :

i R T

where X(t,.) is a model state vector of size M(4K _- i) containing the values of the zonal

wind u, the meridional wind v, the potential temperature 0, the surface pressure Ps, and

the surface humidity q; here _/ is the number of grid points at each level; K is number

of vertical levels, t_. is a given time in the assimilation window; x°bs(t_) is a vector of

observations defined over all grid points on all levels at time t_; W(t_) is an N x N diagonal

weighting matrix. From Navon et al. (1992), we have the following expression

R

where X_(t0) is the initial perturbation, X_(t,.) is the perturbation in the forecast resulting

fi'om the initial perturbation, VJ (X(t0)) is the gradient of the cost function with respect

to the initial conditions.

The tangent linear model of the nonlinear forward model can be symbolically expressed as

X'(t,.) = P,.X'(t0) (101)

where P,. represents the result of applying all the operator matrices in the linear model to

obtain X'(t_.) fronl Xt(t0).

We define the adjoint model as

X_'(t0) = P#X(t_.), r' = 1,...,R, (102)

where (^) represents an adjoint variable. After some algebra we obtain (see Navon et al.

1992) that the expression for the gradient of the cost function with respect to the initial

condition is

R 1 N

VJ

From this analysis, we note that the so called adjoint model operator is just the transpose

of the tangent linear model operator.

4.2 Coding of the Adjoint Model

Since the adjoint model equations consist of the transpose of the linearized version of the

nonlinear forward model, if we view the tangent linear model as the result of the multiplication of a number of operator matrices:

P = A1A2...AN, (104)

25
where each matrix Ai(i = 1,...,N) represents either a subroutine or a single DO loop,

then the adjoint model can be viewed as being a product of adjoint subproblems

pT =ANAN_ T T 1 .

..A T. (105)

Thus, the adjoint model is simply the complex conjugate of all the operations in the tangent

linear model. Each DO loop or each subroutine in the tangent linear model has its adjoint

image DO loop and subroutine, respectively. Therefore, we code the adjoint model directly

from the discrete tangent linear model by rewriting the code of the tangent linear model

statement by statement in the opposite direction. This simplifies not only the complexity of

constructing the adjoint model but also avoids the inconsistency generally arising from the

derivation of the adjoint equations in analytic form followed by the discrete approximation

(due to non-commutativity of discretization and adjoint operations).

4.3 Notational Convention for Variables and Subroutines Used in the

Adjoint Model Code

In a similar way as in the tangent linear model, we employed the same original variable

names used in the nonlinear forward model for the corresponding adjoint variables in the

adjoint model code. For instance, we use ;;U" for ;@", ;;PKHT" for ;;(PKI-HT) '', ;;AKM"

for ;;(AKM)", etc.. As in the TLM, this convention also means that the adjoint control

variables in the adjoint model share the same common structure and same common block

names as the GCM itself. Thus, one needs to pay attention to it when running the adjoint

in conjunction with the original GCM.

We also just append a ;;0" at the end of a variable name (in a similar way as done previously in the tangent linear model) to represent the corresponding basic state variable,

such as using ;;U0" for ;;Ub_sic st_t¢", ;;PICHTO" for ;;(PKHT)b_ic _t_t¢", ;;AKMO" for

_(AKM)b_i_ _t_t¢", etc., needed in the adjoint code.

For naming subroutines, we simply change the letter ;;L" at the beginning of the names of

the tangent linear model subroutines to ;;A" and used them as corresponding adjoint model

subroutine names. We also retain the adjoint subroutine names which do not exceed six

letters to conform with ANSI FORTRAN 77 language. For instance, the original subroutines

of the nonlinear model ;;t_ASG", ;;MOISTIO" and ;;t_NEVP", have corresponding names

in the adjoint model as ;;At_ASG", ;;AMOIST" and ;;At_NEVP", respectively.

4.4 Verification of the Correctness of the Adjoint Model

Integrating the nonlinear model forward in time and its adjoint backwards in time, while

forcing the r.h.s.of the adjoint model with difference between model and observations (see

26
Eq. (103)),onecanobtain valuesfor the gradientof costfunctionwith respectto distributed controlvariables,whichmayconsistof eitherthe initial conditionsor the initial

conditionsplusboundaryconditionsormodelparameters.SincethemoistversionofNASA

GEOS-1C-GridGCM consistsof thousandsof linesof code,anyminorcodingerror may

causethefinalgradientof costfunctionwith respecto thecontrolvariablesto beerroneous.

Therefore,weneedto verifythe correctnessof the linearizationandadjointcodingsegment

by segment.Eachsegmentmayconsistof eithera subroutineor of severalDO loops.For

a detailedderivationof the adjoint modeland verificationof its correctness,seeNavonet

al. (1992).

The correctnessof the adjoint of eachoperatorwas checkedby applyingthe following

identity (Navonet al. 1992)

(io6)

where Q represents the input of the original code, A represents either a single DO loop or

a subroutine. The left hand side involves only the tangent linear code, while the right hand

side involves also adjoint code (A'T). If equality (106) holds, the adjoint code is correct

when compared with the TLM. In practice the identity Eq. (106) holds only up to machine

accuracy. In our verifications of the correctness of each segment of the adjoint model and

the whole adjoint model, the LHS and the RHS of Eq. (106) attained 13 digits of accuracy

which is near the machine precision limit. The test were performed at NASA's Cray C-90

computer which has intrinsic double precision. These results show that our adjoint code

consists of absolutely the exact adjoint operators of the TLM of the moist version of NASA

GEOS-1 C-Grid GCM.

A gradient check (Figs. 1, 2) was then performed to assess accuracy of the discrete adjoint

model. This verification method is described next. First, we chose the cost function J as in

Eq. (99) and the N x N diagonal matrix W = diag(W_,W_,Wo, Wq,Wps), where the

diagonal submatrices are defined as : W_ = 5 × 10-1I s2m -2, W_ = 5 × 10-1I s2m -2,

W0 = 10-3I K -2, Wq = 5 × 10-3I, Wps= 10-3I rob-2. Then, let

J(X + (_h) = J(X) + (_hTVJ(X) + O((_2), (lO7)

be a Taylor expansion of the cost function. Here _ is a small scalar and h is a vector of

unit length (such as h = VJ/IIVJII). Rewriting the formula above we can define a function

of _ as

(P(o_) = J(X + o_h)- J(X)

(_hTVJ(X) = 1 + O((_). (108)

27
7

6

z

c_

4

3

2

/

8

8

PHI WITH RESPECT TO LOG18{£LPH£}

I I I I I I I

_Y

i i i i i i i

i 17 16 15 14 13 12 I/

x

Figure 1: Variation of the @(o 0 with respect to log _ (gradient check of correctness of adjoint

model). Integration period is 6 hours and t -- 6 hours model generated observations were

used. January 1, 1985 00Z DAO_s data was used as t = 0 observations. The first guess is

the shifted 6-hour initial condition. The time integration scheme employed is the leapfrog

scheme.

PHI WITH RESPECT TO LOG18{£LPHA}

2 I I I I I I I

/

o

/

2

z

<

s

i i i i i i

7 8 i 17 16 15 14 13 12 I/

x

Figure 2: As in Figure 1, but for the Matsuno time integration scheme.

28
.5

' 1 @

i

c_

<

1.5

I

G 2.@

o

s

2.5

3.@

3.5

4.@ i

8

LOG/@IPHI[ALPHA) /I TO LOG/@(ALPHA)

I I I I I

R

R

£

/

/

/

£

£

i i i i i i i

17 16 15 14 13 12 I/

x

Figure 3: Variation of the log I_(a) - 11 with respect to log a. Integration period is 6 hours

and t = 6 hours model generated observations were used. January 1, 1985 00Z DAO% data

is used as t = 0 observations. The first guess is the shifted 6-hour initial condition. The

time integration scheme employed is the leapfrog scheme.

LOG/@IPHI[ALPHA) /I TO LOG/@(ALPHA)

I I I I I

Figure 4: As in Figure 3, but for the Matsuno time integration scheme.

29
Forvaluesof _ whicharesmallbut not toocloseto the machinezero,oneshouldexpectto

obtainvaluesfor _(_) whicharecloseto unity. Weobtainedvaluesfor thefunction _(_)

equalunity to a highdegreeof accuracywhenthe parameter_ is variedfrom 10-2 to 10-6.

from the residualof _(_) (figs. 3, 4), we found that the residual tends to zero. The

gradient check verifies that the adjoint model is correct and can be used, for example, to

perform 4-D VDA experiments.

Although the gradient check results are good, comparing the gradient check results with

those for the adiabatic version of the NASA GEOS-1 GCM (Yang and Navon 1996), we

find that including moist processes into the GCM decrease the validity of the tangent linear

approximation. It worsens both in terms of the accuracy of the gradient values of the cost

functional as well as in the range of the perturbations, in which satisfactory values of the

gradient of the cost functional are maintained. The reasons are:

1) The high nonlinearity of the original moist process package including the RAS cumulus

parameterization and large-scale precipitation and evaporation scheme. Since the ArakawaSchubert parameterization is the most complex cumulus convective parameterization scheme

which provides the most complete physics approach and has an inherent ]terative feature,

the nonlinearity of the AS scheme as well as the RAS scheme is much stronger than that

of some other moist parameterization schemes, such as Kuo's scheme.

2) In AS and RAS, the on-off discontinuous effects are more pronounced than in other

types of moist parameterization schemes. This is due to the fact that there are more on-off

switch processes used in Arakawa-Schubert type scheme and the ]terative feature of the

AS parameterization. Obviously these discontinuities worsen the validity of the tangent

linear approximation and cause the values of the calculated gradient of the cost functional

to exhibit jumps.

Several research efforts were presented related to the serious influence of the on-off switch

processes on the validity of the tangent linear approximation. For instance, Vukicevic

and Errico (1993) tested the accuracy of the tangent linear model of a mesoscale model,

compared the "true" perturbation obtained by direct nonlinear integration and concluded

that significant errors may be expected in the regions where the moist diabatic processes

are important for finite perturbations in the initial conditions. Recently, both Bao and

Kuo (1995) and Xu (1996) carried out a detailed study using idealized continuous examples

with delta function mimicking the on-off switches in physical schemes. They indicated that

ignoring the variation of the switch point due to the perturbation in the initial conditions,

i.e., keeping the switching point in the tangent linear model the same as in the basic state,

could cause significant errors in the tangent linear model solution and gradient calculation.

Thus how to deal with the on-off switches used in the moist process parameterization is a

crucial issue.

The essence of the influence of on/off discontinuous processes is that they may cause sudden

jumps in the model integration trajectories, these jumps changing the model trajectories

30
and forcing the value of the cost function and its gradient to undergo a sudden change.

This sudden change is equivalent to introducing a high nonlinearity into the variational

assimilation system. As a consequence, it may cause either a failure or a slow-down of the

minimization processes in 4-D Vat. Simply smoothing the parameterization at the discontinuous points such as tested by (Zupanski 1993; Tsuyuki 1996a, b, c) cannot completely

remove the influence of on/off discontinuous. Also the simple smoothing method may introduce a negative effect: that of changing the character of the original parameterization

system.

3) Truncation errors may worsen the validity of the tangent linear approximation. In the

original code of the RAS scheme, there are some computational steps involving processes

whereby a very small output results from the difference between two very large terms whose

values are of similar magnitude, such as in calculating the rate of change of large-scale

variables by the forcing of cumulus-scale processes. These very small time tendency terms

result from the difference between two corresponding variable terms. Besides, a nonlinear

term involving N dependent variables in the nonlinear forward model will create N terms in

tangent linear computations. These newly created computation processes may increase the

truncation error. We have checked the validity of the tangent linear approximation term

by term and line by line in the tangent linear model, and found that the truncation error

is a contributing factor towards a reduction in the range of validity of the tangent linear

approximation.

Acknowledgments

This work has been supported by NASA-Grant NAG 5-1660 managed by Dr. Ken Bergman,

Section Head, Climate Modeling.

31

References

Arakawa, A., 1969: Parameterization of cumulus convection, Proc. WMO/IUGG Syrup.

Numerical Weather Prediction, Tokyo, 20 November-4 December, 1968, Japan Meteorology Agency, IV, 8, 1-6.

Arakawa, A., 1971: A parameterization of cumulus convection and its application to numerical simulation of the tropical general circulation Paper presented at the 7th Tech.

Conf. on Hurricanes and Tropical Meteorology, Barbados, Amer. Meteor. Soc.

Arakawa, A., 1972: Parameterization of cumulus convection, Design of the UCLA general

circulation model, Numerical simulation of weather and climate, Tech. Rept. 7, Dept.

of Meteorology, University of California, Los Angeles.

Arakawa, A., and M. J. Suarez, 1983: Vertical Differencing of the Primitive Equations in

sigma coordinates, Mort. Wea. Rev., 111, 34-45.

Arakawa, A., and W. H. Schubert, 1974: Interaction of a cumulus cloud ensemble with the

large-scale environment, Part I, d. Atrnos. Sci., 31,674-701.

Bao, J.-W., and Y.-H. Kuo, 1995: On-off switches in the adjoint method: step functions,

Mort. Wea. Rev., 123, 1589-1594.

Betts, A. K., 1973a: Nonprecipitating cumulus convection and its parameterization, Quart.

d. Roy. Meteor'. Soc., 99, 178-196.

Betts, A. K., 1973b: A composite mesoscale cumulonimbus budget, d. Atmos. Sci., 30,

597-610.

Caselski, B. F., 1974: Cumulus convection in weak and strong tropical disturbances, d.

Atmos. Sci., 31, 1241-1255.

Charney, J. G., and A. Eliassen, 1964: On the growth of the hurricane depression, d. Atmos.

Sci., 21, 68-75.

Fraedrich, K., 1973: On the parameterization of cumulus convection by lateral mixing and

compensating subsidence, Part 1, d. Atmos. Sci., 30,408-413.

Gray, W. M., 1972: Cumulus convection and large-scale circulations, Part III. Broad scale

and meso scale considerations, Atmos. Sci. Paper" No. 190, Colorado State University,

80 pp.

33
Hack,J. J.,W. H.Schubert,andP.L.SilvaDias,1984:A spectralcumulusparameterization

for usein numericalmodelsof thetropicalatmosphere,Mo_t.Wea. Rev., 112,704-716.

Kao, C.-Y., and Y. Ogura, 1987: Response of cumulus clouds to large-scale forcing using

the Arakawa-Schubert cumulus parameterization, J. Atmos. Sci., 44, 2437-2458.

Koster, R. D., and M. J. Suarez, 1992: Modeling the land surface boundary in climate

models as a composite of independent vegetation stands, d. Geophys. Res., 97, 2697-

2715.

Koster, R. D., and M. J. Suarez, 1996: Energy and water balance calculations in the

Mosaic LSM, NASA Tectmical Meruora_duru 104606, Vol. 9, Technical Report Series

on Global Modeling and Data Assimilation, Max J. Suarez, Editor.

Kuo, H. L., 1965: On the formation and intensification of tropical cyclones through latent

heat release by cumulus convection, J. Atmos. Sci., 22, 40-63.

Kuo, H. L., 1974: Further studies of the parameterization of the influence of cumulus

convection on large-scale flow, d. Atmos. Sci., 31, 1232-1240.

Krishnamurti, T. N., 1969: An experiment in numerical prediction in equatorial latitudes,

Quart. J. Roy. Meteor'. Soc., 95, 594-620.

Krishnamurti, T. N., H. Pan, C. B. Chang, J. Ploshay, and W. Oodaly, 1979: Numerical

weather prediction for GATE, Quart. J. Roy. Meteor'. Soc., 105, 979-1010.

Krishnamurti, T. N., Y. Ramanathan, H.-L. Pan, R. J. Pasch, and J. Molinari, 1980:

Cumulus parameterization and rainfall rates I, Mor_. Wea. Rev., 108,465-472.

Lilly, D. K., 1960: On the theory of disturbances in a conditionally unstable atmosphere,

Mor_. Wea. Rea., 88, 1-17.

LSpez, R. E., 1972a: Cumulus convection and large scale circulations, Part I. A parameteric

model of cumulus convection, Atmos. Sci. Paper" No. 188, Colorado State University,

100 pp.

LSpez, R. E., 1972b: Cumulus convection and large scale circulations, Part II. Cumulus

and meso scale considerations, Atmos. Sci. Paper" No. 189, Colorado State University,

63 pp.

Lord, S. J., and A. Arakawa, 1980: Interaction of a cumulus cloud ensemble with the

large-scale environment, Part II, d. Atmos. Sci., 37, 2677-2692.

Lord, S. J., 1982: Interaction of a cumulus cloud ensemble with the large-scale environment,

Part III: Semi-prognostic test of the Arakawa-Schubert Cumulus Parameterization, J.

34
Atmos. Sci., 39, 88-103.

Lord, S. J., W. C. Chao, and A. Arakawa, 1982: Interaction of a cumulus cloud ensemble

with the large-scale environment, Part IV: The discrete model, d. Atmos. Sci., 30,

104-113.

Manabe, S., J. Smagorinsky, and R. F. Strickler, 1965: Simulated climatology of a general

circulation model with a hydrologic cycle, Mor_. Wea. Rev., 93, 769-798.

Moorthi, S., and A. Arakawa, 1985: Baroclinic instability with cumulus heating, d. Atmos.

Sci., 42, 2007-2031.

Moorthi, S., and M. J. Suarez, 1992: Relaxed Arakawa-Schubert: A parameterization of

moist convection for general circulation models, Mor_. Wea. Rev., 120, 978-1002.

Navon, I. M., X. Zou, J. Derber, and J. Sela, 1992: Variational data assimilation with an

adiabatic version of the NMC spectral model. Mor_. Wea. Rev., 120, 1433-1446.

Ooyama, K., 1964: A dynamical model for the study of tropical cyclone development,

Geofis. h_st., 4, 187-198.

Ooyama, K., 1969: Numerical simulation of the life cycle of tropical cyclones, J. Atmos.

Sci., 26, 3-40.

Ooyama, K., 1971: A theory on parameterization of cumulus convection, J. Meteor. Soc.

Japar_, 49, 744-756.

Phillips, N. A., 1974: Application of Arakawa's energy conserving layer model to operational numerical weather prediction, Office Note 104, Natior_al Meteorological Cer_ter',

NWS/NOAA, 40 pp.

Ramanathan, Y., 1980: Cumulus Parameterization in a case study of a monsoon depression,

Mor_. Wea. Rev., 108, 313-321.

Vukidevid, T., and R. M. Errico, 1993: Linearization and adjoint of parameterized moist

diabatic process, Tellus, 45A, 493-510.

Schubert, S. D., J. Pfaendtner, and R. Rood, 1993: An assimilated data set for Earth

Science applications, Bull. Ameri. Meteor'. Soc., 74, 2331-2342.

Silva-Dias, P. L., and W. H. Schubert, 1977: Experiments with a spectral cumulus parameterization theory, Atmos. Sci. Pap., No. 275, Colorado State University, 132

pp..

Suarez, M. J., and L. L. Takacs, 1995: Documentation of the ARIES/GEOS Dynamical

35
Core,"Version2, NASA Technical Memorandum 104606, Vol. 5, Technical Report

Series on Global Modeling and Data Assimilation, Max J. Suarez, Editor, April 1995.

Sud, Y., and A. Molod, 1988: The roles of dry convection, cloud-radiation feedback processes

and the influence of recent improvements in the parameterization of convection in the

GLA GCM, Mort. Wea. Rev., 116, 2366-2387.

Sud, Y. C., W. C. Chao, and G. K. Walker, 1991: Contributions to the implementation

of the Arakawa-Schubert cumulus parameterization in the GLA GCM, J. Atmos. Sci.,

48, 1573-1586.

Takacs, L. L., A. Molod, and T. Wang, 1994: Goddard Earth Observing System (GEOS)

General Circulation Model (GCM), "Version 1, NASA Technical Memorandum 104606,

Vol. 1, Technical Report Series on Global Modeling and Data Assimilation, Max J.

Suarez, Editor, September 1994.

Tiedtke, M., 1989: A comprehensive mass flux scheme for cumulus parameterization in

large-scale models, Mort. Wea. Rev., 117, 1779-1800.

Tsuyuki, T., 1996a: Variational data assimilation in the tropics using precipitation data,

part I: Column model, Meteorology and Atmospheric Physics, 60, 87-104.

Tsuyuki, T., 1996b: Variational data assimilation in the tropics using precipitation data,

part II: 3-D model, Mort. Wea. Rev., 124, 2545-2561.

Tsuyuki, T., 1996c: Variational data assimilation in the tropics using precipitation data,

part III: Assimilation of SSM/I precipitation rates, Submitted to Mort. Wea. Rev..

Xu, Q., 1996: Generalized adjoint for physical processes with parameterized discontinuities,

Part I: Basic issues and heuristic examples, J. Atmos. Sci., 53, 1123-1142.

Yanai, M., 1964: Formation of tropical cyclones, Rev. Geophys., 2, 367-414.

Yanai, M., 1971a: A review of recent studies of tropical meteorology relevant to the planning

of GATE, Experimental Design Proposal by the Interim Scientific and Management

Group (ISMG), WMO, Vol.2, Annex I, 1-43.

Yanai, M., 1971b: The mass, heat and moisture budgets and the convective heating within

tropical cloud clusters, Paper presented at the 7th Tech. Cortf. on Harricanes and

Tropical Meteorology, Barbados, Amer. Meteor. Soc..

Yanai, M., S. Esbensen, and J.-H. Chu, 1973: Determination of bulk properties of tropical

cloud clusters from large scale heat and moisture budgets, J. Atmos. Sci., 30,611-627.

Yang, W., and I. M. Navon, 1995: Documentation of the tangent linear model and its

36
adjointof the adiabaticversionof the NASAGEOS-1C-gridGCM- version5.2,Tcchr_ical Report, FSU-SCRI-95T-52, Supercomputer Computations Research Institute,

Florida State University, 60 pp.

Yang, W., and I. M. Navon, 1995: Documentation of the tangent linear model and its

adjoint of the adiabatic version of the NASA GEOS-1 C-grid GCM - version 5.2, NASA

Techr_ical Memorar_dum, No. 104606, Vol. 8, Technical Report Series on Global

Modeling and Data Assimilation, Max J. Suarez, Editor, March 1996; Also available at

http://dao.gsfc.nasa.gov/subpages/tech-reports.html, Volume 8.

Zupanski, D., 1993: The effects of discontinuities in the Betts-Miller cumulus convection

scheme on four-dimensional variational data assimilation, Tellus, 45A, 511-524.

37
Form Approved

REPORT DOCUMENTATION PAGE 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

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

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

May 1997 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

Technical Report Series on Global Modeling and Data Assimilation Volume 11 -

Documentation of the Tangent Linear and Adjoint Models of the Relaxed Arakawa- C NAS5-32332

Schubert Moisture Parameterization of the NASA GEOS-1 GeM (V5.2) Code 910.3/913

6. AUTHOR(S)

Weiyu Yang, I. Michael Navon, and Richardo Todling

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

Climate and Radiation Branch

Data Assimilation Office

Goddard Space Flight Center

Greenbelt, Maryland

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

National Aeronautics and Space Administration

Washington, DC 20546-0001

8. PEFORMING ORGANIZATION

REPORT NUMBER

97B00053

10. SPONSORING / MONITORING

AGENCY REPORT NUMBER

TM-104606, Vol. 11

11. SUPPLEMENTARY NOTES

Yang: Supercomputer Computations Research Institute, Florida State University, Tallahassee, Florida;

Navon: Dept. of Mathematics and SCRI, Florida State University, Tallahassee, Florida;

Todling: General Sciences Corporation, Laurel, Maryland

12a. DISTRIBUTION / AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 46

Availability: NASA CASI (301) 621-0390.

12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

A detailed description of the development of the tangent linear model (TLM) and its adjoint model of the Relaxed

Arakawa-Schubert moisture parameterization package used in the NASA GEOS-1 C-Grid GCM (Version 5.2) is

presented. The notational conventions used in the TLM and its adjoint codes are decribed in detail.

14. SUBJECT TERMS

Tangent linear model, TLM, Relaxed Arakawa-Schubert moisture parameterization package,

GEOS-1 C-Grid GCM (version 5.2), Notational conventions.

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION

OF REPORT OF THIS PAGE

Unclassified Unclassified

NSN 7540-01-280-5500

19. SECURITY CLASSIFICATION

OF ABSTRACT

Unclassified

15. NUMBER OF PAGES

56

16. PRICE CODE

20. LIMITATION OF ABSTRACT

UL

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39.18

298-102