NASA Technical Memorandum 104606, Vol. 8 /J

Technical Report Series on

Global Modeling and Data Assimilation

Max J. Suarez, Editor

Volume 8

Documentation of the Tangent Linear Model

and Its Adjoint of the Adiabatic Version of the

NASA GEOS-1 C-Grid GCM - Version 5.2

Weiyu Yang and I. Michael Navon

March 1996

DATA ASSIMILATION OFFICE

NASA Technical Memorandum 104606, Vol. 8

Technical Report Series on

Global Modeling and Data Assimilation

Max J. Suarez, Editor

Goddard Space Flight Center

Greenbelt, Maryland

Volume 8

Documentation of the Tangent Linear Model

and Its Adjoint of the Adiabatic Version of the

NASA GEOS-1 C-Grid GCM - Version 5.2

Weiyu Yang

I. Michael Navon

Florida State University

Tallahassee, Florida

National Aeronautics and

Space Administration

Goddard Space Flight Center

Greenbelt, Maryland

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

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

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

its adjoint model of the adiabatic version of NASA GEOS-1 C-Grid GCM(Version 5.2)

is presented. The derivations of and the methods for coding the TLM and its adjoint as

well as tile notation conventions used in these two models are described in detail. The

flow charts of the NASA GEOS-1 GCM, its tangent linear model and adjoint model

are provided. The procedures and their results of correctness verification of the TLM

and the adjoint model are presented. Finally, tutorial examples of derivation of adjoint

code from the tangent linear code are provided for the benefit of various users.

111

Contents

List of Figures iv

List of Tables

1 Introduction

2 Description of the NASA GEOS-1 GCM

4

2.1 Basic original atmospheric dynamical equations ................

2.2 Description of discretization methods and structure of the model ......

2.3 Discrete dynamical equations ..........................

2.4 Flow chart of the NASA GEOS-1 GCM ....................

Tangent linear model of the adiabatic version of NASA GEOS-1 C-Grid

GCM

3.1 Linearized discrete dynamical equations ....................

3.2 Coding the tangent linear model ........................

3.3 Notational convention for variables and subroutines used in the tangent linear

model code ....................................

3.4 Flow chart of the dynamical core of the tangent linear model ........

3.5 Verifying correctness of the tangent linear model ...............

Adjoint model of the adiabatic version of the NASA GEOS-1 C-Grid

GCM

4.1 Using adjoint method to calculate the gradient of the cost function .....

4.2 Coding the adjoint model ............................

V

2

2

3

4

5

12

12

14

15

16

19

23

23

24
4.3 Notationalconventionfor variablesandsubroutinesusedin theadjointmodel

code ........................................ 24

4.4 Flowchartof the dynamicalcoreof the adjointmodel ............ 25

4.5 Verificationof thecorrectnessof the adjointmodel .............. 28

4.6 Someexamplesof codingthe discreteadjointmodelfrom the tangentlinear

modelcode .................................... 31

References 59

vi
List of Figures

Correlation coefficient between the N fields and the L fields with respect to

integration periods .................................

Relative error IIDII/IILN with respect to integration periods .........

Variation of the ¢(¢_) with respect to log a ...................

Variation of the log I¢(a)- 11 with respect to logcr ...............

22

22

30

30

vii
List of Tables

Correlation Coefficients Between N Field and L Field ............ 20

Relative Error _ (%) ............... 20

IILI 1 ...............

Relative Error _ (%) of NASA/GLA 3-D SLSI GCM ........... 21

.°.

VII!
1 Introduction

The GEOS-1 C-Grid GCM was developed by Data Assimilation Office (DAO) at Goddard

Laboratory for Atmosphere (GLA), NASA 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 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 NASA GEOS-1 C-grid GCM 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 GEOS-1 C-grid GCM can be used independently

in another GCM. For instance, full physics package of GEOS-1 C-grid GCM has been used

into NASA/GLA Semi-Lagrangian Semi-Implicit (SLSI) GCM. Having developed the tangent linear model(TLM) and adjoint model of the NASA GEOS-1 C-grid GCM contributes

to various applications involving other 4-D variational data assimilation systems with different GCMs.

Tlle earliest predecessor of the GEOS-1 C-Grid GCM was developed in 1989 based on

"plug-compatible" 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 physical

parameterizations together with plug-compatible concept of "Dynamical Core" introduced

by Suarez and Takacs (1994) facilitated development and testing of new algorithms. Together DAO and Climate and Radiation Branch at GLA, NASA have produced a library

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

GCM applications.

In order to obtain a 4-D variational data assimilation system based on the NASA GEOS-1

C-grid GCM, a first prerequisite is to develop the tangent linear and its adjoint model, two

key parts of any four dimensional variational assimilation system. With the tangent linear

and its a djoint model, the 4-D variational data assimilation (VDA) system of the adiabatic

version of the NASA GEOS-1 GCM was employed successfully to carry out a series of 4-D

VDA experiments to research the Hessian precondition methods and to test a new proposed

Hessial_ estimation algorithm (Yang, et al, 1995). This document describes the development

of the tangent linear model and its adjoint model of the adiabatic version of NASA GEOS-1

C,-(b'id (',(_M.

h, So('tion 2 w_' provide a condensed description of NASA GEOS-1 GCM, which includes the

basic ()rigin_l almosl)h('ri(' dynamical equations and their discrete forms, the discrelization

_n(,l h_,(ls as well as lh,, m(,dol structure. Then we present a detailed flow chart of its code,

which sh()uld pr,)ve*()I,_' usol'l_l lo tirst time users. Section 3 describes and do('um(,uIs in
detail the derivationof and codingthe tangentlinearmodel(TLM) of the GEOS-1GCM,

its flowchart andits correctnessverificationagainstthefull nonlinearforwardcode,aswell

asthe notation conventionsusedin the TLM. Section4 describesin detail the derivation

of the adjoint modelcode(for theadiabaticversionof the GEOS-1GCM). A flowchartof

the adjoint codeis providedalongwith notationconventionsandadjointmodelcorrectness

verification procedures. Finally, tutorial examples of derivation of adjoint code from the

tangent linear code are provided for the benefit of various users.

2 Description of the NASA GEOS-1 GCM

2.1 Basic original atmospheric dynamical equations

In NASA GEOS-1 C-Grid GCM, a a vertical coordinate is defined by

P - PT

(7--

7(

(1)

where 7r - Ps - PT, Ps is the surface pressure and PT is a constant prescribed pressure at

the top of the model atmosphere. In the current NASA GEOS-1 GCM version, pT=O.

The basic original atmospheric dynamics equations of the NASA GEOS-1 GCM are as

follows(for obtaining the adiabatic version, just need to delete the terms related to diabatic

processes). The continuity equation is

8_ 8(_&)

8--/= -Vs-(_V) 8_ (2)

where V is the horizontal velocity vector. The state equation is

= -- (0,) (3)

where o_is the specific density, 0 - T/P is the potential temperature, T is the temperature,

P - (p/po) '_, _; = R/cp, R is the gas constant, cp is the specific heat at constant pressure,

and P0 is a reference pressure which be taken as P0 = 1000 hPa.

The hydrostatic equation is

where (I) is the geopotential.

8¢

3--f = -CpO (4)

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

conserving differencing scheme:

0(_e) 8(_&0) _Q

0t - -vs-(.vo) 8_ + cp--_ (5)

2
whereQ is the diabaticheatingperunit mass.

The equationsof tendenciesof an arbitrary numberof atmosphericconstituents,suchas

water vaporandozone,arealsowritten in flux form:

Ot - Oa + (6)

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

of air.

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

Arakawa and Lamb (1981), to facilitate derivation of an energy- and enstrophy-conserving

differencing scheme.

OV .0V _ V_((I) + K) -%OVoP g OT

--ot = -(f + _')k x V - a--_ - --- 7r Oa (7)

where f is the Coriolis parameter, k is the unit vector in the vertical, _"= V_ x V is the

1

vertical component of the vorticity along a surfaces, h" = 7(V. V) is the kinetic energy per

unit mass, g is the acceleration of gravity, and T is the horizontal frictional stress.

2.2 Description of discretization methods and structure of the model

In GEOS-1 GCM, a Lorenz grid is used in the vertical, with both winds and temperatures

defined at the same levels. The atmosphere between a = 0 and _r = 1 is divided into LM

layers. At these LM layers, the velocity, the potential temperature and the specific masses

of all trace constituents are defined. The vertical velocity dr is defined at interfaces between

layers and a.t top and bottom surfaces.

On horizontal grids, the prognostic variables are located on an Arakawa C grid. The

temperature, pressure and all tracers are located at "p-points", which exclude the poles.

The "u-points", at which zonal wind components are defined, are located between "p-points"

and on the same latitude circles, while "v-points" are located between "p-points" and on

the same meridians. The vorticity is defined at "_-points" on the same latitude circles as v

and on the same meridians as u.

The discretization of the momentum equation is carried out with a second-order energy'

and t)otential enstrophy conserving scheme of Sadourny described by Burridge and Haseler

(1977). A simple second-order finite difference scheme is used for discretizing the thermodynamic equation and the continuity equation.

A polar Fourier filter is applied to the tendencies of all the prognostic variables. The

purpose of the polar filter is to avoid linear computational instability due to convergence
of the meridiansnearthe poles. The filter acts polewardof about 45° latitude, and its

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

wavenumbers and the amount by which they are damped.

The time differencing scheme used in NASA GEOS-1 C-Grid GCM is the Brown-Campana

scheme (1978). It is an explicit scheme used in conjunction with a leap-frog differencing

scheme that relaxes somewhat the instability condition for gravity waves. The basic idea of

the Brown-Campana scheme is to average the pressure gradient force over three time levels.

In NASA GEOS-1 C-Grid GCM, a simple strategy is used which assumes the pressure

gradient force is linearized. That is, to average the three time levels of the mass field and

only to compute the pressure gradient once, while the averaged mass field is used only for

the pressure gradient calculations. Asselin (1972) time filter and Shapiro filter are also used

in the dynamic core of NASA GEOS-1 C-Grid GCM.

The NASA GEOS-1 C-Grid GCM has a resolution of 5° x 4 ° longitude-latitude grid points

in horizontal plane and 20 a-coordinate levels. The time step used is five minutes.

2.3 Discrete dynamical equations

The discrete dynamical equations used in the adiabatic version of the NASA GEOS-1 C-Grid

GCM are as follows (for all the definition of the used symbols and the detailed method of

the Arakawa-Lamb C-grid discrete scheme please see Documentation of the ARIES/GEOS

Dynamical Core (Suarez and Takacs, 1994)). The equations are presented here for ease of

reference to corresponding equations in the tangent linear model.

The hydrostatic equations are

(_LM = (_s + CpOLM( [:)LM +I -- PLM )

and

^

¢1 = (I)l+l + epOl+l(Pl+l -- Pt),

The continuity equations are

O_ri,j _ LM 1

o--7-- -

1=1

and

07ri,j l 1

(7r_)i'J'/+l --_ --O'l+l 0t /_1 (A )i,j

The thermodynamic equation is

O( rOt) ,j 1

(s)

fort= 1, LM-1 (9)

--[6iu? + _jv/]i,.i(fia)t for l = 1, LM - 1,

(10)

(11)

(12)
The component forms of the momentum equation are

OUi,j,l 1

-- [Oli,jVi+l,j+l "t- ]_i,j i,j+l "_- "[i,jVi,j -1- _i,jVi+l,j +

0t (A_x)i,j

+u_;ju_,j_1 - ui,j+iu,;j+l - e,+_,ju_*+i,j+ e_,ju_--i,j]t -

1 1

- u,_,)+ - u,)],,j (_)_,j(_t_)_ 2

1 _i(¢z + I(t) + cpOt _ l i,j (A_z)i,j _iTr

Ovi,j,t 1

-- [Oli_l,j_l tt__l,j_ 1 + _i,j-lU_,j-1 "-_ ot (zXvy)_

-_-"[i,jUi:j nI- _i-l,jU__l, j Jc

+_,jv_*+,,j - .__,,jvi'__,j - _,jv,;j+_ + _,j-lV_*j_l]t1 1

[(.),(v, 7-_.,J .-j

-- Vl_l) + (TrO')l+l(Vl+l -- ?)l)]{, j

(_J)_,j(_to)t 2

1 _j(_l + Kl) + cpOl

for j= 1,JM-1

forj=2, JM- 1

(13)

(14)

Here the c_, _, _,/f, e, qo, u and # are linear combinations of neighboring potential vorticities.

For further details concerning this GCM, we refer to Suarez et al., (1994) and Takacs et

al. (1994). The detailed flow chart of the NASA GEOS-1 GCM with full physics package

is in the next subsection (for getting the adiabatic version of it_ one just needs to skip the

computational processes related to the physics packages and the moisture processes).

2.4 Flow chart of the NASA GEOS-1 GCM

At the next six pages we will present the flow chart of the NASA GEOS-1 C-Grid GCM.
START >

CALL SETUP:Initial model set up, set control parameters[

for model execution and output. J

CALL DIAGSIZE: Diagnostic memory allocation. ]

CALL RESTART: GET initial conditions. ]

CALL CLRANAL: Initialize analysis tendencies to zero. ]

CALL DIAGDRVN: Set pointer locations for diagnostics] turned on in NAMELIST output. J

GEOS-11ARIES GCM FLOW CHART, page 1.

Compute number of distinct frequencies]

Call output routines at begining of experiment]

_" Start the main loop to execute the GCM simulatioo_

CALL SETMET: Initialize resolution dependent] terms in model.

CALL GETBCS: Update GCM boundary conditions]

physics flags. /

Set

GEOS-11ARIE$ QCM FLOW CHART, page 2.

6
CALL CTOA: Convert gridded data "C'" forgridded physics datapackages to "A" us /o}

Compute the moist processes, ]

Relaxed Arakawa-Schubert scheme]

and large-scale convection. J

CALL MOISTIO:

CALL SWRIO: Compute the short wave

radiation processes.

CALL LWNIO: Compute the long wave

radiation processes.

CALL TURBIO: Compute parameterizationthe turbulence processes.

GEOS-1/ARIES GCM FLOW CHART, page 3,

CALL ATOC: Convert "A" gridded data to "'C'

gridded data to execute the hydro

dynamical processes.

Calculate the diagnostics, including

1. total diabatic U-tendency;

2. total diabatic V-tendency;

3. total diabatic T-tendency;

4. total diabatic q-tendency;

5. the analysis tendencies increment of U, V, T, q;

6. incident solar radiation;

7. net solar radiation at the ground;

8. solar radiation heating.

_e tStarting the hydrodynamical processes,

ime integration scheme is either the Matsuno .,fl

wo-step scheme or the leapfrog schem_..11"

GEOS-1/ARIES GCM FLOW CHART, page 4.
_p of the

Matsuno scheme,

the predictor. The leapfrog time

integration scheme.

The second step, the corrector.

I

.............. __,y_o_,____ _.y____ __o_./_ ________,?___ ............

[ CALL TICK: Update time information. ]

GEOS-IIARIES GCM FLOW CHART, paga 5,

Compute more diagnostics, including

1. averaged P-field;

2. averaged U-field;

3. averaged V-field;

4. averaged T-field;

5. averaged q-field;

6. averaged QQ-field;

7. precipitable water;

8. temperature and moisture

convergence diagnostics.

Check for pressure diagnostic, ensure

that Ps < Psmax and Ps > Psmin.

[ Update alarm flags. /

GEO_I/ARIE8 GCM FLOW CHART, page 6.
Call all output routines to output prognostic |

1

and diagnostic results. J

CALL RESTART: write out restart and

diagnostic results to

disk.

nd of the main Ioop._

CALL RESTART: write out current

restart and diagnostic

results to disk.

GEOS-11ARIES GCM FLOW CHART, page 7.

THE FLOW CHART OF THE DYNAMICAL CORE PACKAGE

DESCRIBING BOTH THE FIRST STEP (PREDICTOR) AND

THE SECOND STEP (CORRECTOR) OF THE MA TSUNO

TIME INTEGRATION SCHEME AS WELL AS THE

LEAPFROG TIME INTEGRATION SCHEME.

Put controltotal variable diabatic tendency tendencies terms. into the 1I

CALL SHAPIJ: apply the global Shapiro low-pass

filter scheme on U, V, T, q fields to

damp small-scale dispersive waves.

If integrating backward in time, invert the sign of the time [

I

tendency terms to ensure non-reversible effects are positive.j

CALL SETGRID: define parameters for integrating]

at the staggered grid points, j

GEOS-11ARIES GCM FLOW CHART, page 8.
CALL SETDMP: calculate damping coefficient_

for high latitude filter, j

CALL PKAP:

Compute perturbation

Average mass to

calculate Phillips ' 'P**KAPPA'" on C-grid]

geopotential height.]

vorticity points. ]

CAAL SUBI: compute kinetic energy, potential

vorticity, ustar, vstar, which are

used for Arakawa C-grid scheme.

GEOS-11ARIES (_CM FLOW CHART, pago 9.

CALL SUB2:compute which are parameter used for Arakawa alpha, beta,C-gridgamma,scheme.delta

CALL HADVECT: compute tendencies of height

and wind due to the horizonal

advection processes.

CALL HADVCTT: compute temperature

tendencies due to the horizonal

advection processes.

CALL HADVCTT: compute moisture ]

tendencies due to the horizonal

advection processes. J

Compute total pressurethe adiabatic tendencies. pressure and ]

_J_ GEOS-1/ARIES GCM FLOW CHART, pago 10.

10
CALL GETOMEGA: compute omega di'agnostic. }

Compute PI*SIGMADOT]

CALL VADVCT: calculate vertical centered second-order } advection of U, V, T, q.

CALL FFTDDT: apply FFT scheme to filter the

tendencies of U, V, T, q over the

high-latitude region (polar filter).

_J_ GEOS-I/ARIES GCIM FLOW CHART, pmge 11.

Add omega analysis diagnostics.increment to dynamical]

CALL STEP: update prognostic fields one timestep, compute total tendency

diagnostics, check global mean

surface pressure and negative

humidities, bump diagnostic

counters, as well as CALL

TMFILT for applying the

Asselin time filter.

EXIT

I

!

GEOS-11ARIES GCM FLOW CHART, pago 12.

11
3 Tangent linear model of the adiabatic version of NASA

GEOS-1 C-Grid GCM

3.1 Linearized discrete dynamical equations

The linearized discrete dynamical equations of Eqs. (8)- (14) used in the derivation of the

tangent linear model of the adiabatic version of the NASA GEOS-1 GCM are as follows

(we use { } to describe the basic state trajectory terms and ()' to denote the perturbation

variables terms. For all other definitions of the symbols used please see the documentation

of the ARIES/GEOS Dynamical Core (Suarez and Takacs, 1994)).

The linearized equations of the hydrostatic equations (8)- (9) are

((_LM) t _'_ (_s)t + Cp{OLM}(PLM+I -- PLM)' + Cp{/)LM+, -- PLM}(OLM) ' (15)

and

(_)' = (_+_)' + cp{Oz+,}(P_+l- P_)'+cp{P_+_- pt)(O_+_)',

The linearized continuity equations are

forl= 1, LM- 1 (16)

alld

LM 1 .

O(lri,j)' _ __, 2--[6i(u_)' + 6j(v I )']i,3(6a)t (17)

Ot _=1(Ap),,j

O(_j)' (lS)

{Tri,j,l+l}((Ti,j,l+l). t _- {_i'J'lT1)(Tri'J'l+l)t-_ --0"/+1 tOt

t 1

- Z (_)_,j[_,(uT)'+ _j(v;)']_,_(_)t

/=1

Tile linearized form of the thermodynamic equation (12) is

0({.}(0_)')_,_ +

1

[_({_*}(O")'+ ff"}(_')')z +

(z_i),,j

+6j({_'}(OJ)' + {O_}(,*)')&j -

_({_}(o)' + {_0}(o)' + {0_}(_)')

(_la i,j,l

Ot Ot

The linearized component forms of the momentum equation (13)- (14) are

0('_i,j,l)' 1 ), ,, ,

(19)

(20)

12
O(v<j,l)'

Ot

* t * 1.[/_i "it q-{/_i,j}(Yi,j+l) -_ {Vi,j+lJ \ ,3] q-

_-{"/i,j}(Vi:j)' _ - {Vi,j}("fi,j) "_ {(Si,j}(Yi+l,j) _- {Vi+l,ji((Si,j) ! _-

+{"_,_}(_b-_ 1' + {_7,j-1}(-_,j)'- {._,j+l}(_b+, )' - {_b+l }("_,J+,1'-

* * ! * t *

--{Q+I,j}(Ui+Io ')'-- {Ui+I,j}(f-i+I,j ) _- {(i,j}(Ui_l,j) "Jr {?li_l,j}((i,j)']l _-

+ { 1 1 [_i(u/ -- UI-1)-_ (TC(Y)l+l('tl+l .4 - Ul)ji,j] }(#i], , ,i,j -(_)_,j(_)_ 2

1 1 [{7--_., )' {ut (_i) {(Tri)i,j}(_lO.) l _ (_'0")/} (Ul -- UI-1 -_ -- ?t/-I }

1 i [{7----_-_., "Ill)t+ {?/l+l ?/l} ((71"CY=_')i+1) t] {(TCi)i,j}(_lO.) l -_ (TrO')l+l} (?'//+l -- -- --

1 [ { (dP) ii_i(_r),] _

(A'*x)i,j _d'_,+ h'_)' +cp O_\d_r ]_ i,j

1 cp6i{Tr} Ol for j= 1,JM- 1

(/x_x)_,j \ dTr/ t ] i,j

1

(_kVy)j [{0_i--1 ,j--1 }(?Z;-1,j--1 )t .__ {U__I,j_ 1 }(O:i_l,j_ ' )1 q_ (21)

-_-{/_i,j-1}(_Z;,j_l)t q- {U;,j_l}(/_i,j-1)t q-

"_-{Ti,j}(tli,j) "_- {Ui,j)(Ti,j ) -_- {_i_l,j}(Ui_l,j) t "_- {Ui_l,j}(_i_l,j)t-_

_-{#i,j}(Vi*Wl,j )l _-

-{_,j}(v;,j+,)' -

+ (eJ)_,jt_t_)_2

_ 1 1

{(_J)_,j}(azo)t2

_ 1 1

{(_J)_,j}(&.)_2

(zxvy)¢ aJ(_ + I¢_)'+ _ 0,

\d_r'/l ] Ji,j

* I * I * .)! { Yi+l ,j } (t*i,j) - {.__,,,}(v,__,.) - { Yi-1 ,j } ( #i- 1,3

{ Vi:j+ 1 }( Crgi,j )l ___ { _tgi,j_ 1 }( Vi:j_l )t .._ { Vi:j_ 1 }(Ctgi,j_l )! _1_

7:vJ "J I "

[(TrO.)l(v I _ Vl_l ) _1_ (71-o-)1_1_1 (VI+ 1 __ yl)]i,j (_3):,)_

)

(V' -- Vl-1)t q- {Vl -- Vl-1} (-(_) i,j

[{_-)J+l } (Y/+I- Vl) t-_- {Vl+I- Vl} ((71"0":---_")]+l)/J --

i,j

7r ,]i,j

forj=2, JM- 1

tlere the a, /_, 7, a, e, qo, u and # denote linear combinations of neighboring potential

vorticities.

13
To obtain the perturbation temperature output fields (T)', with the definition of 8,

{0}+ (e)' - {r} + (T)'

{P} + (P)' (22)

we have

(T)' = ({O} + (0)')({P} + (P)')- {0}{P}

= [({8} + (8)')({ps} _ + ((ps)') _) - {0}{ps} _] P t (23)

where we take PT = 0 and use Equation (27a) of Suarez and Takacs (1994)'s documentation

of the dynamical core of NASA model,

Pl -- _r'_Pl (24)

where

We may design the TLM version of NASA GEOS-1 GCM from the above TLM discrete

equations. However, to code the TLM conveniently and to avoid coding mistakes, we choose

another way to code the tangent linear model, i.e., we linearize the original adiabatic version

of NASA GEOS-1 GCM code segment by segment. The detailed method is presented in

the next subsection.

3.2 Coding 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 iineaxized nonlinear forward model in the vicinity of a basic

state which is on a model trajectory. For any original code fine, we may write it as

U = f(X) (26)

where

X = (xl, x2, ".., xm) T (27)

where 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 xl, x2, "-', xm (the components

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

14
original modelcontrol variables or of the intermediate variables derived from the original

control variables, m is the number of the required variables.

The corresponding tangent linear code will assume the form:

_U = _Xl /f\(O-h-_} +$x2 if\|O.____} +...+ (28)

kOXl/

where X = Xb_sic state means that in the expression °-1-,ox,i = 1, 2, .. ., m, all the values of

the required variables Xl, x2, ..., Xm are chosen to have the exact same values as those of

the basic state trajectory values 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/fU and bxl, dfx2, ..., _xm are the corresponding perturbation

variables of U and Xl, x2, .'., xm, respectively.

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

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

consists of calculating in parallel the nonlinear model integration trajectory as the basic

state Xbasic state and carrying out the integration of perturbation variables(such as dfU) in

the tangent linear model.

The dynamical core part of the tangent linear model flow chart of the NASA GEOS-1

GCM is presented in Subsection 3.4 (for the adiabatic version one needs to omit some steps

related to the computation of physical processes, moisture processes and some diagnostic

calculations),

3.3 Notational convention for variables and subroutines used in the tangent linear model code

For convenience, the same original names used in the nonlinear forward model are employed

for the corresponding perturbation variables in the tangent linear model code. For instance,

we use "U" for "_fU", "PKHT" for "(5(PKHT)", "USTR2" for "_f(USTR2)", etc.. This

also means that the perturbation control variables in the TLM share the same common

structure and same common block names as the GCM itself. So one needs to pay attention

to this issue when running the TLM in conjunction with the original GCM.

We just add a "0" at the end of a variable name in the original nonlinear forward model

to represent the corresponding basic state variable, such as using "U0" for "Ubasic state",

"PKHTO" for "( PKHT)basic state" , "USTR20" for "(gSTR2)bas_c state", etc..

For naming subroutines in the tangent linear model , we simply add a "L" at the beginning

15
of the original names of subroutines of the nonlinear forward model. To conform to the

general FORTRAN language rule, if the new name of a tangent linear subroutine exceeds

six letters, we just retain its first six letters. For instance, for the subroutines of the

original nonlinear model "SUB1", "VADVCT" and "HADVCTT", the corresponding

names of the subroutines in the tangent linear model code are "LS'UBI", "LVADVC" and

"LH ADVC", respectively.

3.4 Flow chart of the dynamical core of the tangent linear model

The following is the flow chart of the dynamical core of the tangent linear model of NASA

GEOS-1 GCM.

THE FLOW CHART OF THE TANGENT LINEAR MODEL OF

THE DYNAMICAL CORE PACKAGE DESCRIBING

THE TWO STEPS OF THE MA TSUNO TIME INTEGRATION

SCHEME AND THE LEAPFROG TIME INTEGRATION SCHEME.

Put total diabatic tendencies into the ]

/

perturbation variable tendency terms. /

CALL LSI-IAPl:apply the global Shapiro filter

on U, V, T, q perturbation fields to

damp small-scale dispersive waves.

If integrating backward, invert the sign of the time tendency]

/

terms to ensure non-reversible effects are positive. J

CALL SETGRIO: define parameters for integrating] at the staggered grid points.

TLM OF THE DYNAMICAL CORE, pego 1,

16
CALL SETDMP: calculate damping coefficient_

for high latitude filter. J

CALL LPKAP: calculate perturbation "'P**KAPPA'" term.}

Compute perturbation geopotential height.)

Average perturbation mass to vort. points. ]

CAAL LSUBI: compute the perturbation terms,

kinetic energy, potential

vorticity, ustar, vstar, which are

used for Arakawa C-grid scheme

TLM OF THE DYNAMICAL CORE, pago 2.

CALL SUB2:compute parameters alpha, beta, gamma, delta

which are used for Arakawa C-grid scheme.

CALL LHADVE: compute tendencies of perturbations

of height and wind due to the

horizonal advection processes.

CALL LHADVC" compute perturbation temperature:

tendencies due to the horizonal

advection processes.

CALL LHADVC: compute perturbation moisture

tendencies due to the horizonal

advection processes (tk_r the

adiabatic version, skip this step).

[ Compute the adiabatic perturbation pressure] and total perturbation pressure tendencies. J

_/ TLM OF THE DYNAMICAL CORE, page 3.

17
Compute delta(PI*SIGMADOT) ]

CALL LVADVC: calculate centered second-order ]

vertical advection of perturbations of

U, V, T, q.

CALL FFTDDT: apply FFT scheme to filter the

tendencies of the perturbation

variables of U, V, T, q over the

the high-latitude region.

TLM OF THE DYNAMICAL CORE, pBge 4.

Add analysis increment to dynamical

omega diagnostics (for the adiabatic

tangent linear model, it is not

necessary to compute diagnostics.

May skip this step).

CALL LSTEP: update perturbation prognostic fields

one time-step, CALL LTMFIL

for applying Asselin time filter.

(The following computations are

skipped in the adiabatic version of

TLM, i.e., compute total perturbatio[

tendency diagnostics, check global

mean surface pressure and negative

humidities, as well as bump

diagnostic counters.)

EXIT

I

!

TLM OF THE DYNAMICAL CORE. page 8.

18
3.5 Verifying correctness of the tangent linear model

rib verify the correctness of the tangent linear model, we compared the output of each

subroutine of the tangent linear model with its counterpart in the original forward model.

To verify the full tangent linear model, we employed a more quantitative method, described

below.

The evolution of X, the vector of control variables, is given by the integration of the model

M between times to and tn as:

X(tn) = M (tn,to)(X (to)) = M (tn,to) (Xo (to) + 6X (to)) (29)

whereas the first order evolution of the perturbation _X (tn) is the result of the integration

of the tangent linear model R :

(tn) = R (tn, to) (to)

We then compare the total perturbation

N (_X (to)) = M (t_,to) (Xo (to) + _X (to)) - M (tn,to) (Xo (to))

with its linear component

L (_iX (to)) = R(tn,to)(_X (to)

The difference between the two is denoted as

D (SX (to)) = N (SX (to)) - L (_X (to))

(30)

(31)

(32)

(33)

In order to quantify this comparison, we choose a norm whose square is defined by

flXlf2 = xTwx (34)

in actor<lance with the norm used in tire inner product of tile cost function for the variationM

data assimilation problem. The relative difference between the tangent linear model and

lhe nonlinear forward model is then defined as the ratio _[. We first, examine different

components of [LD_]_[ and calculate correlation coefficients between nonlinear output fields IILII

-Y and linear output fields L according to the individual model variables contributions

(_s, v, 7', l's).

The data used to verify tile tangent linear model is the Jaimary 1, 1985 00Z ECMWF data.

As in tlahier and Courtier(1992), we chose zonal average fields as basic state initial condi1.ion, while departure of zonal average tields multiplied by a. scaling factor a serves as the

19
perturbation of initial conditions. Both the basic state initial condition fields (i.e., zonal

average fields) and the perturbed initial condition fields (which is just the original ECMWF

data) are good initial conditions for stable integration of the model and do not cause excitation of gravity waves. The amplitudes of different components of this perturbation are

very large when a = 1: at level a = 0.223, the maximum zonal wind perturbation reaches

63 m/s; at level a = 0.352 , the maximum meridional wind perturbation reaches 65 m/s;

the temperature perturbation at a = 0.029 level is close to 18 K, while the maximum perturbation of the surface pressure is -424 hPa, due to the orography of the Tibetan Plateau.

The tangent linear model check with these strongly perturbed initial conditions can shed

light on possible coding errors. The period of integration, tn, is taken to be 12 hours.

Table 1: Correlation Coefficients Between N Field and L Field:

c_

1.0

10-1

10-2

10-3

10-4

10-5

u v T Ps

0.8919275 0.9094758 0.8628350 0.8714384

0.9981103 0.9986473 0.9971051 0.9981210

0.9999804 0.9999861 0.9999694 0.9999809

0.9999998 0.9999999 0.9999997 0.9999998

1.0000000 1.0000000 1.0000000 1.0000000

1.0000000 1.0000000 1.0000000 1.0000000

1.0

10-1

10-_

10 -3

10-4

lO-S

Table 2: Relative Error _ (%):

u v T ps

45.32 41.99 50.69 49.41

6.18 5.20 7.62 6.14

0.63 0.53 0.79 0.62

6.33 × 10-2 5.27 × 10 -2 6.22 × 10 -_

6.33 x 10 -3 5.27 × 10 -3

7.88 × 10-2

7.89 × 10 -3 6.22 × 10 -3

6.33× 10-4 5.27x 10 -4 7.89× 10 -4 6.22x 10 -4

Table 1 presents correlation between N field and L field for various values of parameter

while Table 2 displays relative error _-_. From these tables we see that all correlation

coefficients between N and L for each of the variable fields exceed 86% (when a=l.0), and

reach values close to unity when a is less or equal to 0.1. As a decreases, the relative

errors decrease to very small values in a linear manner. The correlation coefficients reach

up to 10 digits of accuracy in vicinity of unity while relative error values attain an order of

magnitude of 10-4 when ct is equal to 10 -5. Comparing with similar relative errors analysis

applied to adiabatic version of NASA/GLA Semi-Lagrangian Semi-Implicit (SLSI) GCM

('Fable 3) (Li et al., 1994), the tangent linear model of the adiabatic version of NASA

2O
GEOS-1C-Grid GCM seemsto displaybetter linearity propertythan the tangentlinear

modelof NASA/GLA SLSIGCM. Theseanalysesprovidea reliableindicationabout the

correctnessof the tangentlinearmodelcode.

Table3: RelativeError _ (%) of NASA/GLA 3-DSLSIGCM(fromLi et al., 1994):

O¢

1.0

10-1

10-2

10 -3

10 -4

10-5

U v

34.85 42.09

Z y

T

44.12

(lnps)'

28.77

5.22 10.63 7.63 3.44

6.05 4.81 26.40 3.44

8.85 x 10 -2 0.14 0.11 5.41 × 10 -2

2.77x 10 -2 4.14x 10 -2 2.67x 10 -2 1.49x 10 -2

1.71x 10 -2 2.16x 10 -2 1.38x 10 -2 1.34x10 -2

To assess impact of the length of integration period (or the length of data assimilation

window) on the validity of tangent linear approximation, we carried out a check for different

window lengths, up to 96 hours. For convenience, we averaged correlation coefficients of

four model variables and used the norm (described in Eq. (3.6) in Li et al., 1994) for HD[I

and [[L[[ to calculate relative error. The diagonal component values of weighting matrix

used are W_, = Wv = 10 -3 82rn -2, WT = 10 -1 K -2, Wps=lO -3 hPa -2, respectively.

We chose values of the parameter a=l.O, as representing a strong perturbation, c_=0.1, as

representing a normal perturbation, while a=0.01 represented a small perturbation.

Fig. 1 displays the correlation coefficients between the N and the L fields for three values

of the parameter a with respect to different lengths of integration period. Fig. 2 displays

• oikOA

the relatwe error II.r_llcurves. Considering the ce = 0.1 curve we deduce that the error of the

tangent linear mollie is small when a normal perturbation is used corresponding to perturbations of wind, temperature and surface pressure of lm/s, 1K and lOhPa, respectively.

Even for integration periods of up to 96 hours, the correlation coefficient still exceeds 86%

and the relative error is 57%. These numerical results confirm earlier results of Courtier and

Talagrand (1987), Lacarra and Talagrand (1988), Rabier and Courtier (1992) that tangent

linear model well approximates the nonlinear forward model for up to 4 days when initial

perturbations are not too strong. For the weak perturbation case (a = 0.01), the model

displays a very good linear behavior with a correlation coefficient reaching 99.6% while the

relative error remains below 10% after 4 days of numerical integration. For strongly perturbed initial conditions(a = 1), the validity of the tangent linear model decreases quickly as

the length of integration period is increased. With these strong perturbations, the validity

limit of the tangent linear model is less than a day.

21
I 00

•95

.90

.85

80

.75

70

.65

•55

50

45

.4_

\,

\

l i l I

10 20 30

\

-\

"\

\

"\

-\

"\

p\

l I I f I _t r I i t

40 50 60 70 80 90

T_me i_ Hours

1.00

95

90

.B5

,80

75

,70

.b5

60

.55

.50

45

.40

100

Figure 1: Correlation coefficient between the N fields and the L fields with respect to

integration periods. Curve 1 corresponds to a -- 1, curve 2 corresponds to a = 0.1, and

curve 3 corresponds to a = 0.01.

100 - ,

90

80

70

bO

50

40

30

20

10

0

0

/

/

/

/

/

,<

7-

• _"l_

_.1.s

10 20 30 40 50 60 70

Time in Hours

3 _3----3--"

Llll_

80 90

100

90

80

70

60

50

40

30

20

10

0

100

Figure 2: Relative error IIDII/IINI with respect to integration periods. Curve 1 corresponds

to a = 1, curve 2 corresponds to a = 0.1, and curve 3 corresponds to a = 0.01.

22
4 Adjoint model of the adiabatic version of the NASA GEOS1 C-Grid GCM

4.1 lTsing adjoint method to calculate the gradient of the cost function

The practical determination of the adjoint model of the adiabatic version of the NASA

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

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

carrying out the 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 describing the numerical weather prediction model

as strong constrains.

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

to the observations, it assumes the form :

1 R

E (35)

where X(tr) is a model state vector of size M(3K + 1) containing the values of the zonal

wind u, the meridional wind v, the temperature T and the surface pressure P_. M is the

number of grid points at each level; K is number of vertical levels, tr is a given time in the

assimilation window; x°b_(tr) is a vector of observations defined over all grid points on all

levels at time tr; W(tr) is an N × N diagonal weighting matrix. From Navon et a1.(1992),

we have the following expression

R

xv0): Z

rm0

x'(t_). (36)

where X'(t0) is the initial perturbation, X'(t_) is the perturbation in the forecast resulting

from 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) (37)

where P_ represents the result of applying all the operator matrices in the linear model to

obtain X'(tr) from X'(t0).

23
We definethe adjointmodelas

X_(t0) = PT:K(tr), r-= 1,...,R, (38)

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

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

conditions is

R

VJ (X(t0)) = _ pTw(t_) (X(t_)- X°b'(t_)) (39)

r----O

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

of the tangent linear model operator.

The flow chart of the adjoint model (of the dynamical core of the adiabatic version of NASA

GEOS-1 GCM is presented in Subsection 4.4.

4.2 Coding 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 - AIA2. ". AN, (40)

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 = A TNATN-1 ""AT" (41)

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

linear model. Each of the DO loops and each of the subroutines in the tangent linear

model have their adjoint image DO loop and subroutine. Therefore, we code the adjoint

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

linear model sentence by sentence (i.e., on DO loop by one DO loop, subroutine by subroutine) 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

24
adjointmodelcode.Forinstance,weuse"U" for"/)", "PKHT" for "(PA:HT)", "USTR2"

for "(US(FR2)",etc.. As in the TLM, this conventionalsomeansthat the adjoint control

variablesin the adjoint modelsharethe samecommonstructureandsamecommonblock

namesastheGCM itself. Sooneneedsto payattentionto it whencarryingout four dimensionaldataassimilationexperimentsin whichthe adjointmodelwill berun in conjunction

with the originalGCM.

We also just add a "0" at the end of a variablename(in a similar way as donepreviouslyin the tangentlinear model)to representhe correspondingbasicstate variable,

suchasusing "U0" for "Ub_slcstat¢","PKHTO" for "(PKHT)b_sicstat¢","USTR20" for

"( U ST R2 )basic state", etc..

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 to the general FORTRAN language rule. For instance, for the subroutines of

the original nonlinear model "SUB 1", "VADVCT" and "H ADVCTT", the corresponding

names of the subroutines in the tangent linear model code are "LSUBI", "LI/ADVC" and

"LH ADVC", respectively, and the corresponding adjoint subroutine names are "ASU B 1",

"AVADI/C" and "AHADVC", respectively.

4.4 Flow chart of the dynamical core of the adjoint model

In this subsection, we present the flow chart of the dynamical core of the adjoint model of

NASA GEOS-1 GCM.

25
THE FLOW CHART OF THE ADJOINT MODEL OF THE

DYNAMICAL CORE PACKAGE DESCRIBING THE

TWO STEPS OF THE MA TSUNO "riME INTEGRATION SCHEME

AND THE LEAPFROG TIME INTEGRA T/ON SCHEME. All

COMPUTATIONS RELA TED TO RETRIEVE THE BASIC STATE

TRAJECTORY AND INI'rlALIZA TION OF VARIABLES AND

PARAMETERS ARE NOT INCLUDED.

CALL ASTEP: update adjoint prognostic fields

backward one time-step, and

CALL ATMFIL for applying

the adjoint computations of Asselin

time filter.

CALL AFTDDT:

CALLAVADVC:

[

CALL AHADVC:

|

apply the adjoint computation of the[

FF'F scheme to the tendencies of

the adjoint variables U, V, T, q

over the high-latitude region.

ADJOINT OF THE DYNAMICAL CORE, page 1.

compute the adjoint operations of

the centered second-order vertical

advections of perturbations of U, V,

T,q.

Compute the adjoint of the delta(PI*SIGMADOT) _

J

Compute the adjoint of the total perturbation

pressure tendencies and the adiabatic perturbation

pressure tendencies.

compute the adjoint operations of

the perturbation moisture tendencies

due to the horizonal advection

processes (for the adiabatic version,

skip this step).

ADJOINT OF THE DYNAMICAL CORE, page 2.

26
CALL AHADVC: compute the adjoint operations of

the perturbation temperature

tendencies due to the horizonal

advection processes.

CALL AHADVN: compute the adjoint operations of

the tendencies of perturbations of

height and wind due to the horizonal

advection processes.

CALL ASUB2: compute the adjoint operator

related to parameters alpha,

beta, gamma, delta which are

used in the Arakawa C-grid

scheme.

ADJOINT OF THE DYNAMICAL CORE, pago 3,

CAAL ASUBI: compute adjoint operator related

to perturbation terms, kinetic

energy, potential vorticity,

ustar, vstar, which are used in

the Arakawa C-grid scheme.

perturbation Adjoint computations mass to vorticity of average points. ]

Adjoint perturbation computations geopotential of theheight. i

CALL APKAP: apply the adjoint operations for calculating]

the perturbation "'P**KAPPA" term. j

/

I_ ADJOINT OF THE DYNAMICAL CORE, DegQ 4.

27
when integrating backward (in time)

in TLM, invert the sign of the time

tendency terms in the adjoint model

to ensure that the non-reversible

processe are positive.

CALL ASHAPI: apply the adjoint operations for

the global Shapiro low-pass filter

on U, V, T, q perturbation fields

to damp small-scale dispersive

waves.

( EXIT )

!

!

ADJOINT OF THE DYNAMICAL CORE, page 5.

4.5 Verification of the correctness of the adjoint model

Integrating 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 Eq.

(39)), one can obtain value of gradient of cost function with respect to distributed control

variables, which may consist of either the initial conditions or the initial conditions plus

boundary conditions or model parameters. Since the adiabatic version of NASA GEOS-1

C-Grid GCM consists of thousands of lines of code, any minor coding error may cause

the final gradient of cost function with respect to the control variables to be erroneous.

Therefore, we need to verify the correctness of the linearization and adjoint coding segment

by segment. Each segment may consist of either a subroutine or of several DO loops. For

a detailed derivation of the adjoint model and verification of its correctness, see Navon et

a1.(1992).

The correctness of the adjoint of each operator was checked by applying the following

identity (Navon et al., 1992)

(AQ).T(AQ) = Q.T (A.T(AQ)), (42)

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

a subroutine (See Navon et al. 1992). The left hand side involves only the tangent linear

code, while the right hand side involves also adjoint code (A'T). If Eq. (42) holds, tile

adjoint code is correct when compared with the TLM. In practice identity Eq. (42) holds

only up to machine accuracy. In our verifications of the correctness of each segment of

the ad.joint model and the whole adjoint model, the LIIS and the RIIS of Eq.Eq. (42)

attained 13 digits of accuracy which is near the machine accuracy limit of NASA Charney

C-90 ('.RAY Computer which has intrinsic double precision.

In the subroutine "ASHAP", to ensure the correctness of the adjoint code check, we employ

the "DOUBLE PRECISION" definition on some variables. In practice, to save computational cost, if do not need higher accuracy results of the adjoiut model, one nlay omit use

the "DOUBLE PRECISION" definition.

These results show that our adjoint code consists of absolutely the exact adjoint operators

of the TLM of the adiabatic version of NASA GEOS-1 C-Grid GCM.

A gradient check (Fig. 31) was then performed to assess accuracy of the discrete adjoint

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

follows:

1 R

r_0

where X(t_.) is an N = (M(3/_" + 1)) component vector containing values of In, v, T,

Ps), with which NASA GEOS-1 GCM model is initialized, over all grid points and at all

vertical levels at time t_; M is the number of grid points at each vertical level; /_" is the

number of vertical levels; R is the number of time levels for the analyzed fields in the

assimilation window; Z_ is a certain observation time in the assimilation window; x°b'(l,.)

is the N-component vector of analyzed values of X over all grid points on all levels at

time t,,; and W(tr) is an N × N diagonal weighting matrix, where W_,, W_,, W'r and

Wp,_ are diagonal submatrices consisting of weighting factors for each variable, respectively.

Their respective values (as used in gradient check calculation) were W,, -- 10-3I s_m. -2,

W_, --- 10-:_I s';m -2, Wy = 10-1I I( -2, Wp_= 10-'_I hPa -'2. Then, let

J(X + (_h) = J(X) + (_hTVj(X) + O(a2), (44)

be a Taylor expansion of the cost fimction, ttere c_ is a small scalar and h is a vect.or of

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

of(I 8S

qs(g) = J(X + c_h) - J(X)

c_hTVJ(X ) = 1 + O(g). (45)

29
1.15

1 10

1.05

1 .00

95

90

85

80

II .75

-10 -q =8 -7 -6 -5 -4 3 2

Log t&tpKal

I I ' 1 _1 15

110

1 a5

1 00

_5

911

85

80

i L _ ±__L._ 75

-f 0

Figure 3: Variation of the _b(a) with respect to log a (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 ECMWF data was used as t = 0 observations. The first guess

is the shifted 6-hour initial condition.

Figure 4: Variation of the log I¢(a) - 11 with respect to logo. Integration period is 6 hours

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

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

3O
For valuesof a which are small but not too close to the machine zero, one should expect

to obtain a value for ¢(t_) which is close to unity. We obtained that the value of function

¢(a) equals unity to a high degree of accuracy when parameter a varied from 10 -1 to 10-s,

and obeys a monotonically decreasing rule when a decreased over 12 orders of magnitude.

From the residual of ¢(c_)(Fig. 4), we found that the residual tends linearly to zero. The

gradient check verifies that adjoint model is correct and can be safely used to perform 4-D

VDA experiments.

4.6 Some examples of coding the discrete adjoint model from the tangent

linear model code

In the following we will provide some simple tutoring coding examples and technical methods

which are very helpful for understanding the tangent linear model and the adjoint model

coding techniques used in the adiabatic version of NASA GEOS-1 C-Grid GCM.

• EXAMPLE 1:

In the original nonlinear model

DO 10 I=1, N-1

10 X(I)=A*Y(I+I)

where A is a parameter, X and Y are N dimensional vectors. This DO loop is linear and do

not need to be linearized. So the corresponding tangent linear code will remain identically

the same as the original one,

DO I0 I=l, N-I

i0 X(I)=A*Y(I+I)

The corresponding adjoint code form depends on whether the values of Y(1) will be reused

or not after this DO loop. If tim values of Y(I) will not be reused after this DO loop, the

31
matrix form of this DO loop is

x(1)

x(2) I =

X(N- 1) '

A 0 ... 0 0

0 A ... 0 0

0 0 ... 0 A

Y(2)

Y(3)

Y(N)

The adjoint of the equation (46) can be written out directly as

r(2)

Y(3)

Y(N)

A 0 0 ... 0 O\

J

0 A O ...0 0

..... . .

0 0 O ... O A

x(a) )

x(2)

X(N)

Equation (47) is equivalent to the following code

(46)

(47)

DO I0 I=l, N-I

I0 Y(I+I)=A*X(I)

However, If tile values of Y(I) will be reused after this EXAMPLE 1 DO loop, the corresponding matrix form should modified to be

x(1)

x(2)

X(N- l)

g(2)

Y(a)

Y(N)

A 0 0 --- 0 0

0 A 0 ".. 0 0

, , • • . . ,

0 0 0 ... 0 A

1 0 0 ... 0 0

0 1 0 -.- 0 0

..... ,

0 0 0 ... 0 1

[Y(2) /

Y(3)

Y(N)

(48)

32
Theadjoint of the equation (48)shouldbewritten out as

Y(2)

Y(3)

Y(N)

A

0

X(1)

x(2)

o o ... o o a o o ... o o\

/

A 0 ... 0 0 0 1 0 ... 0 0

............... X(N- 1)

¥(2)

Y(3)

0 0 .-- 0 A 0 0 0 ... 0 1

Y(N)

Thus the adjoint code of Equation (49) will be

(49)

DO I0 I=l, N-I

10 Y(I+I)=Y(I+I)+A*X(I)

Thus this is a different adjoint code compared to the code derived from Equation (/refeq47).

The issue of identifying which variables on the right side of tangent linear code belong either

to reused variables or non-reused variables is very important for adjoint code derivation.

• EXAMPLE 2:

In the original nonlinear model code and the tangent linear model code, a DO loop assumes

the following form

DO I0 I=I, N-I

I0 X(I)=X(I) +A*Y(I+I)

The matrix form of the above DO loop (if the values of Y(/) will not be reused) for the

33
tangentlinearmodelis

x(a)

x(2)

X(N - 1)

1 0 ... 0 0 A 0 -.. 0 0

0 1 ... 0 0 0 A ... 0 0

0 0 ... 0 1 0 0 ... 0 A

X(1)

x(2)

X(N - a)

v(2)

Y(3)

Y(N)

(50)

The adjoint operation matrix equation will be (by transposition, i.e., for tangent linear code

X = Lf'(, the adjoint code is _[ = LTx.)

x(1)

x(2) ,

X(N- 1)

V(2)

Y(3)

Y(N)

1 0 ... 0 0

0 1 ... 0 0

...... i

0 0 ... 0 1

A 0 ... 0 0

0 A ... 0 0

0 0 ... O A

x(x)

x(2)

X(N - 1)

(51)

So the corresponding adjoint code is

DO I0 I=_-, N-_

10 Y(I+I)=A*X(I)

Similarly, we may derive easily the adjoint code for this example when Y(I) will be reused

after this DO loop in the tangent linear model as

34
DO 10 I=l, N-I

10 Y(I+I)=Y(I+I)+A*X(I)

• EXAMPLE 3:

In the nonlinear model, a DO loop is

DO 10 I=1, N

I0 Z(I)--X(I)*Y(I)

where X, Y and Z are N dimensional vector variables. The corresponding tangent linear

model (:ode is

DO 10 I=I, N

10 Z(I)=YO(I)*X(I)+IO(I)*Y(I)

Recall from Subsections 3.3 and 4.3 that Y0 and X0 are basic state variables. If Y(I) and

X(I) will not be reused anymore after this DO loop, the matrix form of this DO loop is

35
Z(1)

z(2)

z(u)

YO(1)

0

0

X(1)

X(2)

X(N)

#

Y(1)

Y(2)

Y(i)

The adjoint of the equation

x(1)

x(2)

X(N)

Y(a)

Y(2)

Y(N)

0

YO(2)

• . 0 X0(1) 0

•. o o xo(2)

• . YO(N) 0 0

(52) can be written out directly as

YO(1) 0

0 YO(2)

0 0

XO(1) 0

o xo(2)

o o

YO(N)

0

0

XO(N)

I

The adjoint DO loop code is

• . 0

• " 0

• .

• . XO(N)

Z(1) /

Z(2)

Z(N)

(52)

(53)

10

DO 10 I=1, N

X(I)=YO(I)*Z(I)

Y(I)=XO(I)*Z(I)

36
If Y(I) and X(I) will be reused in future, the matrix form of the linearized code of EXAMPLE 3 is

z(_)

z(2)

Z(N)

X(a)

X(2)

x(N)

r(2)

_'(N)

(YO(1) 0

0 YO(2)

0 0

1 0

0 1

0

0

0

0

X(1)

X(2)

X(N)

r(1)

Y(2)

Y(N)

0

o xo(_) o

o o xo(2)

YO(N) 0 0

0 0 0

0 0 0

1 0 0

0 1 0

• . 0 0 1

• ° °

• " " 0 0

0

0

XO(N)

0

0

0

0

"* 0

• ° °

• ' " 1

(54)

37
Thenthe adjointmatrix formis

X(1) ro(1)

x(2) o

X(N)

r(1)

Y(2)

Y(N)

0

X0(1)

0

0

Z(1)

z(2)

0

Y0(2)

0

0

X0(2)

1 0 -.. 0 0 0 .-- 0

0 1 --. 0 0 0 ... 0

Y0(N) 0 0 --- 1 0 0 "" 0

o o o ... o a o ... o

0 0 0 ... 0 0 1 ... 0

XO(N) 0 0 ... 0 0 0 ... I

Z(N)

X(1)

X(2)

X(N)

Y(_)

Y(2)

Y(N)

Thus the adjoint code can be obtained from above equation as

(55)

DO i0 I=l, N

X(I)=X(I)+YO(I)*Z(I)

I0 y(I)=Y(I)+XO(I)*Z(I)

• EXAMPLE 4:

38
In last three examples, all the calculations do not have recursion relation. Here we present

a simple example which exhibits a recursion relation. The DO loop code of the original

nonlinear model and its tangent linear model is

DO 10 I=l, N-I

10 Z(I+I)=Z(I)+Y0(I)*X(I)+XO(I),Yy

where Z(I), X(I) and YY are control variables and will be reused in future, Y0 and X0

are two parameter vectors. This recursion loop can be expressed in matrix form as

z(1)

z(2)

z(:{)

z(N- _)

Z(N)

X(1)

1 0 0 ... 0 0 0 0 ... 0 0 0

0 1 0 -.. 0 0 0 0 ..- 0 0 0

0 0 ] --. 0 0 0 0 ... 0 0 0

0 0 0 ... 0 1 0 0 --. 0 0 0

0 0 0 .-. 0 1 0 0 ---0 YO(N-1) XO(N-1)

0 0 0 .-. 0 0 1 0 ... 0 0 0

0 0 0 ... 0 0 0 t ... 0 0 0

0 0 0 ... 0 0 0 0 ... 0 1 0

0 0 0 --. 0 0 0 0 ... 0 0 1

39
1 0 0 ---0 0 0 0 ...

0 1 0 ... 0 0 0 0 .-.

0 0 1 ... 0 0 0 0 ..-

0 0 0

0 0 0

0 0 0

0 0 0 -.- 1 0 0 0 ---Y0(N-2) 0 X0(N-2)

0 0 0 ... 0 0 0 0 -.. 0 0 0

0 0 0 ... 0 0 1 0 --- 0 0 0

0 0 0 ... 0 0 0 1 .-- 0 0 0

0 0 0 .-. 0 0 0 0 .-- 0 1 0

0 0 0 .--0 0 0 0 ..- 0 0 1

1 0 0 .-. 0 0 0 0 ... 0 0 0

1 0 0 --. 0 0 YO(1) 0 ---0 0 XO(1)

0 0 0 --- 0 0 0

0 0 0 ... 0 0 0

0 0 0 ... 0 0 0

0 0 0 -.. 0 0 1

0 0 0 ... 0 0 0

0 0 0 --. 0 0 0

0 0 0 -.- 0 0 0

0 ... 0 0 0

0 --- 0 0 0

0 ... 0 0 0

0 ..- 0 0 0

1 ..- 0 0 0

0 --- 0 1 0

0--- 0 0 1

:_ °°°_

z(x)

z(2)

Z(N - 2)

Z(N- 1)

X(1)

X(2)

X(N- 1)

YY

(56)

40
Thecorrespondingadjointmatrix equationis

z(1)

z(2)

Z(N - 2)

Z(N - 1)

X(1)

X(2)

X(N- 1

Y Y

1 1

0 0

0 0

0 0

0 0

0 ... 0 0 0 0 --- 0 0

0 -.. 0 0 0 0 ... 0 0

0 ... 0 0 0 0 ... 0 0

0 ... 0 0 0 0 ... 0 0

0 ..- 0 0 0 0 --- 0 0

o Yo(1) o ... o o 1 o ... o o

0 0 0 ... 0 0 0 1 ... 0 0

0 0 0 ... 0 0 0 0 ... 0 0

0 0 0 "" 0 0 0 0 ... 1 0

0 XO(1) 0 -.-0 0 0 0 ---0 1

1 0 0

0 1 1

0 0 0

0 0 0

0 0 0

0 0 0

• .. 0 0 0 0 ... 0 0

• .. 0 0 0 0 ... 0 0

• .. 0 0 0 0 ... 0 0

• .. 0 0 0 0 ... 0 0

• -- 0 0 0 0 -.. 0 0

• .. 0 0 1 0 ... 0 0

o o vo(2) ... o o o 1 ... o o

0 0 0

0 0 0

• .. 0 0 0 0 ... 0 0

• .. 0 0 0 0 ... 1 0

o o xo(2).., o o o o ... o 1

41
100 ... 0 0

010 ... 0 0

00 1 ... 0 0

.... .° •

.... o ° °

.... o . •

000 ... 0 0

000 -.. 1 1

• ""* 000 ... 0 0

000 ... 0 0

.... • . °

000 ... 0 0

000 ... 0 YO(N-1)

0 0 0 ... 0 XO(N-1)

00 ... 00

00 ... 00

00 --. 00

00 -.. 00

00 .-- 00

10..-00

0 1 -.. 00

00 ... 00

00-.. 10

00 ... 01

Finally, the adjoint code for this recursion DO loop is

z(2)

z(3)

Z(N- 1)

Z(N)

X(1)

x(2)

X(N- 1)

YY

(57)

DO 19 I=N-I-, I-,-IZ(I)=Z(I)+Z(I+I9

X (I)--X(I)+YO (I)*Z(I+ I)

YY=YY+XO (I)*Z (I+19

10 Z(I+_)=0.0

For this example, the DO loop variable I must evolve in the opposite direction of the

variation of the DO loop variable in the tangent linear model code. As a matter of fact,

the adjoint model is always integrated in the opposite direction of the tangent linear model.

Thus, to avoid unexpected mistakes in the adjoint code and to conveniently identify and

code the adjoint model, we follow the convention that the DO loop variables in all DO loops

of adjoint model code should evolve in the opposite direction of the corresponding tangent

linear model DO loops.

• EXAMPLE 5:

Here we present a simple example for a subroutine, the original code of which is

42
SUBROUTINE SIM(X, Y, N)

DIMENSION X(N), Y(N), Z(N*2)

Y(N)=O .0

DO I0 1=1-, N-I

I0 Y(I)=X(I+I)**2

DO 20 I=l, N

20 Z(I)=X(I)**3

DO 30 I=l. N

30 Z(I+N)=X(I)**4

DO 40 I---1-, N

40 Y(I) =Y(I) +I0. O*Z (I) +8. O*Z (I+N)

CALL OTHER(Y, N)

RETUKN

END

In this subroutine, the input variable is "X", and the output variables are "X" and "Y".

Its tangent linear code is

I0

2O

SUBROUTINE LSIM(X, Y, N, XO, YO)

DIMENSION X(N), Y(N), Z(N*2), XO(N), YO(N),

YCN)=O.O

YO(N)=O.O

DO 10 I=I-; N-_

Y(I):2.0*XO(I÷I)*X(I÷I_

YO(I)=XO(I+I)**2

DO 20 I=I-; N

Z(I)=3.0*XO(I)**2*X(I)

ZO(I)=XO(I)**3

DO 3O I=l, N

Z(I÷N)=4.0*XO(I)**3*X(I)

ZO(N*2)

43
30 ZO(I+N)=XO(I)**4

DO 40 I=I, N

y(1)=Y(1)+iO.O*Z(I)+8.0*Z(I+N)

40 YO(I)=YO(I)+IO.O*ZO(I)+8.0*ZO(I+N)

CALL LOTHER(Y, N, YO)

RETURN

END

The corresponding adjoint subroutine code(assume X(I) will be reused in future) is

SUBROUTINE ASIM(X, Y, N, XO, YO)

DIMENSION X(N), Y(N), Z(N*2), XO(N), YO(N), ZO(N*2)

CALL AOTHER(Y, N, YO)

DO 40 I=N, I, -I

Z(I)=IO.O*Y(I)

40 Z(I+N)=8.0*Y(I)

DO 30 I=N, 1, -I

30 X(I)=X(I)+4.0*XO(I)**3*Z(I+N)

DO 20 I=N, 1, -I

20 X(I)=X(I)+3.0*XO(I)**2*Z(I)

DO 10 I=N-I, I, -1

I0 X(I+I)=X(I+I)+2.0*XO(I+I)*Y(I)

RETURN

END

• EXAMPLE 6:

Finally, we provide here a real subroutine in the adiabatic version of the NASA GEOS-1

GCM. This subroutine relates the horizontal advection calculations of momentum equations

using the Arakawa-Lamb c-grid energy-conserving form.

44
-----The original nonlinear forward model subroutine code is----

SUBROUTINE HADVECT ( DUDT, DVDT, CONV,

FACTU, FACTV, FACTH, A,B,C,D,

ZKE, Q, PHI, PKZ,

USTAR,VSTAR,TBARU,TBARV, IM,JM,

UDOT,VDGT )

DIMENSION DUDT(IM,JM),

DIMENSION A(IM,JM), B(IM,JM),

DIMENSION ZKE(IM,JM), Q(IM,JM),

DIMENSION USTAK(IM,JM), VSTAR(IM,JM)

DIMENSION TBARU(IM,JM), TBARV(IM,JM)

DIMENSION UDOT(IM,JM), VDOT(IM,JM)

DIMENSION FACTU(IM,JM)

DIMENSION FACTV(IM,JM)

DIMENSION FACTH(IM,JM)

DVDT(IM,JM), CDNV(IM,JM)

C(IM,JM), D(IM,JM)

PHI(IM,JM), PKZ(IM,JM)

PARAMETER ( ZERO = 0.00 )

PARAMETER ( AHALF = O.SO )

PARAMETER ( ONE = 1.00 )

PARAMETER ( TWO = 2.00 )

PARAMETER ( THREE = 3.00 )

PARAMETER ( FOUR = 4.00 )

IMJM = IM* JM

IMJMMI = IM*(JM-I)

IMJMM2 = IM*(JM-2)

IMJMM3 = IM*(JM-3)

IMJMM4 = IM*(JM-4)

CP = GETCON('CP')

J = 2

DO I=i, IMJMM2

CONV(I,J) = - (USTAR(I,J) - USTAR(I-i,J)

+ VSTAR(I,J) - VSTAR(I,J-I) )

45
ENDDO

C FIX LONGITUDINAL BOUNDARIES

C

IMI = IM

I = I

DO J=2,JM-I

CONV(I,J) = - (USTAR(I,J) - USTAR(IMI,J)

+ VSTAR(I,J) - VSTAR(I,J-I) )

ENDDO

J = 2

DO I=I,IMJMM2

UDOT(I,J) = A (I ,J ) * VSTAR (I+1,J )

+ B (I ,J ) * VSTAR (I ,J )

+ C (I ,J ) * VSTAR (I ,J-l)

+ D (I ,J ) * VSTAR (I+l,J-l)

- (ZKE(I+Z,J ) + PHI (I+l,J ) )

+ (ZKE(I ,J ) + PHI (I ,J ) )

-CP*TBARU(I,J)* (PKZ(I+I,J) - PKZ (I ,J ) )

ENDDO

C FIX LONGITUDINAL BOUNDARIES

C

DO J=2,JM-I

IMI = IM-I

I = IM

IP_ = 1

UDOT(I,J)

-CP*TBARU(I,J)* (PKZ(IPI,J

= A (I ,J ) * VSTAR (IPl,J )

+ B (I ,J ) • VSTAR (I ,J )

+ C (I ,J ) * VSTAR (I ,J-I)

+ D (I ,J ) * VSTAR (IPi,J-l)

- (ZKE(IPI,J) + PHI (IPI,J) )

+ (ZKE(I ,J ) + PHI (I ,J ) )

) - PKZ (I ,J ) )

IMI = IM

I = I

IPI = 2

46
UDOT(I,J)

-CP*TBARU(I,J), (PKZ(IP_J ) - PKZ

ENDDO

= A (I ,J ) * VSTAR(IPI.J)

+ B (I ,J ) * VSTAR(I ,J )

+ C (I ,J ) * VSTAR(I ,J-l$

+ D (I ,J ) * VSTAR(IPl, J-l)

- (ZKE(IPI,J) + PHI (IPi,J ) )

+ (ZKE(I ,J ) + PHI (I ,J ) )

(I ,J ) )

J=JM-I

DOI=I, IM

VDOT(I, J)=ZER0

ENDDO

J=2

DO I=I, IMJMM3

VDOT(I,J)

ENDDO

= - C (I ,J+l) * USTAR(I ,J+l)

- D (l-l,J+l) * USTAR(l-l,J+l)

- A (I-1,J ) * USTAR (I-1,J )

- B (I ,J ) * USTAR (I ,J )

- ( ZKE(I ,J+l) + PHI (I ,J+l) )

+ ( ZKE(I ,J ) + PHI (I ,J ) )

-CP*TBARV(I,J)*( PKZ(I ,J+l) - PKZ (I ,J ) )

DO J=2,JM-2

IMI = IM

I = 1

IPI =2

VDOT(I,J)

-CP*TBARV(I,J)*( PKZ(I

= - C (I ,J+l) * USTAR (I

- D (IMI,J+I) * USTAR (IMi

- A (IMI,J ) * USTAR (IMI

- B (I ,J ) * USTAR (I

- ( ZKE(I ,J+l) + PHI (I

+ ( ZKE(I ,J ) + PHI (I

,J+l) - PKZ (I

J+i)

J+l)

J )

J )

J+l) )

J ) )

J ) )

IMI = IM-I

I = IM

4?
IPI = 1

VDOT(I,J) = - C (I

- D (IMi

- A (IMI

-B (I

- ( ZKE(I

+ ( ZKE(I

-CP,TBARV (I, J)* ( PKZ(I

ENDDO

J+l) * USTAR (I ,J+l)

J+l) * USTAR (IMi,J+l)

J ) * USTAR (IMI,J )

J ) * USTAR (I ,J )

J+l) + PHI (I ,J+l) )

J ) + PHI (I ,J ) )

J+l) - PKZ (I ,J ) )

DO I=I,IMJMM2

CONV(I,2) = CONV(I,2) * FACTH(I,2)

DUDT(I,2) = DUDT(I,2) + UDOT(I,2) * FACTU(I,2)

DVDT(I,2) = DVDT(I,2) + VDOT(I,2) * FACTV(I,2)

ENDDO

DO I=I, IMJM

UDOT(I, I)=ZERO

VDOT(I, 1)=ZERO

END DO

RETURN

END

48
--Its corresponding TLM subroutine code is accordingly--

SUBROUTINE LHADVE ( DUDT, DVDT, CONV,

FACTU, FACTV, FACTH, A,B,C,D,

ZKE, PHI, PKZ,

USTAR,VSTAR,TBARU,TBARV, IM,JM,

UDOT,VDOT,

DUDTO, DVDTO, CONVO,

AO,BO,CO,DO,

ZKEO, PHIO, PKZO,

USTARO,VSTAR0,TBARU0,TBARVO,

UDOTO,VDOT0)

DIMENSION

DIMENSION

DIMENSION

DIMENSION

DIMENSION

DIMENSION

DUDT(IM,JM),

DUDTO(IM,JM),

A(IM,JM),

AO(IM,JM),

ZKE(IM,JM),

ZKEO(IM,JM),

DVDT(IM,JM), CONV(IM,JM)

DVDTO(IM,JM), CONVO(IM,JM)

B(IM,JM), C(IM,JM),

BO(IM,JM), CO(IM,JM),

PHI(IM,JM), PKZ(IM,JM)

PHIO(IM,JM),PKZO(IM,JM)

DIMENSION USTAR(IM,JM), VSTAR(IM,JM)

DIMENSION USTARO(IM,JM), VSTARO(IM,JM)

DIMENSION TBARU(IM,JM), TBARV(IM,JM)

DIMENSION TBARUO(IM,JM), TBARVO(IM,JM)

DIMENSION UDOT(IM,JM), VDOT(IM,JM)

DIMENSION UDOTO(IM,JM), VDOTO(IM,JM)

DIMENSION FACTU(IM,JM)

DIMENSION FACTV(IM,JM)

DIMENSION FACTH(IM,JM)

PARAMETER ( ZER0

PARAMETER ( AHALF

PARAMETER ( ONE

PARAMETER ( TWO

PARAMETER ( THREE

PARAMETER ( FOUR

= o.oo )

: o.5o )

: t.oo )

: 2.00 )

= 3.00 )

= 4.00 )

IMJM = IM* JM

IMJMM2 = IM*(JM-2)

CP = GETCON('CP')

DO I=I, IMJM

D(IM,JM)

DO(IM,JM)

49
UDOT(I, I)=ZER0

VDOT(I, I)=ZER0

ENDD0

DOJ = 2, JM-I

D0 I=2,1M

CONV(I,J) = - (USTAR(I,J) - USTAR(I-I,J)

+ VSTAR(I,J) - VSTAR(I,J-I) )

CONVO(I,J)= - (USTARO(I,J) - USTARO(I-I,J)

+ VSTAR0(I,J)- VSTAR0(I,J-I) )

ENDD0

ENDD0

J=JM-I

DOI=I, IM

VDOT(I, J)=ZERO

ENDDO

C FIX LONGITUDINALBOUNDARIES

C

IMI = IM

I = 1

DOJ=2,JM-I

CONV(I,J) = - (USTAR(I,J) - USTAR(IMI,J)

+ VSTAR(I,J) - VSTAR(I,J-I) )

CONVO(I,J)= - (USTARO(I,J) - USTAR0(IMI,J)

+ VSTAR0(I,J)- VSTAR0(I,J-_) )

ENDD0

DOJ=2,JM-I

DO I=_,IM-_

UDOT(I,J)

A (I

+ AO(I

+ B (I

,J ) * VSTARO(I+I,J )

,J ) * VSTAR (l+l,J )

,J ) * VSTARO(I ,J )

50
+ BO(I ,J ) * VSTAR (I ,J )

+ C (I ,J ) * VSTARO(I ,J-l)

+ CO(I ,J ) * VSTAR (I ,J-l)

+ D (I ,J ) * VSTARO(I+I,J-I)

+ DO(I ,J ) * VSTAR (I+l,J-1)

- (ZKE(I+I,J ) + PHI (I+l,J ) )

+ (ZKE(I ,J ) + PHI (I ,J ) )

-CP*(TBARU(I,J)* (PKZO(I+I,J) - PKZO (I ,J ) )

+TBARUO(I,J)* (PKZ(I+l,J ) - PKZ (I ,J ) ))

UDOTO(I,J) = AO(I ,J ) • VSTARO(I+I,J )

+ BO(I ,J ) • VSTARO(I ,J )

+ CO(I ,J ) * VSTARO(I ,J-l)

+ DO(I ,J ) * VSTARO(I+I,J-I)

- (ZKEO(I+I,J ) + PHIO (l+l,J ) )

+ (ZKEO(I ,J ) + PHIO (I ,J ) )

-CP*TBARUO(I,J)* (PKZO(I+i,J ) - PKZO (I ,J ) )

ENDDO

ENDDO

C FIX LONGITUDINAL BOUNDARIES

C ...........................

I = IM

IPI = i

DO J=2,JM-I

UDOT(I ,J)

A (I ,J ) • VSTARO(IPI,J )

+ AO(I ,J ) • VSTAR (IPI,J )

+ B (I ,J ) • VSTARO(I ,J )

+ BO(I ,J ) • VSTAR (I ,J )

+ C (I ,J ) • VSTARO(I ,J-l)

+ CO(I ,J ) * VSTAR (I ,J-l)

+ D (I ,J ) * VSTARO(IPI,J-I)

+ DO(I ,J ) * VSTAR (IPI,J-I)

- (ZKE(IPI,J ) + PHI (IPI,J

+ (ZKE(I ,J ) + PHI (I ,J

-CP,(TBARU(I,J)* (PKZO(IPI,3 ) - PKZO (I

+TBARUO(I,J)* (PKZ(IPI,J ) - PKZ

UDOTO(I,J) = AO(I ,J

+ BO(I ,J

+ CO(I ,J

+ DO(I ,J

- (ZKEO(IPI,J

+ (ZKEO(I ,J

) )

))

,J ) )

(i ,J ) ))

) • VSTARO(IPI,J )

) • VSTARO(I ,J )

) • VSTARO(I ,J-l)

) • VSTARO(IP1,J-I)

) + PHIO(IPI,J ) )

) + PHIO (I ,J ) )

51
-CP.TBARUO(I,J)* (PKZO(IPI,J) - PKZO (I ,J) )

ENDDO

DO 3=2, JM-2

DO I=2,IM

VDOT(I,J)

ENDDO

ENDDO

- CO(I ,J+l) * USTAR (I ,J+IO

- C (I ,J+I9 * USTARO(I ,J+19

- DO(I-_,J+I) * USTAR (I-I.J+l)

- D (I-_,J+l) * USTARO(I-I%J+I)

- AO(I-I-.J ) * USTAR (I-_,J)

- A (I-I.J) * USTARO(I-I_J )

- BO(I ,J ) * USTAR (I ,J )

- B (I ,J ) • USTARO(I ,J )

- ( ZKE(I ,J+1) + PHI (I ,J+19 )

+ ( ZKE(I ,J ) + PHI (I ,J ) )

-CP.(TBARVO(I,J)*( PKZ(I ,J+l) - PKZ (I ,J ) )

+TBARV(I,J)*( PKZO(I ,J+t) - PKZO (I ,J ) ))

VDOTO(I,J) = - CO(l ,J+l) * USTAR0(I ,J+19

- DO(I-r,J+I) * USTARO(I-I.J+I)

- AO(l-r,J ) * USTARO(I-_,J )

- BO(I ,J ) • USTARO(I ,J )

- ( ZKEO(I ,3+19 + PHIO (I ,J+I9 )

+ ( ZKEO(I ,J ) + PHIO (I ,J ) )

-CP*TBARVO(I,J)*( PKZO(I ,J+19 - PKZO (I ,J ) )

IMr = IM

I = 1-

DO J--2,JM-2

VDOT(I ,J)

- CO(l ,J+IO * USTAR (I ,J÷l$

- C (I ,J+l) * USTARO(I ,J+l-)

- DO(IM_,J+IO * USTAR (IM_,J+IO

- D (IMI_J+IO * USTARO(IM_,J+I9

- AO(IMr, J ) * USTAR (IM_,J)

- A (IMr, J ) * USTARO(IMr, J )

- SO(I ,J ) • USTAR (I ,J )

52
ENDDO

- B (I ,J ) * USTARO(I ,J )

- ( ZKE(I ,J+l) + PHI (I ,J+l) )

÷ ( ZKE(I ,J ) + PHI (I ,J ) )

-CP*(TBARVO(I,J)*( PKZ(I ,J+19 - PKZ (I ,J ) )

+TBARV(I,J)*( PKZO(l ,J+l) - PKZO (I ,J ) ))

VDOTO(I,J) = - CO(I ,J+l) * USTARO(I ,J+l)

- DO(IMI,J+I) * USTARO(IMI,J+I)

- AO(IMI,J ) * USTARO(IMI,J )

- BO(I ,J ) * USTARO(I ,J )

- ( ZKEO(I ,J+l) + PHIO (I ,J+1) )

+ ( ZKEO(I ,J ) + PHIO (I ,J ) )

-CP*TBARVO(I,J)*( PKZO(I ,J+l) - PKZO (I ,J ) )

DO I=I,IMJMM2

CONV(I,2) = CONV(I,2) * FACTH(I,2)

CONVO(I,2) = CONVO(I,2) * FACTH(I,2)

DUDT(I,2) = DUDT(I,2) + UDOT(I,2) * FACTU(I,2)

DUDTO(I,2) = DUDTO(I,2) + UDOTO(I,2) * FACTU(I,2)

DVDT(I,2) = DVDT(I,2) + VDOT(I,2) * FACTV(I,2)

DVDTO(I,2) = DVDTO(I,2) + VDOTO(I,2) * FACTV(I,2)

ENDDO

DD I=I, IMJM

UDOT(I, 1)=ZERO

VDOT(I, I)=ZERO

END D0

RETURN

END

53
----.The corresponding adjoint subroutine code is accordingly--

SUBROUTINE AHADVE ( DUDT, DVDT, CONV,

FACTU, FACTV, FACTH, A,B,C,D,

ZKE, PHI, PKZ,

USTAR,VSTAR,TBARU,TBARV, IM,JM,

UDOT,VDOT,

AO,BO,CO,DO,

PKZO,

USTARO,VSTARO,TBARUO,TBARVO)

DIMENSION DUDT(IM,JM),

DIMENSION A(IM,JM),

DIMENSION AO(IM,JM),

DIMENSION ZKE(IM,JM),

DIMENSION PKZO(IM,JM)

DVDT(IM,JM), CONV(IM,JM)

B(IM,JM), C(IM,JM),

BO(IM,JM), CO(IM,JM),

PHI(IM,JM), PKZ(IM,JM)

DIMENSION USTAR(IM,JM), VSTAR(IM,JM)

DIMENSION USTARO(IM,JM), VSTARO(IM,JM)

DIMENSION TBARU(IM,JM), TBARV(IM,JM)

DIMENSION TBARUO(IM,JM), TBARVO(IM,JM)

DIMENSION UDOT(IM,JM), VDOT(IM,JM)

DIMENSION FACTU(IM,JM)

DIMENSION FACTV(IM,JM)

DIMENSION FACTH(IM,JM)

D(IM,JM)

DO(IM,JM)

PARAMETER ( ZERO = 0.00 )

PARAMETER ( AHALF = 0.50 )

PARAMETER ( ONE = 1.00 )

PARAMETER ( TWO = 2.00 )

PARAMETER ( THREE = 3.00 )

PARAMETER ( FOUR = 4.00 )

IMJM = IM* JM

IMJMM2 = IM*(JM-2)

CP = GETCON('CP')

C **** APPLY SCALE FACTORS ****

54
DOI=I, IMJM

UDOT(I, 1)=ZERO

VDOT(I, 1)=ZERO

ENDDO

DOI=IMJMM2,I, -I

VDOT(I,2)=DVDT(I,2)• FACTV(I,2)

UDOT(I,2)=DUDT(I,2)• FACTU(I,2)

CONV(I,2)=CONV(I,2), FACTH(I,2)

ENDDO

DOI=l, IMJM

ZKE(I, I)=ZER0

PHI(l, I)=ZER0

A(I, 1)=ZERO

B(I, 1)=ZERO

C(I, i)=ZERO

m(I, I)=ZERO

ENDDO

IMI = IM

I = I

DO J=JM-2, 2, -I

AMIO=-CP*VDOT(I,J)

AM20=AMIO*TBARVO(I,J)

TBARV(I,J)=AMIO,(PKZO(I ,J+l) - PKZO (I ,J ) )

PKZ(I ,J+I)=PKZ(I ,J+I)+AM20

PKZ (I ,J )=PKZ (I ,J )-AM20

ZKE(I, J)=ZKE(I, J)+VDOT(I,J)

ZKE(I ,J+I)=ZKE(I ,J+I)-VDOT(I,J)

PHI (I ,J )=PHI (I ,J )+VDOT(I,J)

PHI (I ,J+I)=PHI (I ,J+I)-VDOT(I,J)

B (I ,J )=B (I ,J )-VDOT(I,J) , USTARO(I,J )

USTAR(I ,J )=USTAR(I ,J )-BO(I ,J ) * VDOT(I,J)

A (IMI,J )=A (IMI,J )-VDOT(I,J) • USTARO(IMI,J)

USTAR(IMI,J )=USTAR(IMI,J )-A0(IMI,J ) • VDOT(I,J)

D (IMI,J+I)=D (IMI,J+I)-VDOT(I,J) • USTARO(IMI,J+I)

55
USTAR(IMI,J+I)=USTAR(IMI,J+_-DO(IMI,J+I) * VDOT(I,J)

C (I ,J+l)=C (I ,J+I)-VDOT(I,J) * USTARO(I ,J+l)

USTAR (I ,J+I)=USTAR (I ,J+l)-CO(l ,J+l) * VDOT(I,J)

END D0

DO J=JM-2, 2, -I

DO I=IM, 2, -I

AMIO=-CP*VDOT(I,J)

AM20=AMIO*TBARVO(I,J)

TBARV(I,J)=AMIO,( PKZO(I ,J+l) - PKZO (I ,J ) )

PKZ(I ,J+I)=PKZ(I ,J+I)+AM20

PKZ (I ,J )=PKZ (I ,J )-AM20

ZKE(I, J)=ZKE(I, J)+VDOT(I,J)

ZKE(I ,J+I)=ZKE(I ,J+I)-VDOT(I,J)

PHI (I ,J )=PHI (I ,J )+VDOT(I,J)

PHI (I ,J+I)=PHI (I ,J+I)-VDOT(I,J)

B (I ,J )=B (I ,J )-VDOT(I,J) * USTARO(I ,J )

USTAR (I ,J )=USTAR (I ,J )-BO(I ,J ) • VDOT(I,J)

A (l-l,J)=A (l-l.J)-VDOT(I,J) * USTARO(I-I,J )

USTAR (l-l,J)=USTAR (I-1,J)-AO(I-I.J ) * VDOT(I,J)

D (I-I,J+I)=D (I-I,J+I)-VDOT(I,J) * USTARO(I-I,J+I)

USTAR (I-I.J+I)=USTAR (I-I.J+I)-DO(I-I,J+I) * VDOT(I,J)

C (I ,J+l)=C (I ,J+t)-VDOT(I,J) * USTARO(I ,J+l)

USTAR (I ,J+I)=USTAR (I ,J+l)-CO(I ,J+l) * VDOT(I,J)

END D0

END D0

C FIX LONGITUDINAL BOUNDARIES

C ....

I = IM

IPI = 1

DO J=JM-I, 2, -I

AMIO=-CP*UDOT(I,J)

AM20=AMIO*TBARUO(I,J)

PKZ(IPI,J )=PKZ(IPr, J )+AM20

PKZ (I ,J )=PKZ (I ,J )-AM20

TBARU(I,J)=AMIO*(PKZO(IPI;J ) - PKZO (I ,J ) )

ZKE(I ,J )=ZKE(I ,J )+UDOT(I,J)

PHI (I ,J )=PHI (I ,J )+UDOT(I,J)

ZKE(IPI,J )=ZKE(IP_,J )-UDOT(I,J)

PHI (IPI,J)=PHI (IPl,J)-UDOT(I,J)

56
VSTAR (IPI.J-I)=VSTAR (IPI,J-I)+DO(I ,J ) * UDOT(I,J)

D (I ,J )=D (I ,J )+VSTARO(IP1,J-I) • UDOT(I,J)

VSTAR (I ,J-_)=VSTAR (I ,J-I)+CO(I ,J ) • UDOT(I,J)

C (I ,J )=C (I ,J )+VSTARO(I ,J-l) * UDOT(I,J)

VSTAR (I ,J )=VSTAR (I ,J )+BO(I ,J ) * UDOT(I,J)

B (I ,J )=B (I ,J )+VSTARO(I ,J ) * UDOT(I,J)

VSTAR (IPI;J)=VSTAR (IPI,J)+AO(I ,J ) * UDOT(I,J)

A (I ,J )=A (I ,J )+VSTARO(IPI,J ) * UDOT(I,J)

END DO

C **********************************************************************

C   COMPUTE U-WIND TENDENCIES ** **

DO J=JM-I, 2, -i

DO I=IM-f, T, -I

AMIO=-CP*UDOT(I,J)

AM20=AMIO*TBARUO(I,J)

PKZ(I+I.J )=PKZ(I+I.J )+AM20

PKZ (I ,J )=PKZ (I ,J )-AM20

TBARU(I,J)=AMIO*(PKZO(I+I,J ) - PKZO (I ,J ) )

ZKE(I ,J )=ZKE(I ,J )+UDOT(I,J)

PHI (I ,J )=PHI (I ,J )÷UDOT(I,J)

ZKE(I+I,J )=ZKE(I+I.J )-UDOT(I,J)

PHI (I+l,J )=PHI (I+l,J )-UDOT(I,J)

VSTAR (I+I,J-I)=VSTAR (I+I,J-I)+DO(I ,J ) * UDOT(I,J)

D (I ,J )=m (I ,J )+VSTARO(I+I,J-I) * UDOT(I,J)

VSTAR (I ,J-I)=VSTAR (I ,J-I)+CO(I ,J ) * UDOT(I,J)

C (I ,J )=C (I ,J )+VSTAmO(I ,J-i) * UDOT(I,J)

VSTAR (I ,J )=VSTAR (I ,J )+BO(I ,J ) * UDOT(I,J)

B (I ,J )=B (I ,J )÷VSTARO(I ,J ) * UDOT(I,J)

VSTAR (I+_,J )=VSTAR (I+1,J )+AO(I ,J ) * UDOT(I,J)

A (I ,J )=A (I ,J )+VSTARO(I+t,J ) * UDOT(I,J)

END DO

END DO

C FIX LONGITUDINAL BOUNDARIES

C ...........................

IMI = IM

I = 1

DO J=JM-1, 2, -iUSTAR(I, J)=USTAR (I, J) -CONV (I, J)

USTAR (IMI. J) =USTAR (IMi, J) +CONV (I, J)

VSTAR (I, J) =VSTAR (I, J) -CONV(I, J)

57
VSTAR(I,J-19=VSTAR(I,J-19+CONV(I,J)

ENDDO

C **********************************************************************

C   COMPUTEHEIGHTTENDENCIES ** _*

DOJ=JM-_, 2, -_

DOI=IM, 2, -_

USTAR(I,J)=USTAR(I,J)-CONV(I,J)

USTAR(I-I_J)=USTAR(I-_,J)+CONV(I,J)

VSTAR(I,J)=VSTAR(I,J)-CONV(I,J)

VSTAR(I,J-19=VSTAR(I,J-19+CONV(I,J)

ENDDO

ENDDO

D_ I=_, IMJM

UDOT(I, 19=ZERO

VDOT(I, _=ZEmO

END DO

RETURN

END

58
References

Arakawa, A., and V. R. Lamb, 1981: A potential enstrophy and energy conserving scheme

for the shallow water equations, Mon. Wen. Rev., 109, 18-36.

Asselin, R., 1972: Frequency filter for time integrations, Mon. Wea. Rev., 100, 487-490.

Brown, J. A., and K. Campana, 1978: An economical time-differencing system for numerical

weather prediction, Mon. Wen. Rev., 106, 1125-1136.

Burridge, D. M., and J. Haseler, 1977: A model for medium range weather forecasting

- adiabatic formulation, Tech. Report, No. _, European Center for Medium Range

Weather Forecasts, Bracknell, Berkshire, UK.

Courtier, P., and O. Talagrand, 1987: Variational assimilation of meteorological observations with the adjoint equation - Part I. Numerical results. Q. J. R. Meteorol. Soc.,

113, 1329-1347.

Fox-Rabinovitz, M., H. M. Helfand, A. Hou, L. L. Takacs, and A. Molod, 1991: Numerical

experiments on forecasting, climate simulation and data assimilation with the new 17

layer GLA GCM, Proceedings of the AMS Ninth Conference on Numerical Weather

Prediction, 21-25 October 1991, Denver, CO, 506-509.

Gates, W. L., 1992: AMIP: The Atmospheric Model Intercomparison Project, Bull. Ameri.

Meteor. Soc., 73, 1962-1970.

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 NumericaZ Weather Prediction, 21-25 October 1991, Denver, CO,

514-517.

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. Ameri. Meteor. Soc., 70,620-622.

Lacarra, J.-F., and O. Talagrand, 1988: Short-range evolution of small perturbations in a

barotropic model, Tellus, 40A, 81-95.

Li, Y., I. M. Navon, J. R. Bates, S. Moorthi, R. W. Higgins, X. Zou and W. Yang, 1993: Development of the Adjoint model of the Adiabatic Version of the 3-D NASA/GLA SemiLagrangian Semi-implicit GCM, Technical Report, FSU-SCRI-93-49, April 1993.

59
Li, Y., I. M. Navon,W. "fang,X. Zou,J. R. Bates,S. Moorthi,and R. W. Higgins,1994:

Four-DimensionalVariationalDataAssimilationExperimentswith a Multilevel semiLagrangiansemi-ImplicitGCM. Mon. Wea. Rev., 122, 5,966-983.

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

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

Rabier, F., and P. Courtier, 1992: Four-Dimensional assimilation in the presence of baroclinic instability, Q. J. R. Meteorol. Soc., 118, 649-672.

Sadourny, R., 1975: The dynamics of finite difference models of the shallow water equations,

J. Atmos. Sci., 32,680-689.

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

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

Schubert, S., M. Suarez, C. K. Park, and S. Moorthi, 1993: GCM simulations of interseasonal variability in the Pacific/North American region, J. Atrnos. Sci., 50, 1991-2007.

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

Core, Version 2, May 19, 1994.

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

General Circulation Model (GCM), Version 1, May, 17, 1994.

Yang, Weiyu, I. Michael Navon, 1995: A New Hessian Preconditioning Method Applied to

Variational Data Assimilation Experiments using Adiabatic Version of NASA/GEOS-1

GCM, Technical Report, FSU-SCRI-95-03, Supercomputer Computations Research

Institute, Florida State University, 43pp.

Yang, Weiyu, I. Michael Navon, 1995: Extending a Hessian Preconditioning Method to

4-D Variational Data Assimilation Experiments using the Adiabatic Version of the

NASA/GEOS-1 C-Grid GCM, WMO International Symposium on Assimilation of Observations in Meteorology and Oceanography, Tokyo, Japan, March 13-17, 1995, also

in the WMO International Workshop on Imbalances of Slowly Varying Components o/

Predictable Atmospheric Motions, Beijing, China, March 7-10, 1995.

Yang, Weiyu, I. Michael Navon, and Philippe Courtier, 1995: A New Hessian Preconditioning Method Applied to Variational Data Assimilation Experiments using the Adiabatic

Version of NASA/GEOS-1 GCM, Submitted to Mon. Wea. Rev..

6O
Previous Volumes in This Series

Volume 1

September 1994

Documentation of the Goddard Earth Observing System

(GEOS) general circulation model - Version 1

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

Volume 2

October 1994

Direct solution of the implicit formulation of fourth order

horizontal diffusion for gridpoint models on the sphere

Y. Li, S. Moorthi, and J.R. Bates

Volume 3

December 1994

An ancient thermal infrared radiation parameterization for

use in general circulation models

M.-D. Chou and M.J. Suarez

Volume 4

January 1995

Documentation of the Goddard Earth Observing System

(GEOS) Data Assimilation System - Version 1

James Pfaendtner, Stephen Bloom, David Lamich,

Michael Seablom, Meta Sienkiewicz, James Stobie,

and Arlindo da Silva

Volume 5

April 1995

Documentation of the Aries-GEOS dynamical core: Ve.rsion 2

Max J. Suarez and Lawrence L. Takacs

Volume 6

April 1995

A Multiyear Assimilation with the GEOS-1 System:

Overview and Results

Siegfried Schubert, Chung-Kyu Park, Chung-Yu Wu,

Wayne Itiggins, Yelena Kondratyeva, Andrea Molod,

Lawrence Takaes, Michael Seablom, and Richard

Rood

Volume 7

September 1995

Proceedings of the Workshop on the GEOS-1 Five-Year

Assimilation

Siegfried D. Schubert and Richard B. Rood

61


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 ot information. Send comments regarding this burden estimate or any other aspect of this

1. AGENCY USE ONLY (Leave b/ank) | 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

I

March 1996 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

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

Documentation of the Tangent Linear Model and Its Adjoint of the Adiabatic

Version of the NASA GEOS-I C-Grid GCM - Version 5.2 C - NAS5-32332

6. AUTHOR(S)

Weiyu Yang and I. Michael Navon Code 910.3

Max J. Suarez, Editor

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

Data Assimilation Office

Goddard Space Flight Center

Greenbelt, Maryland

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

National Aeronautics and Space Administration

Washington, DC 20546-0001

8. PEFORMING ORGANIZATION

REPORT NUMBER

96B00050

10. SPONSORING I MONITORING

ADGENCY REPORT NUMBER

TM-104606, Vol. 8

11. SUPPLEMENTARY NOTES

Yang and Navon: Supercomputer Computations Research Institute, Florida State University,

Tallahassee, Florida; Suarez: Goddard Space Flight Center, Greenbelt, Maryland

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

Unclassified - Unlimited

Subject Category 46

This publication is available from the NASA Center for AeroSpace

Information, 800 Elkridge Landing Road, Linthicum Heights, MD

21090-2934, (301)621-0390.

13. ABSTRACT (Maximum 200 words)

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

of the adiabatic version of the NASA GEOS-1 C-Grid GCM (Version 5.2) is presented. The derivations of and the methods for coding the TLM and its adjoint, as well as the notation conventions

used in these two models, are described in detail. The flow charts of the NASA GEOS-1 GCM, its

tangent linear model, and adjoint model are provided. The procedures and their results of correctness verification of the TLM and the adjoint model are presented. Finally, tutorial examples of

derivation of adjoint code from the tangent linear code are provided for the benefit of various users.

14. SUBJECT TERMS

Tangent linear model, TLM, NASA GEOS-1 C-Grid GCM (Version 5.2)

adjomt model, adiabatic version, derivations, methods for coding, notation

conventions, correctness verification tutorial

17. SECURITY CLASSIRCATION 18. SECURITY CLASSIRCATION

OF REPORT OF THIS PAGE

Unclassified Unclassified

NSN 7540-01-280-5500

19. SECURITY CLASSIRCATION

OF ABSTRACT

Unclassified

15. NUMBER OF PAGES

76

16. PRICE CODE

20. LIMITATION OF ABSTRACT

UL

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39.18

298-102

National Aeronautics and

Space Administration

Goddard Space Flight Center

Greenbelt, Maryland 20771

Official Business

Penalty for Private Use, $300

SPECIAL FOUR1H-Q.A_ RATE

POSTAGE & FEES PAID

NASA

PERMIT No. G27

POSTMASTER: If Undeliverable (Section 158.

Postal Manual) Do Not i:lmum