NCAR/TN-398 + STR

NCAR TECHNICAL NOTE

June 1995

A Description of the

Fifth-Generation Penn State/NCAR

Mesoscale Model (MM5)

Georg A. Grell1

Jimy Dudhia2

David R. Stauffer3

MESOSCALE AND MICROSCALE METEOROLOGY DIVISION

2NATIONAL CENTER FOR ATMOSPHERIC RESEARCH

BOULDER, COLORADO

1FORECAST SYSTEMS LABORATORY,

NATIONAL OCEANIC AND ATMOSPHERIC ADMINISTRATION

BOULDER, COLORADO

3DEPARTMENT OF METEOROLOGY,

THE PENNSYLVANIA STATE UNIVERSITY

UNIVERSITY PARK, PENNSYLVANIA
Table of Contents

Preface . . . . . . . . . . . . . . . . . . . . . . vii

Acknowledgments . . . . . . . . . . . . . . . . . ix

Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . 1

Chapter 2 Governing equations and numerical algorithms . . 1

2.1 The hydrostatic equations . . . . . . . . . . . . . . . 1

2.2 The nonhydrostatic equations . . . . . . . . . . . . . 3

2.2.1 Complete Coriolis force option . . . . . . . . . . . . . 5

2.3 Nonhydrostatic nite dierence algorithms . . . . . . . 6

2.4 hydrostatic nite dierence algorithms . . . . . . . . . 10

2.5 Time splitting . . . . . . . . . . . . . . . . . . . . 11

2.5.1 The nonhydrostatic semi-implicit scheme . . . . . . . . 11

2.5.2 The hydrostatic split-explicit scheme . . . . . . . . . . 13

2.6 Lateral boundary conditions . . . . . . . . . . . . . . 16

2.6.1 Sponge boundary conditions . . . . . . . . . . . . . . 16

2.6.2 Nudging boundary conditions . . . . . . . . . . . . . 16

2.6.3 Moisture variables . . . . . . . . . . . . . . . . . . 17

2.7 Upper radiative boundary condition . . . . . . . . . . 17

Chapter 3 The mesh-renement scheme . . . . . . . . . . . . 19

3.1 The monotone interpolation routines . . . . . . . . . . 19

3.2 Overlapping and moving grids . . . . . . . . . . . . . 23

iii
3.3 The feedback . . . . . . . . . . . . . . . . . . . . 23

3.3.1 A nine-point averager . . . . . . . . . . . . . . . . 24

3.3.2 A smoother-desmoother . . . . . . . . . . . . . . . . 24

Chapter 4 Four-dimensional data assimilation . . . . . . . . . 25

4.1 Analysis nudging . . . . . . . . . . . . . . . . . . . 26

4.2 Observational nudging . . . . . . . . . . . . . . . . 32

Chapter 5 Physical parameterizations . . . . . . . . . . . . . 37

5.1 Horizontal diusion . . . . . . . . . . . . . . . . . . 37

5.2 Dry convective adjustment . . . . . . . . . . . . . . 37

5.3 Precipitation physics . . . . . . . . . . . . . . . . . 38

5.3.1 Resolvable scale precipitation processes . . . . . . . . . 39

5.3.1.1 Explicit treatment of cloudwater, rainwater, snow, and ice . 39

5.3.1.2 Mixed-phase ice scheme . . . . . . . . . . . . . . . . 45

5.3.2 Implicit cumulus parameterization schemes . . . . . . . 46

5.3.2.1 The Kuo scheme . . . . . . . . . . . . . . . . . . . 46

5.3.2.2 A modied Arakawa-Schubert scheme . . . . . . . . . . 50

5.3.2.3 The Grell scheme . . . . . . . . . . . . . . . . . . . 64

5.3.3 Parameterization of shallow convection . . . . . . . . . 70

5.4 Planetary boundary layer parameterizations . . . . . . . 73

5.4.1 Surface-energy equation . . . . . . . . . . . . . . . . 73

5.4.1.1 Net Radiative 
ux Rn . . . . . . . . . . . . . . . . . 74

5.4.1.2 Heat Flow into the Substrate Hm . . . . . . . . . . . . 78

5.4.1.3 Sensible-Heat Flux Hs and Surface Moisture Flux Es . . . 78

5.4.2 Bulk-aerodynamic parameterization . . . . . . . . . . 78

5.4.3 Blackadar high-resolution model . . . . . . . . . . . . 80

5.4.4 Vertical diusion . . . . . . . . . . . . . . . . . . . 84

5.4.5 Moist vertical diusion . . . . . . . . . . . . . . . . 86

iv
5.5 Atmospheric radiation parameterization . . . . . . . . . 87

5.5.1 Longwave radiation . . . . . . . . . . . . . . . . . . 87

5.5.2 Shortwave radiative scheme . . . . . . . . . . . . . . 90

Appendix 1 Glossary of symbols . . . . . . . . . . . . . . . . 92

Appendix 2 Look-up table for transmissivity . . . . . . . . . . 111

Appendix 3 Map Pro jections . . . . . . . . . . . . . . . . . . 112

Appendix 4 Land-Use categories . . . . . . . . . . . . . . . . 116

References . . . . . . . . . . . . . . . . . . . . . 117

v
PREFACE

This technical report describes the fth generation Penn State/NCAR Mesoscale Model

Version 1, or MM5V1. It is intended to provide scientic and technical documentation

of the model for users. Source code documentation is available as a separate Technical

Note (NCAR/TN-392) by Haagenson et al. (1994). Comments and suggestions for

improvements or corrections, are welcome and should be sent to the authors.

(Last revision: June 1996)

vi
ACKNOWLEDGMENTS

A large number of scientists, students, and sponsors, too numerous to mention, have

made signicant contributions towards developing and improving this model. We would

like to thank everyone who directly implemented changes and corrections into the code, as

well as everyone who helped to identify errors and problems.

We would also like to extend our gratitude to the two reviewers (Prof. Tom Warner,

and Dr. Jian-Wen Bao), whose comments improved the manuscript signicantly.

Computing for this model development was undertaken on the CRAY YMP at NCAR,

which is supported by the National Science Foundation.

The writing of this technical document was in part supported by the Department

of Energy through grant DEA105-90ER61070, and the California Air Resources Board

through subcontract AG/990-TS05-S02 from Alpine Geophysics.

ix
1. Introduction

This technical report is a description of the fth-generation Penn State/NCAR

Mesoscale Model (MM5). It is based on the original version described by Anthes and

Warner (1978). Although a few of the following details of this model are well represented

in Anthes et al. (1987), extensive changes and increases in options have occurred. For

completeness, those parts that have changed little or none will also be represented here.

The document structure is as follows. In section 2 we will describe the governing equations,

algorithms, and boundary conditions. This will include the nite dierence algorithms

and time splitting techniques of both the hydrostatic and the nonhydrostatic equations

of motion (hydrostatic and nonhydrostatic solver). All subsequent sections will describe

features common to both solvers. Section 3 will discuss the mesh-renement scheme,

section 4 the four-dimensional data-assimilation technique, and section 5 will focus on the

various physics options.

2. Governing equations and numerical algorithms

2.1 Hydrostatic model equations

The vertical -coordinate is dened in terms of pressure.

 =

p  pt

ps  pt

;

where ps and pt are the surface and top pressures respectively of the model, where pt is a

constant.

The model equations are given by the following, where p = ps  pt .

Horizontal momentum;

@pu

@ t =  m

2 

@puu=m

@x +

@pvu=m

@y  

@pu_

@

 mp  



@p

@x +

@

@x  + p

f v + Du (2:1:1)

@pv

@ t =  m2 

@puv=m

@x +

@pvv=m

@y  

@pv_

@

 mp  



@p

@y +

@

@y   p

f u + Dv (2:1:2)

1
Temperature;

@pT

@ t =  m

2 

@puT =m

@x +

@pvT =m

@y  

@pT _

@

+ p

 !

cp

+ p

 Q_

cp

+ DT ; (2:1:3)

where the D terms represent the vertical and horizontal diusion terms and vertical mixing

due to the planetary boundary layer turbulence or dry convective adjustment. The heat

capacity for moist air at constant pressure is given by cp = cpd(1 + 0:8qv), where qv is

the mixing ratio for water vapor and cpd is the heat capacity for dry air.

Surface pressure is computed from

@p

@ t = m2 

@pu=m

@x +

@pv=m

@y  

@p_

@ ; (2:1:4)

which is used in its vertically integrated form

@p

@ t = m

2 Z 1

0



@pu=m

@x +

@pv=m

@y 

d: (2:1:5)

Then the vertical velocity in -coordinates, _ , is computed from (2.1.4) by vertical

integration. Thus

_ = 

1

p

Z 

0



@p

@ t + m

2 

@pu=m

@x +

@pv=m

@y 

d0 ; (2:1:6)

where 0 is a dummy variable of integration and _( = 0) = 0.

In the thermodynamic equation, (2.1.3), ! = dp

dt and is calculated from

! = p

_ + 

dp

dt ; (2:1:7)

where

dp

dt =

@p

@ t + m



u

@p

@x + v

@p

@y 

: (2:1:8)

The hydrostatic equation is used to compute the geopotential heights from the virtual

temperature, Tv :

@

@ln( + pt=p) = RTv 

1 +

qc + qr

1 + qv

1

; (2:1:9)

2
where Tv is given by Tv = T (1 + 0:608qv), and qc and qr are the mixing ratios of cloud

water and rain water.

2.2 Nonhydrostatic model equations

For the nonhydrostatic model we dene a constant reference state and perturbations

from it, as follows:

p(x; y; z; t) = p0(z) + p

0(x; y; z; t);

T (x; y; z; t) = T0(z) + T 0(x; y; z; t);

(x; y; z; t) = 0(z) + 

0(x; y; z; t):

Typically the temperature prole for the reference state may be an analytic function that

ts the mean tropospheric temperature prole.

The vertical -coordinate is then dened entirely from the reference pressure.

 =

p0  pt

ps  pt

;

where ps and pt are the surface and top pressures respectively of the reference state and

are independent of time. The total pressure at a grid point is therefore given by

p = p

 + pt + p0 ;

where p(x; y) = ps(x; y)  pt . The three-dimensional pressure perturbation, p0 , is a

predicted quantity.

The model equations (Dudhia 1993) are then given by the following:

Horizontal momentum;

@pu

@ t =  m2 

@puu=m

@x +

@pvu=m

@y  

@pu_

@ + uDIV

 mp





@p0

@x 



p @p@x

@p0

@  + p

f v + Du (2:2:1)

@pv

@ t =  m2 

@puv=m

@x +

@pvv=m

@y  

@pv_

@ + vDIV

 mp





@p0

@y 



p @p@y

@p0

@   p

f u + Dv (2:2:2)

3
Vertical momentum;

@pw

@ t =  m

2 

@puw=m

@x +

@pvw=m

@y  

@pw_

@ + wDIV

+ p

g

0





1

p @p0@ +

T 0vT T0p0

T p0

  p

g [(qc + qr )] + Dw (2:2:3)

Pressure;

@pp0

@ t =  m

2 

@pup0=m

@x +

@pvp0=m

@y  

@pp0_

@ + p

0DIV

 m

2

p


p 

@u=m

@x 



mp @p@x

@u

@ +

@v=m

@y 



mp @p@y

@v

@ 

+ 0g
p

@w

@ + p

0gw (2:2:4)

Temperature;

@pT

@ t =  m

2 

@puT =m

@x +

@pvT =m

@y  

@pT _

@ + T DIV

+

1

cp



p

 Dp0

Dt  0gpw  Dp0 + p

 Q_

cp

+ DT ; (2:2:5)

where

DIV = m2 

@pu=m

@x +

@pv=m

@y  +

@p_

@ ; (2:2:6)

and

_ = 

0g

p w  mp @p@x

u  m

p @p@y

v: (2:2:7)

The DIV terms are not in the hydrostatic equations and arise because pis now

constant in time. Thus the hydrostatic continuity equation no longer applies, leaving

the right hand side terms in (2.2.6) uncancelled by the surface pressure tendency. The

equations are thus in advective form.

Equation (2.2.4) can be derived from the fully compressible mass continuity relation

and the perfect gas law. The only term neglected in equations (2.2.1)-(2.2.5) is a diabatic

term contributing to the perturbation pressure tendency in (2.2.4). This term is negligible

in normal meteorological regimes since it only forces a small divergence (i.e. expansion)

4
in regions of heating. The Dp0 term in (2.2.5) is a small correction to DT allowing for

horizontal pressure dierences in thermal diusion.

2.2.1 Complete Coriolis force option

In the nonhydrostatic model it is possible to include the other components of

the Coriolis force that are neglected in the traditional approximation. The full

Coriolis force leads to a small upward/downward acceleration on westerly/easterly 
ows

and an westward/eastward acceleration on upward/downward 
ows in addition to the

rightward/leftward de
ection of horizontal 
ows in the northern/southern hemisphere.

To determine this force, two additional parameters are dened. We will refer to the

other component of the Coriolis parameter as e = 2
 cos , where 
 is the angular velocity

of the earth and  is the latitude. The other new parameter is , which is the angular

dierence between the y-axis of the grid and true north. It is found from

tan  = cos 

@=@y

@=@y ; (2:2:8)

where  is longitude. A special provision is made for the dateline. Thus  is positive

if north is rotated clockwise from the y-axis.

The momentum equations are then given by the following:

Horizontal momentum

@pu

@ t =  m2 

@puu=m

@x +

@pvu=m

@y  

@pu_

@ + uDIV

 mp





@p0

@x 



p @p@x

@p0

@  + p

f v  pew cos  + Du (2:2:9)

@pv

@ t =  m

2 

@puv=m

@x +

@pvv=m

@y  

@pv_

@ + vDIV

 mp





@p0

@y 



p @p@y

@p0

@   p

f u + pew sin  + Dv (2:2:10)

Vertical momentum;

@pw

@ t =  m

2 

@puw=m

@x +

@pvw=m

@y  

@pw_

@ + wDIV

+ p

g

0





1

p @p0@ +

T 0vT T0p0

T p0

  p

g [(qc + qr )]

+ p

e(u cos   v sin ) + Dw: (2:2:11)

5
2.3 Nonhydrostatic Finite Dierence Algorithms

The B-grid staggering of horizontal velocity variables with respect to the other elds

is shown in Fig. 2.1. Vertical velocity is staggered vertically. Noting that the j index

increments in the x direction, and i in the y direction, the conventional notation will be

as follows.

ax = (ai;j+ 1

2  ai;j 12 )=x: (2:3:1)

a

x =

1

2

(ai;j+ 1

2 + ai;j 12 ); (2:3:2a)

Multiple averaging terms such as axyy can also be dened as successive averages where the

order of superscripts does not matter, e.g.,

a

xyy = axy y

:

Averaging vertically allows for non-uniform grid-lengths and nonlinearly varying elds,

such as temperature and water vapor, by suitably weighting the values.

Thus for half-level elds averaged to full levels

a

 =

ak+ 1

2 (k  k 12 ) + ak 12 (k+ 12  k )

(k+ 1

2  k 12 )

; (2:3:2b)

while averaging full-level elds to half levels uses an equation similar to (2.3.2a). For

temperature, a is the potential temperature, and for water vapor, a is log qv.

The spatial dierencing of the terms in the horizontal momentum prediction equations

is [including the map-scale factor m(x; y)],

@pdu

@ t =  m

2 2

4 

ux pdu

m

xyy!

x +

 

uy pdv

m

xyx!

y

3

5  (pdu

_ xy)

+ uDIV xy  mpd

xy "

p0x

y  (p

)x

y p0

p

xy#

+ p

df v + D(pdu); (2:3:3)

6
y(I)

x(J)

I+1,J-1

I ,J-1

I-1,J-1

I+1, J I+1,J+1

I,J I ,J+1

I-1,J I-1,J+1

I,J-1 I,J

I-1,J-1 I-1,J

(U,V)

Fig. 2.1 Horizontal grid structure in the model.

p

*

( T q Φ σ ̇ ω qc qr ,,, ,, , , , R )

7
@pdv

@ t =  m

2 2

4 

v

x pdu

m

xyy!

x +

 

v

y pdv

m

xyx!

y

3

5  (pdv

_ xy )

+ vDIV xy  mpd



xy "

p0y x  (p)y x p0

p

xy#

 p

df u + D(pdv); (2:3:4)

where pd = pxy, and DIV , the mass divergence term, is given by

DIV = m

2 2

4 

pdu

m

y!

x +

 

pdv

m

x!

y

3

5 + p

_  : (2:3:5)

The triple averaging in the horizontal momentum advection terms follows that of the

hydrostatic model as discussed by Anthes (1972). The subgrid-scale and diusion operators

are represented by D(a) = Khx2(axxxx + ayyyy )+(Kvaz )z+ (PBL tendencies), where

the fourth-order scheme is modied to second-order near the boundaries.

The coordinate vertical velocity, _, is obtained from

_ = 

0g

p w  mp px

xuxy  m

p py

yv

xy

; (2:3:6)

and the vertical momentum equation is

@pw

@ t =  m

2 "

wx pu

m

y

x +



wy pv

m

x

y

#

 (pw

_ 

)

+ wDIV  + p

g

0



 "

1

p p

0 

1




p0T0

p0T

#

+ p

g

0



 "

T 0vT

 

R

cp

p0T0

p0T

#

 p

g(qc + qr ) + D(p

w): (2:3:7)

The pressure tendency equation, neglecting diabatic terms, is given by

@pp0

@ t =  m

2 "

p0 x pu

m

y

x +



p0 y pv

m

x

y

#

 (pp0 _ )

+ p

0DIV + p0gw  m2p


p

u

m

y

x  (px

)x 1

mp u

xy

8
+



v

m

x

y  (p y

)y

1

mp v

xy 

0g

m2p w

; (2:3:8)

and temperature tendency is dierenced as

@pT

@ t =  m

2 "

T

x pu

m

y

x +



T

y pv

m

x

y

#

 (pT

_ )

+ T DIV +

1

cp



p

 Dp0

Dt  0gpw

  D(p

p0)

+ p

 Q_

cp

+ D(p

T ); (2:3:9)

where Dp0=Dt is dierenced like the corresponding terms in (2.3.8). Moisture variables

have similar advection forms to those in (2.3.8) and (2.3.9) except when using the upstream

option where qx

is replaced by the upstream value alone.

9
2.4 Hydrostatic Finite Dierence Algorithms

The hydrostatic nite dierencing of advection, Coriolis and heating follows (2.3.3),

(2.3.4) and (2.3.9) without the DIV terms. The pressure gradient terms in (2.3.3) become

P G = mpdx

y  mRTv xy

(1 + pt=pd) px

y

; (2:4:1)

and likewise for the y-gradient in (2.3.4). The surface pressure tendency is found from the

integration over all (KMAX) layers of thickness (k),

@p

@ t = m2

KMAX

X

k=1

2

4 

pdu

m

y!

x +

 

pdv

m

x!

y

3

5 (k): (2:4:2)

Then _ is found from downward integration,

_ (k + 1) = _ (k) 

@p

@ t

(k)

p  m

2 2

4 

pdu

m

y!

x +

 

pdv

m

x!

y

3

5

(k)

p

; (2:4:3)

using the upper boundary condition that _(k = 1) = 0. The adiabatic term in (2.3.9),

represented by the second set of terms in square brackets, becomes p! in the hydrostatic

model, where ! is dened by

! =

dp

dt = p

_  +  

@p

@ t + muxypx

x + mv

xyp

y

y

: (2:4:4)

The integration of the hydrostatic equation to obtain geopotential height, , in the

hydrostatic model is done as follows.

 =  RTvL



ln( + pt=p); (2:4:5)

where

L =



1 +

qc + qr

1 + qv

1

;

and allows for water loading when the explicit moisture scheme is used. Because  is

required on the velocity levels (half-levels), it has to be integrated rst between the surface,

where  = 1 and  = gh (h is the terrain height above sea-level), and the lowest half-level

using (2.4.5) with just the lowest-level values Tv ; qv ; qc; qr . At all other levels (2.4.5) uses

vertical averaging between two levels.

10
The temporal dierencing in the hydrostatic and nonhydrostatic models consists of

leapfrog steps with an Asselin lter. With this time lter, splitting of the solution often

associated with the leapfrog scheme is avoided. It is applied to all variables as

^

t = (1  2)t + (t+1 + ^t1

); (2:4:6)

where ^ is the ltered variable. The coecient  in the model is 0.1 for all variables. For

stability, diusion terms are evaluated on the variables at time t  1, as are the terms

associated with the moist physical processes.

2.5 Time splitting

In both the nonhydrostatic as well as the hydrostatic numerics, a time splitting scheme

is applied to increase eciency. Because the nonhydrostatic equations above are fully

compressible, they permit sound waves. These are fast and require a short time step for

numerical stability. For the hydrostatic equations, fast moving external gravity waves are

the limiting factor. The techniques described next are designed to split these fast moving

waves from the rest of the solution.

2.5.1 The nonhydrostatic semi-implicit scheme

For the nonhydrostatic equations it is possible to separate terms directly involved with

acoustic waves from comparatively slowly varying terms, and to handle the former with

shorter time steps while updating the slow terms less frequently. The reduced equation

set for the short time step makes the model more ecient. The separated equations only

contain interactions between momentum and pressure and can be written as:

Horizontal momentum;

@u

@ t +

m





@p0

@x 



p @p@x

@p0

@  = Su (2:5:1:1)

@v

@ t +

m





@p0

@y 



p @p@y

@p0

@  = Sv (2:5:1:2)

Vertical momentum;

@w

@ t 

0



g

p @p0@ +

g




p0

p = Sw (2:5:1:3)

Pressure;

@p0

@ t + m

2


p 

@u=m

@x 



mp @p@x

@u

@ +

@v=m

@y 



mp @p@y

@v

@ 

11
0g
p

p @w@  0gw = Sp0 ; (2:5:1:4)

where the S terms contain advection, diusion, buoyancy and Coriolis tendencies. These

are kept constant during the sub-steps. Note that only part of the p0=p term is in (2.5.1.3),

where the rest has been absorbed in the buoyancy term that contributes to Sw.

The method of solution follows the semi-implicit scheme of Klemp and Wilhelmson

(1978) for the short time step. Starting with u; v; w; p0 known at time  , rst the two

horizontal momentum equations are stepped forward to give u+1 and v+1 which are then

used in the pressure equation, giving a time-centered explicit treatment of horizontally

propagating sound waves. Vertical propagation of sound waves is treated implicitly by

making w+1 and p0+1 depend upon time-averaged values of p0 and w respectively in

(2.5.1.3) and (2.5.1.4). For instance, where p0 appears in (2.5.1.3) it is represented by

p0 =

1

2

(1 + )p

0+1 +

1

2

(1  )p

0

;

and similarly for w in (2.5.1.4). The parameter  determines the time-weighting, where

zero gives a time-centered average and positive values give a bias towards the future time

step that can be used for acoustic damping. In practice, values of  = 0:2  0:4 are used.

With second-order vertical spatial derivatives the nite dierence forms of equations

(2.5.1.3) and (2.5.1.4) can be combined, eliminating p0+1, into a nite dierence equation

for w+1, which is solvable by direct recursion on a tri-diagonal matrix.

The implicit vertical dierencing scheme allows the short time step to be independent

of the vertical resolution of the model, which is important for eciency, and thus the

step only depends upon the horizontal grid length. Additionally, the divergence damping

technique of Skamarock and Klemp (1992) is used to control horizontally propagating

sound waves. This method is similar to using time-extrapolated pressure terms in (2.5.1.1)

and (2.5.1.2), where in practice the extrapolation is about 0.1  .

Temperature and moisture are predicted using the normal leapfrog step, t, because

they have no high-frequency terms contributing to acoustic waves. The slow terms for

momentum and pressure contained in the S-terms above are also evaluated on these

leapfrog steps, but for these variables the march from tt to t+t is split into typically

four steps of length  during which momentum and pressure are continually updated.

12
2.5.2 The hydrostatic split-explicit scheme

When numerically solving the hydrostatic equations of motion, the stability criterion

is severely limited by external gravity waves. These are very fast moving gravity waves

that are small in amplitude (quasi-linear) and contain only a small fraction of the total

energy. Hence they change slowly over the time scale of the Rossby waves. Because of this,

splitting methods have been developed to split these fast waves from the solution (similar

also to the above method for the nonhydrostatic equations to split sound-waves). From

all the existing dierent options, we have chosen a method developed by Madala (1981).

This scheme separates the terms governing the gravity modes from those governing the

Rossby modes. The term \split" here refers to the separation of the motion in terms of

eigenmodes. Similar to the nonhydrostatic method, the equations are rewritten in nite

dierence form as

@Psu

@ t + x = Au; (2:5:2:1)

@Psv

@ t + y = Av ; (2:5:2:2)

@PsT

@ t + M2  D = AT ; (2:5:2:3)

@Ps

@ t + N1  D = 0; and (2:5:2:4)

 = M1  T : (2:5:2:5)

where the right hand sides change slowly over the time scale of the Rossby-waves. Matrices

M1, M2, and vector N1 are independent of x, y, and t. Notice the similarity to the

nonhydrostatic splitting method (equations 2.5.1-2.5.4). However, rather then integrating

the \fast" terms on a small time-step directly, the method described below only computes

correction terms to the equations, making this process extremely ecient. To illustrate

this, we follow Madala (1981). From the governing equations he derives equations for the

mass divergence D and the generalized geopotential . They are

@D

@ t + [

2x + 2

y ] = xAu + yAv (2:5:2:6)

and

@

@ t + M3  D = M1  AT : (2:5:2:7)

13
Integrating equations (2.5.2.1-2.5.2.3) from t  t to t + t, where t is the time step of

the slow Rossby modes, one gets

psu(t + t)  psu(t  t) + 2tx = 2 ~tAu(t); (2:5:2:8)

psv(t + t)  psv(t  t) + 2ty = 2 ~tAv(t); (2:5:2:9)

psT (t + t)  psT (t  t) + 2tM2 = 2 ~ tAT (t); (2:5:2:10)

where the operator (~) for the split-explicit scheme is dened as

~ = 

tXm

n=1

(t  t + n );

where m =  t

. Denoting with superscript ex solutions computed using only the explicit

time integration over 2t, equations (2.5.2.8-2.5.2.10) can be written as

psu(t + t) + 2tx[~  (t)] = psuex(t + t); (2:5:2:11)

psv(t + t) + 2tx[~  (t)] = psv

ex(t + t); (2:5:2:12)

psT (t + t) + 2tM2[D~  D(t)] = psT ex(t + t): (2:5:2:13)

Here (t) and D(t) have been computed using the explicit time integration over 2t.

Similar, for the pressure tendency we can write

Ps(t + t) + 2tN1  [D~  D(t)] = P ex(t + t): (2:5:2:14)

To nd equations for the correction terms on the left hand side of equations (2.5.2.11-

2.5.2.13), the divergence and geopotential equations (2.5.2.6-2.5.2.7) are then solved over

the the small time-steps using

[D(t + (n + 1) )  D(t)]  [D(t + (n  1) )  D(t)]

+ 2 (

2x + 2

y )[(t + n )  (t)]

=

1

mi

[Dex(t + t)  D(t  t)]

(2:5:2:15)

14
and

[(t + (n + 1) )  (t)]  [(t + (n  1) )  (t)]

+ 2M3[D(t + n )  D(t)]

=

1

mi

[ex(t + t)  (t  t)]

: (2:5:2:16)

The correction terms themselves are integrated in equations (2.5.2.15)and (2.5.2.16), and

then added to equations (2.5.2.11-2.5.2.14).

 , the timestep of the fast modes, of course varies with the mode. For a clean

separation of the modes, a vertical normal mode initialization developed and applied to the

MM4/MM5 system by Errico (1986) is used at the beginning of the model run to calculate

the vertical modes. In MM5, only the external and the fastest internal mode are being

considered with dierent time steps. This allows the time-steps of the slow tendencies to

be twice as large as they were with the previously used Brown-Campana (1978) algorithm,

and they are comparable to the ones used in the nonhydrostatic numerics.

15
2.6 Lateral Boundary conditions for the coarsest mesh domain

2.6.1 Sponge Boundary Conditions

The sponge boundary condition is given by



@

@ t 

n = w(n) 

@

@ t 

MC

+ (1  w(n)) 

@

@ t 

LS

; (2:6:1)

where n = 1; 2; 3; 4 for cross-point variables, n = 1; 2; 3; 4; 5 for dot-point variables, 

represents any variable, MC denotes the model calculated tendency, LS the large-scale

tendency which is obtained either from observations or large-scale model simulations (oneway nesting), and n is the displacement in grid-points from the nearest boundary (n = 1

on the boundary). The weighting coecients w(n) for cross point variables (counting from

the boundary points inward) are 0.0, 0.4, 0.7, and 0.9, while for dot-point variables they are

equal to 0.0, 0.2, 0.55, 0.8, and 0.95. All other points in the coarse domain have w(n) = 1.

The above method cannot be used for the nonhydrostatic part of the model.

2.6.2 Nudging Boundary Conditions

The relaxation boundary condition involves \relaxing" or \nudging" the modelpredicted variables toward a large-scale analysis. The method includes Newtonian and

diusion terms



@

@ t 

n = F (n)F1(LS  MC)  F (n)F22(LS  MC): n = 2; 3; 4 (2:6:2)

F decreases linearly from the lateral boundary, such that

F (n) = 

5  n

3 

n = 2; 3; 4; (2:6:3)

F (n)=0 n > 4; (2:6:4);

where F1 and F2 are given by

F1 =

1

10t (2:6:5)

and

F2 = s

2

50t

: (2:6:6)

This method is also used for the nonhydrostatic part of the model to nudge the

pressure perturbation to the observations or larger-scale model simulations. However, for

16
the nonhydrostatic solver the vertical velocity is not nudged. It can vary freely, except

for the outermost rows and columns, where zero gradient conditions are specied. For the

velocity components, the values at the in
ow points are specied in a manner similar to

the specication of temperature and pressure. The values at the out
ow boundaries are

obtained by extrapolation from the interior points. These boundary values are required

only in the computation of the nonlinear horizontal momentum 
ux divergence terms;

They are not required in the computation of the horizontal divergence.

2.6.3 Moisture variables

Cloud water, rain water, snow, and ice are considered zero on in
ow and zero gradient

on out
ow. There is an option to specify the boundary values in the same way as for the

other variables (e.g., these variables may be known in a one-way nesting application).

2.7 Upper radiative boundary condition

An option in the nonhydrostatic model is the upper radiative boundary condition.

Klemp and Durran (1983) and Bougeault (1983) have developed an upper boundary

condition that allows wave energy to pass through unre
ected. It can be expressed for

hydrostatic waves as

p^ =

N

K

w; ^ (2:7:1)

where p^ and w^ are horizontal Fourier components of pressure and vertical velocity

respectively,  and N are the density and buoyancy frequency near the model top, and

K is the total horizontal wavenumber of the Fourier component. This expression should

be enforced for all components if the energy transport is to be purely upward with no

re
ection.

The upper boundary condition is combined with the implicit pressure/vertical

momentum calculation. Before either value at time n + 1 is known, the values at the

top model level (w1 is staggered half a grid length above p1) can be expressed as

p

n+1

1 = b + awn+11 ; (2:7:2)

where the coecient a(x; y; t) is dependent upon the thermodynamic structure and the

bottom boundary condition on w in the model column. It varies within only 5 per cent

of a constant value even with high terrain, and is also not strongly time-dependent. The

17
value of b(x; y; t) depends on pressure and most of the pressure tendency terms, and both

a and b are known at this stage. So transforming, assuming a varies little about a non-zero

constant and taking a mean value a

p^ = ^b + aw: ^ (2:7:3)

Combining (2.7.3) with the radiative condition (2.7.1) for wavenumber K = 2=, taking

N at the top of the model, and eliminating ^p, gives

w^ = K^

b

N  aK

: (2:7:4)

Using a limited-area 2D cosine transform, the forward transform, multiplication and

backward transform can be combined into a single operator on the b eld to give wn+1

1 .

Hence

wIJ = X

I+6

i=I6

J

X

+6

j=J6

ij bij ; (2:7:5)

where we have localized the transform to 13 13 points, and array  can be precalculated

and kept constant for the time integration. The elements of  are found from

ij = X

6

k=0

X

6

l=0

ij k l

36

cos

2ki

12

cos

2lj

12

f (K);

(2:7:6)

with f (K) = K

NaK

and K = (k^2 + ^l2)

1

2 .  = 1 except for limits of summations where

 = 1

2

.

Following the suggestion of Klemp and Durran, the nite dierencing of pressure

gradients and divergences should be taken into account in dening the eective

wavenumbers. For a B-grid staggering, the eective wavenumbers can be expressed in

terms of the dimensionless wavenumbers, k and l, where

k^ =

2

x

sin k

12

cos

l

12

; (2:7:7a)

^

l =

2

x

sin

l

12

cos

k

12

; (2:7:7b)

and x is the grid length.

The scheme is summarized as follows; by the precalculation of parameters a and N

for the model domain, use of (2.7.6) to precalculate coecients , then implementation of

(2.7.5) during the simulation.

18
3. The Mesh renement scheme

The 2-way interactive mesh renement scheme is constructed to allow for an arbitrary

number of overlapping and translating rectangular grids with an arbitrary number of

renement levels. The grids must be aligned with the model coordinates (no rotating

meshes), and the mesh renement ratio of the temporal and spatial grid increments is

common for all meshes, and currently set to three. Vital parts of this renement scheme

are the interpolation routines (Smolarkiewicz and Grell, 1992), which are used upon

initialization of new nests as well as for dening the boundaries of the ne meshes. If

the user can supply his own analysis for the ner grids (or a ner grid), the interpolated

elds can be overwritten. In the following section we describe the heart of the scheme, the

monotone interpolation routines.

3.1 The monotone interpolation routines

The most vital element of any mesh renement scheme is an accurate and ecient

interpolation procedure. A complaint about traditional polynomial-tting methods used

for interpolating scalar elds dened on a discrete mesh is that they often lead to spurious

numerical oscillations in regions of steep gradients of the interpolated variables. In order

to suppress computational noise, which is characteristic of quadratic and higher-order

interpolation schemes, one often implements a smoothing procedure or increased diusion.

These, however, introduce excessive numerical diusion that smears out sharp features

of interpolated elds. A more advanced approach invokes the so-called shape-preserving

interpolation, which incorporates appropriate constraints on the derivative estimates used

in the interpolation schemes (see Rasch and Williamson (1990) for a review). In MM5

we consider as an alternate approach a class of schemes derived from monotone advection

algorithms (Smolarkiewicz and Grell, 1992). Smolarkiewicz and Grell (1992) show that

the interpolation problem becomes identical to the advection problem, when the distance

vector is replaced by the velocity vector (see also the end of this section). Here we will

describe the implementation of the advection algorithm used in MM5. The interested

reader is referred to Smolarkiewicz and Grell (1992) for a detailed derivation of the

\advection-interpolation" equivalence.

Since shape preservation and monotonicity are important in the interpolation problem,

19
we chose the Flux Corrected Transport (FCT) scheme that uses the high-order accurate

constant-grid-
ux dissipative algorithms developed by Tremback et al. (1987). We will

rst describe, in abbreviated form, a general FCT algorithm, as used in MM5. Given

the exactness of the alternate-direction representation of the interpolation algorithm,

it is sucient to consider only one-dimensional FCT schemes. Starting with the onedimensional advection equation in 
ux-form

@

@ t = 

@u

@x ; (3:1:1)

where  is a scalar variable advected with a 
ow eld u(x; t), an FCT advection scheme

may be compactly written as

n+1

i = n+1i  (A~

i+1=2  A~i1=2); (3:1:2)

where  denotes a low-order, monotone solution to (3.1.1), and A~ is the antidiusive 
ux,

limited such as to ensure that the solution (3.1.2) is free of local extrema absent in the

low-order solution. Note that

A~

i+1=2 = min 

1; #

i

; "

i+1

Ai+1=2+ + min 

1; "

i

; #

i+1

Ai+1=2

; (3:1:3)

where

Ai+1=2  F Hi+1=2  F Li+1=2; (3:1:4)

with F H and F L denoting 
uxes from a high-order and a low-order advection scheme,

respectively.[ ]+  max(0; ) and [ ]  min(0; ) are the positive- and the negative-part

operators, respectively, and



"

i 

MAX

i  n+1

i

AIN

i + "

; 

#

i 

n+1

i  MIN

i

AOUT

i + "

; (3:1:5a; b)

where AIN

i

, AOUT

i are the absolute values of the total incoming and outgoing A-
uxes,

(3.1.4), from the i-th grid box, respectively. " is a small value, e.g.  1015, and allows

for ecient coding of -ratios when AIN

i or AOUTi vanish. The limiters MAX

i and MINi

dene monotonicity of the scheme (i.e., by design MIN

i  n+1  MAXi ), and they are

traditionally specied (Zalesak 1979) as

MAX

i = max ni1; ni

; n

i+1; n+1

i1 ; n+1i

; n+1

i1 

(3:1:6a)

20
MIN

i = min ni1; ni

; n

i+1; n+1

i1 ; n+1i

; n+1

i1 

: (3:1:6b)

A shape-preserving interpolation scheme requires less restrictive monotonicity

constraints than a conservative advection scheme. The minima over  ratios appearing

in (3.1.3) ensure that the antidiusive 
ux attributed to the i + 1=2 position on the grid

does not contribute to the generation of spurious extrema, either in gridbox i or in gridbox

i + 1. However, monotonicity of the interpolation scheme only requires that n+1

i = (xo )

is free of spurious extrema. Consequently, equation (3.1.3) may be replaced by

A~

i+1=2 = min 

1; #

i 

Ai+1=2+ + min 

1; "

i 

Ai+1=2

: (3:1:30)

Furthermore, since the eective 
ow eld is constant, and therefore incompressible, the

limiters in (3.1.6) may be simplied to

MAX

i = max ni1; ni

; n

i+1; n+1

i

; ; MIN

i = min ni1; ni

; n

i+1; n+1

i 

; (3:1:60a; b)

where the redundant dependence of the limiters on n+1

i has been retained to ensure

strictly nonnegative values of the  ratios in (3.1.5) (cf., Section 3.1 in Smolarkiewicz and

Grabowski, 1990). Since the low-order solution may always be written as an old value,

minus the divergence of 
uxes from the low-order scheme, the entire algorithm consisting

of (3.1.2), (3.1.30), (3.1.4), (3.1.5), and (3.1.60) is in the form of a general advection scheme.

The advection schemes used to calculate the high- and low-order 
uxes for the above

equations are from Tremback et al. (1987). They derive as follows. Starting with the 
ux

form of the one-dimensional advection equation (3.1.1) in nite dierence form

n+1

i = + n

i +

t

x

[Fi+1=2  Fi1=2]

= n

i +

x

t

"X

m

bmn

i+1+m X

m

bmn

i+m

#

; (3:1:7)

where

Fi+1=2 = X

m

bmn

i+1+m (3:1:8)

and

Fi1=2 = X

m

bmn

i+m (3:1:9)

21
were used. Following Tremback et al. (1987), the solutions for the even-order schemes

which are used in the mesh renement scheme are then given by

Fi+1=2

t

x

= +



2

(i  i+1)+

+

2

2

2(i + i+1)

(3:1:10)

for second order accuracy;

Fi+1=2

t

x

= +



12

(i1  7i  7i+1 + i+2)

+

2

24

(i1  15i + 15i+1  i+2)

+

3

12

(i1 + i + i+1  i+2)

+

4

24

(i1 + 3i  3i+1 + i+2)

(3:1:11)

and for fourth order accuracy;

Fi+1=2

t

x

= +



60

(i2 + 8i1  37i  37i+1 + 8i+2  i+3)

+

2

360

(2i2 + 25i1  245i + 245i+1  25i+2 + 2i+3)

+

3

48

(i2  7i1 + 6i + 6i+1  7i+2 + i+3)

+

4

144

(i2  11i1 + 28i  28i+1 + 11i+2  i+3)

+

5

240

(i2 + 3i1  2i  2i+1 + 3i+2  i+3)

+

6

720 (i2 + 5i1  10i + 10i+1  5i+2 + i+3)

; (3:1:12)

for sixth order accuracy;  is dened as

 = U

t

x

: (3:1:13)

In MM5, equations (3.1.10 - 3.1.12) are used together with (3.1.1), (3:1:30), (3.1.4), 3.1.5),

and (3:1:60) to solve the interpolation problem. Note that the velocity vector, is replaced by

the distance vector, Xd, which, with a mesh-renement ratio of three, simply becomes 1=3

or 2=3. For interpolating boundary conditions to the ner meshes, fourth order accuracy

22
is used, while for new nest initialization, sixth order accuracy is used. While the new

nest initialization covers the whole domain, boundary interpolation is performed for the

outermost 2 rows and columns of the nest. Two rows were necessary to ensure that the

same operators were applied to each nested grid-point (including fourth-order diusion).

3.2 Overlapping and moving grids

The mesh-renement scheme allows for overlapping grids on the same nest-level. To

ensure numerical stability, the solution in the overlap region has to be identical. To

accomplish this, after each time-step of the grids in question, the boundary conditions

in the overlap regions are provided by the overlapping mesh. It is very important that this

procedure be performed at every timestep.

Nests can also be moved at any time in the forecast. This can be done many times,

and for any distance (integer number of grid points). The user may also move the nests

automatically if he supplies an algorithm to do so. Upon a move, a new nest initialization

is performed rst. Then all high-resolution information from the previous location of the

mesh is used to overwrite the elds of the newly initialized mesh. Therefore, to ensure

best use of high resolution information, it is better to move a nest more often and for a

smaller distance.

3.3 The feedback to the coarser grids

Since the mesh renement ratio in MM5 is set to three, a higher resolution mesh has

to be integrated three times as often as its \Mother Domain"(MD), where MD means the

coarser domain from which it gets its boundary conditions. To keep the solutions in a 2-way

interaction from diverging, the meteorological elds have to be fed back from the higherresolution mesh to its MD. This is done at the end of the three time-step integration.

Naturally, when this is done without smoothing or averaging, the solution on the MD

will appear somewhat noisy (diluted with small-scale information). To avoid numerical

instability, the following methods are supported in MM5 to remove non-resolvable noise

from the MD. Note that these smoothers are only applied over an interior area that is

completely determined by the higher resolution domain. It is important that input into

the nest, and feedback back to the MD does not overlap. The smoother that is used by

23
the MM5 system in various forms was dened by Shapiro (1970) as

(i; j) =(i; j)

+



2

(1  )((i + 1; j) + (i  1; j) + (i; j + 1) + (i; j  1) + 4(i; j))

+

2

4

((i + 1; j + 1) + (i + 1; j  1) + (i  1; j + 1) + (i  1; j  1)  4(i; j))

(3:3:1)

3.3.1 A Nine-point averager

This method was in the original MM4 nested version (Zhang et al. 1986). It is a

feedback method that averages information for a whole MD grid box (surrounding and

centering on the nested grid point). However, it does not take out all non-resolvable

information on the MD. It also imposes a severe restriction on the terrain for the hydrostatic

model. In case of overlapping and moving nests on several nest levels, it is very elaborate

and complicated to apply. It is still an option in the model, because it may be useful

for simpler applications (like one coarse and one nested domain). However, care must be

taken to create a terrain data set that is consistent with this method. The operator that

is applied to the nested grid-points (note that nothing is done to the MD) is dened by

using  = 0:5 in (3.3.1).

3.3.2 A Smoother-Desmoother

The smoother-desmoother is a lter that removes 2x waves and damps short waves,

but leaves long waves almost unaected. It is much more selective than diusive smoothers.

It is applied to the \coarser grid" only in the area where the coarse grid values are

overwritten with the nested grid values.

A single pass of the smoother-desmoother involves two steps. Equation (3.3.1) is

used rst to smooth the elds, then to desmooth the elds. 1 = 0:50 is used for the

smoothing coecient, and 2 = 0:52 for the desmoothing coecient. The rst step

strongly smoothes the eld, completely removing the 2x wave, and the second step

attempts to restore the other waves to their original amplitudes. There are two passes of

the smoother-desmoother applied in the model.

24
4. Four Dimensional Data Assimilation (FDDA)

The concept of combining current and past data in an explicit dynamical model

such that the model's prognostic equations provide time continuity and dynamic coupling

among the various elds has become known as four-dimensional data assimilation (FDDA).

Current interest in the use of FDDA in mesoscale models, for either model initialization

(dynamic initialization) or for use of the model as an analysis/research tool (dynamic

analysis), is a logical extension of the traditional link between ob jective analysis methods

and dynamic relationships.

Currently, two ma jor types of FDDA are used operationally and in research. The rst

is an intermittent process of initializing an explicit prediction model, using the subsequent

forecast (typically 3-12 h) as a rst guess in a static three-dimensional ob jective analysis

step, and then repeating the process for another forecast cycle. The second is a continuous

(every model time step) dynamical assimilation where forcing functions are added to the

governing model equations to gradually \nudge" the model state toward the observations.

This continuous nudging type of FDDA is used in the PSU/NCAR modeling system.

Nudging was rst developed and tested at Penn State by Kistler (1974), under Prof. J.

Hovermale, and by Anthes (1974), and Hoke and Anthes (1976). See Stauer and Seaman

(1990) for an historical overview of the technique.

Nudging or Newtonian relaxation is a relatively simple but very 
exible technique:

the data used for nudging can be of any type, measured or derived, analyzed to a grid for

assimilation into the model or inserted as individual observations. Gridded analyses of the

observations that are assimilated can be obtained by successive correction, variational, or

statistical optimal interpolation (OI) techniques, and the weights used when nudging to

individual observations can be simple Cressman-type (distance-weighted) functions or more

complicated statistical functions based on OI. It can be shown that successive corrections,

OI, and variational approaches are all fundamentally related to the \idealized analysis"

and thus to each other. In fact, nudging is basically a successive-correction technique which

uses a numerical model to include the time dimension.

The method of Newtonian relaxation or nudging relaxes the model state toward

the observed state by adding, to one or more of the prognostic equations, articial

25
tendency terms based on the dierence between the two states. The model solution can be

nudged toward either gridded analyses or individual observations during a period of time

surrounding the observations. These two techniques, hereafter referred to as \analysis

nudging" and \obs nudging", respectively, can be used individually or simultaneously on

any MM5 model grid. For guidance in selecting which nudging technique(s) to use for

your particular application, as well as suggested parameter specications, see Stauer and

Seaman (1990), Stauer et al. (1991) and Stauer and Seaman (1993).

4.1 Analysis Nudging

The analysis-nudging term for a given variable is proportional to the dierence between

the model simulation and an analysis of observations calculated at every grid point. The

general form for the predictive equation of variable (x; t) is written in 
ux form as

@p

@ t = F (; x; t)+G W(x; t)  (x)  p



( ^0 )+Gp Wp  p (x)  (^p



0 p



) (4:1:1)

All of the model's physical forcing terms (advection, Coriolis eects, etc.) are

represented by F , where  are the model's dependent variables, x are the independent

spatial variables, and t is time. The second and third terms on the right of (4.1.1) are

similar in form and represent the nudging terms for  and p*, respectively. Due to the 
ux

form of the predictive equation, the third term appears in (4.1.1) when nudging pressure in

the continuity equation of the hydrostatic version of MM5. (Note that this term is zero in

the nonhydrostatic version of MM5 because pis computed from the hydrostatic reference

state and is constant in time.)

With Gp = 0, or in the nonhydrostatic version of MM5, (4.1.1) simplies to

@p

@ t = F (; x; t) + G  W  (x)  p



( ^0  ) (4:1:2)

The nudging factor G determines the magnitude of the term relative to all the other

model processes in F . Its spatial and temporal variation is determined mostly by the

four-dimensional weighting function, W, which species the horizontal, vertical and time

weighting applied to the analysis, where W = wxywwt. The analysis quality factor, ,

which ranges between 0 and 1, is based on the quality and distribution of the data used to

26
produce the gridded analysis. The estimate of the observation for  analyzed to the grid

is ^0.

The nudging factor G is dened based on scaling arguments. Because the nudging

contribution is articial, it must not be a dominant term in the governing equations and

should be scaled by the slowest physical adjustment process in the model (inertial eects).

Thus the G is usually dened to be similar in magnitude to the Coriolis parameter, and

it must also satisfy the numerical stability criterion, G  1

t

. Typical values of G are

104s1to 103s1for meteorological systems, where values of G = 3 104s1to

6 104s1 are usually \large enough". A value of G which is too large will force the

model state too strongly toward the observations. This is undesirable because (a) the

ability of the model equations to resolve mass-momentum imbalances will be decreased;

and (b) the ability of the model to generate its own mesoscale meteorological structures

(e.g. fronts, squall lines) will be impaired by heavy insertion of the observed analyses.

Such problems arise because the analyses may not resolve these mesoscale structures or

may be contaminated by observational and analysis errors. On the other hand, if G is

too small, the observations will have minimal eect on the evolution of the model state,

allowing phase and amplitude errors to grow.

For simplicity, if we drop the physical forcing terms F from (4.1.2), and assume

W(x; t)=1; @p

@t = 0 and the observational analysis is perfect and time invariant, then

@

@ t = G( ^0  ) (4:1:3)

which has the solution

 = ^0 + (i  ^0)eGt (4:1:4)

where i is the initial value of  at the start of the nudging period. Therefore, the model

state approaches the observed state exponentially with an e-folding time of TG = 1

G

, which is about 0:93 h for G = 3 104s1. This implies that very high frequency


uctuations in the data, as might be available from wind prolers or Doppler radars (say,

every 5 min), would not be retained well unless G were much greater; but then the

nudging term may not be small compared to some terms of F.

27
Nudging the vorticity is a alternative method whereby the model's divergent wind is

allowed to freely respond in the model's geostrophic adjustment process. Equation 4.1.3

can be modied for u and v such that

@u

@ t = Gu  @@y (^

0   ) (4:1:5)

@v

@ t = Gv  @@x (^

0   ); (4:1:6)

where  is the model vorticity and ^0 is the analyzed observed vorticity. Letting a constant

G = Gu = Gv and forming the vorticity equation from (4.1.5) and (4.1.6) we get

@

@ t = Gr2

(^0   ) (4:1:7)

Therefore, the model vorticity is diused toward the observed vorticity. However, the

Laplacian in (4.1.7) introduces a scale dependence when nudging vorticity. The model

vorticity will be diused more strongly to small-scale features in the observed vorticity

analysis. Thus, when nudging toward large-scale vorticity elds, small-scale features in the

model vorticity will be selectively damped. After all factors are considered, it is generally

advisable to nudge the u- and v-components of the winds directly whenever possible (see

Stauer and Seaman, 1990).

Although this analysis-nudging technique has been traditionally used to assimilate 3-

D analyses based on rawinsonde observations, it can also be used to assimilate 2-D surface

analyses within the model PBL (Stauer et al., 1991).

While the twice daily rawinsonde observations are spaced at about 400 km and number

approximately 100 over the U.S., the spatial distribution of the surface data, available at

3-h intervals from the NCAR data archives, is considerably more dense. Depending on

the time of day, there are roughly 500- 1200 surface data sites with an average spacing

of about 100 km. The greater horizontal and temporal resolution of conventional surface

data is especially attractive for mesoscale data assimilation. The surface-level analyses are

assimilated in the same manner as the 3-D analyses, except for the vertical extent of their

in
uence. These surface analyses, which are assumed to be representative at 10 m above

ground level (AGL), should represent only those scales resolved by the model grid onto

28
which they are assimilated. This is necessary to avoid diculties related to the insertion

of small-scale divergence patterns (near 2x) which might interact adversely with the

model's parameterizations (e.g., the moisture-convergence criterion used to determine the

existence and intensity of Anthes-Kuo convection).

Eective assimilation of single-level data depends on the equivalent depth over which

the data are inserted into the model. Benecial eects on numerical forecasts can be

achieved by distributing single-level data throughout several model layers. This approach

requires that the data be applied in a consistent manner such that they are assured to

be representative of those layers. For example, the homogenizing eect of vertical mixing

during free convective conditions allows us to assume that surface-layer wind and mixing

ratio (q) observations can be applied throughout the model PBL according to a conceptual

model of boundary-layer structure. Such an idealized conceptual model is given by Garratt

et al. (1982) and is based on typical moderate to large instabilities observed at Wangara

and Minnesota (Fig. 4.1). However, in this same situation the frequent presence of

a shallow superadiabatic layer near the surface makes surface temperature or potential

temperature data poorly representative of the boundary layer as a whole, and similarity

relationships describing the potential temperature prole become inaccurate. The same is

true during nocturnal inversion conditions. These and other factors make assimilation of

single-level surface temperature observations unattractive for dening the temperature of

the mixed layer above the surface layer. For example, nudging towards an inaccurately

diagnosed mixed-layer temperature can cause serious errors in the simulated PBL depth or

even lead to a sudden spurious collapse of the unstable PBL structure. This can result from

assimilating a surface temperature observation which is cooler than the model-simulated

value by only a few tenths of one degree. In general, surface temperature data should not

be directly assimilated into the model (see Stauer et al., 1991).

Figure 4.1 shows the unstable lower troposphere comprised of three distinct layers:

a surface layer extending to height hs,awell-mixed layer from hs to h1 and a transition

layer extending from h1 to h2. With the x-axis dened parallel to the mean wind, Fig. 4.1

suggests that the v component of the wind is zero and there is thus no directional shear

through the lower two layers. The surface wind speed analysis, discussed above and

29
h2

h1

HEIGHT

v u θ

MEAN

WIND SPEED

MEAN POTENTIAL

TEMPERATURE

0

hs

Figure 4.1 Schematic representation of mean wind speed and potential temperature profiles in an

idealized conceptual model of the unstable atmospheric boundary layer.

-0.1 0.1 0.3 0.5 0.7 0.9 1.1

1.0

1.1

1.2

1.3

1.4

1.5

L = -100 m

L = -10 m

average

fit

u40/u10=1.000+0.320*z0**0.200:

u40/u10=1.169+0.315*z0:

0<z0<0.2

z0>0.2 _

z0 (m)

u40/u10

Figure 4.2 Schematic representing the relationship used to adapt the 10-m wind to the 40-m level

as a function of roughness length (z0) and Monin length (L).

30
assumed to apply at 10 m AGL, can be modied via similarity theory to apply at

hs. This modied wind for hs is representative of the mean mixed layer up to h1. These

assumptions allow the surface-layer wind information to be applied throughout the model's

multilayer PBL rather than at a single level. Similarly, the surface-layer mixing ratio (not

shown) is also assumed to be representative of the entire mixed layer, and may be applied

over several model layers during free-convective conditions.

The surface-layer wind analyses, however, must be adapted for the depth of the surface

layer, hs, which is assumed to be the height of the lowest model layer under unstable

conditions (40 m AGL for  = 0:995). Similarity relationships that assume a logarithmic

prole of wind with height are used to adjust the 10-m surface wind analysis to the lowest

model-layer height. Figure 4.2 shows that this wind adjustment is more strongly dependent

on roughness (z0) than on stability as measured by the Monin length (L); therefore, a \best-

t" relationship is determined for average stability conditions for roughness lengths ranging

from 0.0 to 1.0 m. Thus, during free-convective conditions, the mean wind assimilated

throughout the model PBL is based on the surface wind analysis modied to account for

the model surface-layer height and the grid-box roughness. The surface analysis of mixing

ratio is assumed to be representative for the lowest model layer and throughout the mixed

layer.

During unstable conditions, the Blackadar scheme may develop a PBL which extends

through several model layers. The \nudging correction" to the wind eld at the lowest

model layer is applied throughout the model PBL to simulate the conceptual model.

Because this idealization is also closely reproduced by the PBL scheme (without nudging),

the nudging strategy within the PBL is compatible with the boundary-layer physics. The

observed surface mixing ratio applied throughout the model mixed layer is adjusted if

necessary to remove any supersaturation with respect to the current model-simulated

temperature and moisture conditions for a given level and grid point. This prevents the

moisture assimilation term from becoming an articial source of precipitation. During

stable conditions, on the other hand, the height of the PBL is dened to be that of the

lowest model layer since the boundary layer is largely decoupled from the free troposphere

above.

31
Therefore, the 3-hourly surface-analysis nudging is also given by (4.1.2), but the

vertical extent of the nudging is controlled by the model- simulated PBL depth, with

^0 for wind and moisture adjusted as previously discussed above. The analysis condence

factor, , for the 3-h surface analyses, is functionally dependent on the spatial distribution

of the surface observations used to produce the analysis. Over land it varies from unity

at grid boxes within one-half the prescribed radius of in
uence of a surface observation to

0.2 for grid boxes outside the prescribed radius.

The vertical weighting factor, w, is dened as

w = w

R

 + w

S

  1: (4:1:8)

where wR

 and wS represent w for assimilation of 3-D rawinsonde and 2-D surface data,

respectively, and wS

 depends on the model-simulated PBL depth. The surface data are

assimilated with full strength (wS

 = 1:0) within the layers dening the PBL and with

reduced strength (wS

 = 0:9) one layer above (in the transition layer). The vertical

weighting function used to assimilate 3-D rawinsonde data is dened such that wR

 = 0:0 in

the PBL, 0.1 in the transition layer and 1.0 aloft. During stable conditions, therefore, the

surface data are applied with full strength only in the lowest model layer and with reduced

strength one layer above. Both types of analysis nudging generally assimilate temporally

interpolated gridded analyses; that is, ^0 in (4.1.2) is interpolated in time, for example,

from either 12-h 3-D analyses or 3-h 2-D surface analyses. Therefore, wt is usually set to

unity, except when decreasing the nudging at the end of a dynamic-initialization period.

4.2 Observational Nudging

This alternative scheme does not require gridded analyses of observations throughout

the case study period, and may be better suited for situations with high-frequency

asynoptic data (e.g., prolers), especially on the subalpha scales. Its form is similar to

(4.1.2) and it uses only those observations which fall within a predetermined time window

that is centered about each model time step. The set of dierences between the model and

the observed state is computed at the observation locations, and analyzed back to the grid

in a region surrounding the observations. The tendency for (x; t) with Gp = 0 is given

32
by

@p

@ t = F (; x; t) + G  p

 PN

i=1 W2

i (x; t)  
i  (o  ^)i

PN

i=1 Wi(x; t)

; (4:2:1)

where F and G are as dened earlier, subscript i denotes the i th observation of a total of

N which are within a preset radius of a given grid point, o is the locally observed value of

, and ^ is the model's prognostic variable interpolated to the observation location in three

dimensions. The observational quality factor,
 , which ranges from 0 to 1, accounts for

characteristic errors in measurement systems and representativeness. The four-dimensional

weighting function accounts for the spatial and temporal separation of the ith observation

from a given grid point at a given time step.

The four-dimensional weighting function for each observation i in (4.2.1) is rewritten

as

W(x; t) = wxy  w  wt (4:2:2)

The horizontal weighting function, wxy, is a Cressman-type spatial weighting function

dened by

wxy =

R2  D2

R2 + D2 0  D  R (4:2:2)

and

wxy = 0 D > R; (4:2:3)

where R is the radius of in
uence and D is the distance from the i th observation location

to the grid point. The vertical weighting function, w is also a distance-weighted function

dened by

w = 1 

jobs  j

R

jobs  j  R (4:2:4a)

w = 0 jobs  j > R; (4:2:4b)

where R is the vertical radius of in
uence and obs is the vertical position of the i th

observation. The temporal weighting function is given by

wt = 1 jt  t0j < =2 (4:2:5)

wt =

  jt  t0j

 =2

 =2  jt  t0j   (4:2:6)

33
wt = 0 jt  t0j > ; (4:2:7)

where t is the model-relative time, t0 is the model-relative time of the ith observation,

and  is the half-period of a predetermined time window over which an observation will

in
uence the model simulation.

For economy, the multi-level observations (soundings) used for obs nudging are usually

vertically interpolated to the model sigma surfaces at each observation location prior to

each simulation. Although the vertical component of the weighting function, w (4.2.4),

is also a distance-weighted function, the vertical radius of in
uence, R , can be dened

to be small (less than the spacing of the model layers) so that each observation above the

model surface layer in
uences only one sigma layer at a given location.

Figure 4.3 illustrates schematically the horizontal and temporal components of W

used for nudging to observations. The horizontal weighting function, wxy , is the Cressman

function given by (4.2.2) and (4.2.3). As shown in the top of the gure, the horizontal

radius of in
uence varies linearly in the vertical with pressure, from Rs at the surface to

the preset value R0 at a pressure level p0 representing the free atmosphere, where terrain

in
uences are assumed to be small. At pressures less than or equal to this user-dened

value, dened by default as 500 mb, the horizontal radius of in
uence is dened by default

as twice the value used in the surface layer, Rs. For example, if Rs = 100 km, R0 = 200

km. For certain situations, such as with upward propagating mountain-induced gravity

waves, the assumption of negligible terrain in
uence within the troposphere is invalid and

should be avoided.

As shown in the top of Fig. 4.3, the corrections computed at a given observation site

and vertical level above the surface layer (lowest model layer) are spread laterally along a

constant pressure level and thus across several sigma layers in regions of sloping terrain.

That is, for any given grid point within the horizontal radius of in
uence, the obs-nudging

correction in the horizontal direction is applied to the sigma layer which has a pressure

value closest to that of the observation.

Observations within the model surface layer are spread along constant sigma surfaces,

but with a modied Cressman function (dashed contours in the middle of Fig. 4.3) which

reduces the in
uence of an observation as a function of the surface pressure (terrain).

34
* *

*

*

*

***

*** ***

*

*

*

* *

** *

**

R

Rs

A

σ1

σ2

σ3

σ4

σ5

Z

X

P4 (

P3

P2

P1

B

P

*)

(

(

( )

)

)

Wxy

Y

X

A

z

z+∆z

z+∆z z+∆z z-∆z

z z

z

B

Rs

1.0

Wt

0.0

t0-τ t0 t0+τ

Figure 4.3: Schematic showing the horizontal weighting function, wxy, and the temporal weighting function, wt

, used for obs nudging.

35
Thus spreading the in
uence of surface-layer observations along the lowest sigma

ensures that the FDDA forcing near the surface in uneven terrain is continuous, and not

like a pebble skipping across a pond. For observations in the surface layer, the distance

factor D in (4.2.2) is replaced with Dm,

Dm = D + RSC1 m jps0  psj; (4:2:8)

where D is as dened above, Cm is a constant, and ps0 and ps are the surface pressures

at the observation location and the grid point, respectively. For example, Cm is typically

dened as 75 mb, and Rs is the surface-layer value for the horizontal radius of in
uence.

As the dierence in surface pressure between the observation location and the grid point

approaches Cm, the second term in (4.2.8) approaches R, and wxy tends to zero faster

for a given D. Therefore, the eect of assimilating surface-layer observations in the valley

(mountains) on grid point locations in the mountains (valley) will be much reduced. This

minimizes the possibility that observations in complex terrain will in
uence the model

solution in areas where they may not be representative. Also, the vertical weighting factor,

w, for these surface-layer observations is dened so that the vertical in
uence of the

surface-layer observations decreases linearly through the lowest 3 or so model layers (about

250 m AGL). As mentioned earlier, single-level data are retained better by the model if

assimilated through several vertical layers.

The temporal weighting function, wt (4.2.5-4.2.7), shown in the bottom of Fig. 4.3, is

nonzero during a preset time window centered about the observation time, t0. It determines

the time period over which the ith observation can in
uence the model simulation via

(4.2.1). In general, this time window can also be dened as a function of the pressure level

of the observation similar to the eect of the horizontal radius of in
uence, R, in (4.2.2).

Thus the nal correction to the model solution via obs nudging re
ects a weighted average

of all observations during the preset time window about the current time step and within

some three-dimensional neighborhood of each grid point.

36
5. Treatment of physical processes

5.1 Horizontal diusion

Two types of diusions are used to control nonlinear instability and aliasing. These

are a second-order diusion of the form

FH2 = p

KHr2; (5:1:1)

where  is any prognostic variable, and a more scale-selective fourth-order form

FH4 = p

K0Hr4; (5:1:2)

where

K0H = s

2KH (5:1:3);

The second order diusion is only used in the coarsest domain for the row and column of

the grid points next to the lateral boundaries, while the fourth-order form is used in the

interior of the coarsest domain as well as in the entire domain of any renement mesh.

The horizontal diusion coecient KH consists of a background value KH0 and a

term proportional to the deformation D

KH = KH0 + :5k

2s2D (5:1:4)

where k is the von Karman constant and D is given by (Smagorinski et al. 1965)

D =

"

@u

@x 

@v

@y 2

+



@v

@x +

@u

@y 2# 1

2

: (5:1:5)

A background value of KH is a function of grid size and time step, where

KH0 = 3: 103 x2

t

: (5:1:6)

Note that the horizontal operators r4 and r2 are applied on constant sigma surfaces. To

ensure computational stability, an upper limit of x2

64t

is imposed on KH

5.2 Dry Convective Adjustment

There may be situations in which super-adiabatic layers are produced in the model

atmosphere. When this happens, and there is no call to the Blackadar planetary boundary

37
layer parameterization, a simple scheme is used to remove any unstable layers. The scheme

operates on the entire sounding at once and conserves the vertical integral of internal and

potential energy. When the model lapse rate of potential temperature @

@p exceeds a critical

value 

@

@p 

c

, the sounding is adjusted so that (1) mass-weighted mean temperature is

unchanged, and (2) the potential temperature lapse rate after adjustment equals 

@

@p 

c

.

Given n layers in which the model potential temperature lapse rate exceeds the critical

value, the rst constraint gives

(Tn + Tn)n + (Tn1 + Tn1)n1 + ::: + (T1 + T1)1 = TXn

i1

i ; (5:2:1)

where Ti are the adjustments to be added to the temperature at layer i, Ti and i are

the temperature and thickness of the sigma layers, and Tis the mass weighted mean

temperature. The second constraint gives

(Ti + Ti)i  (Ti1 + Ti1)i1 =



@

@p

(pi  pi1) i = 2; :::; n; (5:2:2)

where i is the Exner function. There are n equations that can be solved for n variables

Ti. The Gaussian elimination method is used to solve the n n matrix system. After

adjustment, the entire sounding is rechecked for unstable layers.

The moisture in the adjusted layers is assumed constant in the vertical, i.e.,

qvi = qv = Pn

i=1

P

qvii n

i=1 i

(5:1:3)

5.3 Precipitation physics

MM5 has many dierent choices to treat precipitation physics. They are usually

divided into two dierent groups: explicit and implicit schemes. Explicit schemes treat

resolved precipitation physics while implicit schemes treat the non-resolved precipitation

physics. Both may be operating at a grid-point at the same time. A commonly used

terminology of \convective" versus \stable" precipitation is generally not acceptable on

ner grid-resolutions, where convective precipitation is quite often resolved. Hence in the

following subsections we will use resolved/non-resolved and explicit/implicit as common

terminologies. As two additional options, MM5 allows for dry runs, where moisture is

38
treated as a passive variable (no explicit and implicit schemes are applied). Another

option is a \fake dry run", where only the eects of the latent heat release are removed.

These 2 options do not require any further description and will not be discussed in the

following subsections.

5.3.1 Resolvable scale precipitation processes

These schemes are usually activated whenever grid-scale saturation is reached.

In other words, they treat resolved precipitation processes. The most simple way

that sometimes is still used on larger-scales, is to simply remove super-saturation as

precipitation and add the latent heat to the thermodynamic equation. More sophisticated

schemes carry additional variables such as cloud and rainwater (subsection 5.3.1.1), or

even ice and snow (subsection 5.3.1.2). Both schemes described next are enhancements of

MM4's original scheme (Hsie 1984).

5.3.1.1 Explicit treatment of cloudwater, rainwater, snow, and ice

This scheme optionally allows for ice-phase processes below 0 C, where cloud water

is treated as cloud ice and rain is treated as snow (Dudhia 1989). The equations for water

vapor, cloud water (ice) and rain water (snow) mixing ratios are given by the following

@pqv

@ t =  m2 

@puqv=m

@x +

@pvqv=m

@y  

@pqv_

@ + nhqvDIV

+ p



(  PRE  PCON  PI I  PID) + Dqv ; (5:3:1:1:1)

@pqc

@ t =  m

2 

@puqc=m

@x +

@pvqc=m

@y  

@pqc_

@ + nhqcDIV

+ p



(PID + PI I  PRC  PRA + PCON ) + Dqc; (5:3:1:1:2)

@pqr

@ t =  m

2 

@puqr=m

@x +

@pvqr=m

@y  

@p qr_

@ + nhqrDIV



@Vf gqr

@ + p



(PRE + PRC + PRA) + Dqr ; (5:3:1:1:3)

where PCON is condensation (and freezing for T < 0 C) of water vapor into cloud (ice) at

water saturation, PRA is accretion of cloud by rain (ice by snow), PRC is conversion of cloud

to rain (ice to snow) and PRE is evaporation (sublimation) of rain (snow). Additional ice

processes are PI I , the initiation of ice crystals, and PID sublimation/deposition of cloud

39
ice (Fig. 5.1). The fall speed of rain or snow is Vf . The term nh is 1 for nonhydrostatic

and 0 for hydrostatic simulations.

In all the relevant processes, Marshall-Palmer size distributions are assumed for rain

and snow and droplet fall speeds are taken to be of the form V (D) = aDb, where D is

the diameter. For rain, the Marshall-Palmer intercept parameter is N0 = 8 106 m4,

a = 841:99667 and b = 0:8 for V in m s1and D in meters, while for snow N0 = 2 107

m4, a = 11:72 and b = 0:41.

The saturated vapor pressure over water (in mb) is taken to be

esw = 6:112exp 17:67 

T  273:15

T  29:65  ; (5:3:1:1:4)

and for ice

esi = 6:11exp 22:514 

6150

T 

: (5:3:1:1:5)

The saturated water vapor mixing ratio is then given by

qs =

0:622es

p  es

:

PRC, the autoconversion term is represented by

PRC = max[k1(qc  qcrit); 0]; (5:3:1:1:6a)

for cloud to rain and by

PRC = max[(qc  Mmaxnc)=t; 0]; (5:3:1:1:6b)

for ice to snow, where k1 = 103s1, qcrit = 0:5 g kg1, Mmax = 9:4 1010 kg and nc

is given by Fletcher's (1962) formula for the number concentration of ice nuclei (kg1),

nc = 102exp[0:6(273:15  T )]=:

PI I , the initiation of ice crystals is given by

PI I = max[(M0nc  qc)=t; 0]; (5:3:1:1:7)

as long as sucient supersaturation over ice exists, where M0 = 1012 kg.

40
ice cloud

vapor

snow rain

PRF

PRA

PRC

PRA

PRC

PID

PII PCON

PRM

PMF

PRE PRE

Figure 5.1 Box diagram illustrating the processes in the moisture scheme for ice

(crystals), cloud(liquid), snow and rain. PCON, condensation/evaporation of

cloud; PRA, accretion; PRC, conversion; PID, deposition onto ice crystals; PRE,

evaporation for rain and deposition/sublimation for snow; PMF, melting/freezing

due to advection; PII, initiation of ice crystals; and PRM, melting of snow due to

fall.

41
PRA, the accretion rate is given by

PRA =

1

4 aqcEN0

(3 + b)

3+b

; (5:3:1:1:8)

where  is the gamma-function, E is the collection eciency (1 for rain and 0.1 for snow)

and  is given by

 =



N0w

qr 1=4

:

Here w is the mean density of rain or snow particles (1000 and 100 kg m3, respectively.)

PID, the deposition onto or sublimation of ice particles is found from

PID =

4Di(Si  1)nc A + B

; (5:3:1:1:9)

where

Si = qv=qsi;

A =

L2

s

KaRvT 2

; B =

1

qsi

:

Ls is the latent heat of sublimation, Ka is the thermal conductivity of air, Rv is the gas

constant for water vapor, and  is the diusivity of vapor in air. The mean diameter of ice

crystals, Di, is found from the mean mass, Mi = qc=nc, and the mass-diameter relation

for hexagonal plates from Rutledge and Hobbs (1983), Di = 16:3M1=2

i meters.

PRE, the evaporation of rain and sublimation/deposition of snow can be determined

from

PRE =

2N0(S  1)

A + B

"

f1

2 + f2 

a

 1=2

S1=3

c (5=2 + b=2)

5=2+b=2

#

; (5:3:1:1:10)

with the relevant N0, a, and b chosen for rain or snow, and S = Sw or Si. The denition

of A and B also change from the above for rain, substituting Lv for Ls and qsw for qsi . For

snow, 2 is replaced by 4. The values of f1 and f2 are 0.78 and 0.32 for rain and 0.65 and

0.44 for snow. The term in brackets represents a distribution-integrated ventilation factor,

F = f1 + f2S

1=3

c Re1=2

, with Sc = =, the Schmidt number, and Re = V (D)D=, the

Reynolds number, and  is the dynamic viscosity of air.

42
PCON , the condensation is determined as follows. Temperature, water vapor mixing

ratio and cloud water are forecast rst: these preliminary forecast values are designated

by T , q



v and q



c . We dene

M = q



v  q



vs;

where q



vs is the saturated mixing ratio at temperature T 

,

(1) if M > 0 (supersaturation),

PCON =

r1M

t

; (5:3:1:1:11a)

where

r1 =

1

1 + L2v q

vs Rv cpmT 2

;

(2) if M < 0 and qc > 0 (evaporation),

PCON = min 

r1M

t

; q



c t



; (5:3:1:1:11b)

(3) if M < 0 and qc = 0,

PCON = 0: (5:3:1:1:11c)

The PCON term is computed diagnostically so no iteration is needed.

Additionally, as snow falls through the 0 C level, it immediately melts to rain. This

process is given by

PRM = 

gVf qr p

: (5:3:1:1:12)

Advection of ice or snow downwards or of rain or cloud upwards through this level also

melts or freezes the particles, where

PMF = 

!(qc + qr )

p

: (5:3:1:1:13)

In both cases, the 0 C isotherm is taken to be at a full model level boundary. Melting

occurs at the level immediately below this boundary and freezing above it.

The latent heating is thus

Q_ = L(PRE + PID + PI I + PCON ) + Lm(PRM + PMF ); (5:3:1:1:14)

43
where L = Lv for T > 0 C and L = Ls for T < 0 C, and Lm = Ls  Lv .

The fall speed is mass-weighted and so is determined from

Vf = a

(4 + b)

6

b: (5:3:1:1:15)

The fall term in (5.3.1.1.3), the rain and snow prediction equation, may be calculated on

split time-steps, t

0 , in the explicit moisture routine. This ensures that Vft0=z < 1,

which is required for numerical stability. The size of t

0 is determined independently in

each model column based on the maximum value of Vft=z in the column, where t is

the model time step.

44
5.3.1.2 Mixed-Phase Ice Scheme

This scheme is based on the simple ice phase scheme described in the previous

subsection, but it does not immediately freeze or melt water and ice. Supercooled water

can exist below 0C in this scheme, as can unmelted snow exist above 0C. Separate arrays

are used to store vapor, cloud, rain, cloud ice and snow.

Homogeneous freezing of cloud water to cloud ice occurs immediately below -40C

and cloud ice melts immediately above 0C. Snow melts according to

PSM = 

2N0s Lf

Ka(T  T0) "

f1

2 + f2 

a

 1=2

S1=3

c (5=2 + b=2)

5=2+b=2

#

; (5:3:1:2:1)

where f1 = 0:78 and f2 = 0:31 (Rutledge and Hobbs 1984), and the other constants are

the ones relevant to snow in subsection (a). Evaporation of melting snow is modied to

use the values of A and B for rain as in (5.3.1.1.10).

Heterogeneous freezing of cloud water to cloud ice is also included following Bigg

(1953),

PCI = B0fexp[A0(T0  T )]  1g q2

c

wNc (5:3:1:2:2)

where A0 = 0:66K1, B0 = 100m3s1 and the number concentration of cloud droplets

per unit volume of air, Nc = 1010m3.

Sekhon and Srivastava (1970) determined that better comparison against observed

snow distributions can be obtained in theoretical studies if the slope intercept value for

the size distribution is expressed as

N0s(m4)=1:05R0:94 (5:3:1:2:3)

where, N0s is the slope intercept and R (m s1) is the snow fall rate. Thus a variable

intercept parameter replaces the constant N0s used in the simple ice scheme.

This can be expressed in terms of snow mixing ratio, qS , as

N0s = 8

<

:

1:05 "

1

qS 

S

qS b4 #0:949

=

;

4

0:94b+4

(5:3:1:2:4)

where, 1 = 6w

a(4+b)

.

45
5.3.2 Implicit cumulus parameterization schemes

5.3.2.1 The Kuo scheme

In this scheme, the amount of convection is determined by the vertically integrated

moisture convergence. The feedback to the larger scale (the vertical distribution of

heating and moistening), is determined with the help of the normalized vertical proles of

convective heating (Nh()) and moistening (Nm()), and a vertical eddy-
ux divergence

of water vapor associated with cumulus convection Vqf (). Therefore, equations (2.1.3),

(2.2.5) and (5.3.1.1.1) can be rewritten to include the convective-scale 
uxes as

@pT

@ t =  m

2 

@puT =m

@x +

@pvT =m

@y  

@pT _

@

+ p

 !

cp

+ p

 Lv

cpm

Nh()(1  b)gMt + DT ; (5:3:2:1:1)

@pT

@ t =  m2 

@puT =m

@x +

@pvT =m

@y  

@pT _

@ + T :DIV

+

1

cp



p

 Dp0

Dt  0gpw  Dp0  + p

 Lv

cpm

Nh()(1  b)gMt + DT ; (5:3:2:1:2)

@pqv

@ t =  m2 

@puqv=m

@x +

@pvqv=m

@y  + nhqvDIV

+ p



(  PRE  PCON  PI I  PID) +p

 bgMtNm()+pVqf () + Dqv ; (5:3:2:1:3)

where the vertically integrated moisture convergence Mt is

Mt =

m2

g  Z 1

0

rpV q _ v m

d: (5:3:2:1:4)

A portion (1  b) of Mt is assumed to condense and precipitate, where the remaining

fraction b is assumed to moisten the grid column. Following Anthes(1977), b is a function

of the mean relative humidity RH of the column, where

b = 2(1  RH) (5:3:2:1:5)

for RH  0:5, and b = 1 otherwise.

The vertical proles of heating and moistening

46
The normalized, nondimensional functions for the vertical proles of heating and

moistening and the divergence of the vertical eddy 
ux of water vapor are sub ject to the

constraints

Z 1

0

Nh()d = 1; (5:3:2:1:6)

Z 1

0

Nm()d = 1; (5:3:2:1:7)

Z 1

0

Vqf ()d = 0: (5:3:2:1:8)

Anthes et al. (1987) assume simple relationships for these functions, which are

derived from budget studies. For the convective heating prole, Nh, they observe that

the convective heating often has a parabolic shape with a maximum in the upper half of

the cloud. Hence

Nh() = a1x2 + a2x + a3; (5:3:2:1:9)

where

x = ln (5:3:2:1:10)

with the boundary conditions:

Nh()=0; at xb = lnb; and xu = lnu (5:3:2:1:11)

at cloudbase (b) and cloud top (u), and

N0h() = @Nh()

@

= 0 (5:3:2:1:12)

at x, which is dened as

x =

xu + xb

2

; (5:3:2:1:13)

where subscripts u and b stand for the top and the base of the cloud. Using (5.3.2.1.6), a1

can be shown to be

a1 =

2

x3u  x3

b + x2uxb  xux2b

: (5:3:2:1:14)

The vertical moistening prole, Nm(), is simply given following Anthes (1977) as

Nm() = (1  RH())qv s()

R 1ktop (1  RH(0 ))qvs (0)d0 : (5:3:2:1:15)

47
Divergence of the Vertical Eddy Flux of Water Vapor Vqf ()

The divergence of the vertical eddy 
ux of water vapor is dened as

Vqf () = @_

0 q0v

@ : (5:3:2:1:16)

If one assumes a small fraction of convective cloud cover, and the cloud vertical motion _c

is much larger than the larger-scale vertical motion, _ (5.3.2.1.16) can be rewritten as

Vqf () = a

1  a

@

@ [ _c(qvc  qv )]; (5:3:2:1:17)

where qvc is the mixing ratio in the cloud.

According to Anthes (1977), the fractional coverage a is calculate using

a =

(1  b)gMt

R

p

0 !c @qvc @p + @qvc @te 

dp

; (5:3:2:1:18)

which is the ratio between the grid-average condensation rate and that of a single cloud.

The term @qvc @te represents the contribution to the rate of change of cloud-mixing ratio

by entrainment (Anthes 1977). Anthes et al. (1987) assume a typical value for the

denominator of approximately 4:3 103cb s1 and then rewrite (5.3.2.1.17) as

Vqf () = (1  b)gMt

4:3 103

@

@ [ _c(qvc  qv )]: (5:3:2:1:19)

For further simplication, Anthes et al. (1987) next assume that _c also has a parabolic

shape and can be expressed as

_ c = c1x

2 + c2x + c3; (5:3:2:1:20)

where x = lnp, and _c = 0 at cloud-top and base. Furthermore, qvc  qv is assumed to

have a parabolic prole with pressure

qvc  qv = b1x

2 + b2x + b3 (5:3:2:1:21)

with

x = ln[(1  )(100  pt) + pt ]: (5:3:2:1:22)

48
The procedure

The simple procedure can be summarized as follows:

1. Compute Mt from (5.3.2.1.4)

2. Check whether Mt  3: 105kg m2 s1, a critical threshold value.

3. Check the model sounding for convective instability to see if convection is possible.

4. Determine cloud top and base from sounding.

5. Check whether cloud-depth is larger than a critical value (  :3)

6. Calculate the normalized vertical prole functions

7. Calculate _c on the full  levels from (5.3.2.1.20)

8. Compute qvc  qv from (5.3.2.1.21)

9. Calculate Vqf from (5.3.2.1.19)

49
5.3.2.2 A modied Arakawa-Schubert scheme

The version of the Arakawa-Schubert scheme used here was developed by Grell (1993).

In contrast to the original scheme (Arakawa and Schubert 1974, AS), it includes moist

convective-scale downdrafts. Other changes have been implemented to also allow the

scheme to be used successfully in mesoscale models in mid-latitudes (Grell et al. 1991).

To simplify the description we have adapted a terminology originally introduced by Betts

(1974), which splits the parameterization problem from the modeling view in three parts:

static control, dynamic control, and feedback. The static control includes usually a

cloud-model and calculates cloud thermodynamic properties, the dynamic control is what

determines the amount and location of the convection, and the feedback determines the

vertical distribution of the integrated heating and moistening.

Static control

As with all commonly used one dimensional steady state cloud models (plumes,

bubbles, or jets), our AS scheme makes use of the assumption that entrainment occurs

over the depth of the buoyant element according to the entrainment hypothesis

 =

1

m(z)

@m(z)

@z



:2

r

; (5:3:2:2:1)

where  is the total net fractional entrainment rate of the buoyant element, m its mass 
ux

(mu for updraft, md for downdraft), and r its radius. Following AS, the dependence on the

radius is not explicitly used. However, implicitly, the radius of the cloud is assumed to be

constant. Detrainment was originally only assumed to happen at the cloud top, but this

assumption may easily be varied (Houze et al. 1979, Lord 1978) by dening a fractional

detrainment rate, ud, and rewriting (5.3.2.2.1) for the updraft of cloud type  as

u = ue  ud =

1

mu(z)

@mu(z)

@z

=

1

mu(; z) @mu(; z)@z 

ent 



@mu(; z)

@z 

det

 (5:3:2:2:2)

where ue is the gross fractional entrainment rate, and u is the total net fractional

entrainment rate of the updraft. Subscripts ent and det indicate changes due to

entrainment and detrainment, respectively. Looking at the budget of a thermodynamic

50
variable in an innitesimal layer of the updraft we get

@muu

@z =



@mu

@z 

e ~ 



@mu

@z 

d

u + Su: (5:3:2:2:3)

Together with (5.3.2.2.2) this leads to the steady state plume equation

@u(; z)

@z = ue( ~(z)  u(; z)) + Su (5:3:2:2:4)

where  is a thermodynamic variable, the tilde denotes an environmental value, and

subscript u denotes an updraft property. S stands for sources or sinks. Similarly, for

the downdraft, we can rewrite equations (5.3.2.2.2) and (5.3.2.2.4) as

d = de  dd = 

1

md(z)

@md(z)

@z

= 

1

md(z) @md(z)@z 

ent 



@md(z)

@z 

det

 (5:3:2:2:5)

and

@d(z)

@z = de( ~(z)  d(z)) + S; (5:3:2:2:6)

where subscript d denotes a downdraft property. For moist static energy

~h(z) = CpT~(z) + gz + Lq~(z); (5:3:2:2:7)

equations (5.3.2.2.4) and (5.3.2.2.6) simply become

@hu(; z)

@z = ue(h~

(z)  hu(; z)) (5:3:2:2:8)

and

@hd(z)

@z = de[h~

(z)  hd(z)]: (5:3:2:2:9)

Next, for the moisture budget of the updraft, we use

u = qu(; z) + ql(; z) (5:3:2:2:10)

and

Su = c0mu(; z)ql (; z): (5:3:2:2:11)

51
Here Su is the total water that is rained out, c0 is a rainfall conversion parameter and

could be a function of cloud size or wind shear, ql is the suspended liquid water content of

the cloud, and qu is the water vapor mixing ratio inside the updraft. Equation (5.3.2.2.4)

can then be rewritten as

@(qu(; z) + ql(; z))

@z = ue(~q(z)  qu(; z)  ql(; z)) + Su: (5:3:2:2:12)

For the downdraft, the equation for the water vapor reads

@qd(z)

@z = de[~q(z)  qd(z)] + Sd: (5:3:2:2:13)

Sd here is a source; namely the evaporation of rain. Assuming saturation in the updraft

and downdraft, we can make use of the approximate equation

qc(; z)=~q

 +




1 + 


1

L

[hc(; z)  ~h(z)]; (5:3:2:2:14)

where


 =

L

cp 

@q~

@T 

P

(5:3:2:2:15)

the asterisk denotes a saturated value, and hc here stands for the moist static energy in the

cloud (updraft or downdraft), if saturation is assumed. Next, to arrive at a usable closure,

the up- and down-draft mass 
uxes are normalized by the updraft base (mb()) mass 
ux,

and the downdraft base m0() mass 
ux of a subensemble. Hence, for the updraft,

mu(; z) = mb()u(; z) (5:3:2:2:16)

and

ue  ud =

1

u(z)

@u(z)

@z : (5:3:2:2:17)

Equivalently, for the downdraft we may write

md(z) = m0()d(; z) (5:3:2:2:18)

and

de  dd =

1

d(z)

@d(z)

@z : (5:3:2:2:19)

52
Here, m0 is the mass 
ux at the originating level and d, as u in equation (5.3.2.2.16), is

the normalized mass 
ux prole.

To leave only one unknown variable, we follow Houze et al. (1979) and make the

originating mass 
ux of the downdraft a function of the updraft mass 
ux and reevaporation

of convective condensate. Therefore, the condensate in the updraft

Cu()d = mbd Z zT

zB

u(; z)Sudz)  I1mbd (5:3:2:2:20)

is apportioned according to

Cu()d = (Rc() + Ed())d = (() + ())Cu()d; (5:3:2:2:21)

where  +  = 1 and Ed, the evaporation of condensate in the downdraft for cloud type

, can be written as

Edd = m0()d Z z0

0

d(; z)Sddz  I2m0d: (5:3:2:2:22)

From equations (5.3.2.2.20-5.3.2.2.22) we see that

Edd = Cud = I1mbd = I2m0d (5:3:2:2:23)

and hence

m0() = ()I1mb()

I2() = ()mb(): (5:3:2:2:24)

Here 1 is the precipitation eciency. Following Fritsch and Chappell (1980), it is made

dependent on the windshear.

To solve the above equations we need to specify boundary conditions as well as make

some arbitrary assumptions. For the updraft we assume

hu(zb) = MAX(h~(z)) (5:3:2:2:25)

with

z  zb

and

hu(; zT ) = h~(zT ); (5:3:2:2:26)

53
where the asterisk denotes a saturation value. Similarly, for the downdraft,

hd(z0) = MIN(h~(z)): (5:3:2:2:27)

Physically, for both updraft and downdraft, we allow for maximum buoyancy. The

boundary conditions for the updraft are dierent than in the original scheme, which had

a rigid dependence on the planetary boundary layer height. In the original scheme, the

mixed layer was assumed to be well mixed, and the cloud base was located on top of the

mixed layer. In semi-prognostic tests (Grell et al. 1991) large variations of moist static

energy proles were found in very low levels of the troposphere. This was caused by cold

downdraft out
ow. Naturally, the in
ow to an updraft will not be a mixture of downdraft

air and the more buoyant air; it is more likely the air with high moist static energy from

the layer above the downdraft out
ow. Furthermore, compensatory subsidence should only

continue to the level from which the updraft draws its air. Compensatory uplifting may

be required in very low layers of the troposphere because of the downdraft mass 
ux.

Feedback

The feedback to the larger-scale environment is expressed in a convenient form as



@s

@ t 

cu =

1:



@

@z FsLl ; (5:3:2:2:28)



@q

@ t 

cu = 

1:



@

@z Fq+l  R; (5:3:2:2:29)

where s is the dry static energy (s = cpT + gz). The convective-scale 
uxes within a grid

box are dened as

FsLl  Fs  LFl (5:3:2:2:30)

Fq+l  Fq + Fl (5:3:2:2:31)

where Fs is the 
ux of dry static energy, Fq is the 
ux of water vapor, and Fl is the 
ux

of suspended cloud liquid water. These are dened as

Fs(z)  +

Z



u(; z)[su(; z)  s(z)]mb()d



Z



d(; z)[sd (; z)  s(z)]m0()d

(5:3:2:2:32)

54
Fq (z)  +

Z



u(; z)[qu(; z)  q(z)]mb()d



Z



d(; z)[qd(; z)  q(z)]m0()d

(5:3:2:2:33)

Fl(z) 

Z



u(; z)l(; z)mb ()d (5:3:2:2:34)

The rainfall (convective-scale sink of cloud water) is dened as

R(z)  +

Z



u(; z)c0()l(; z)mb ()d



Z



d(; z)qe (; z)m0()d

(5:3:2:2:35)

Here qe is the amount of moisture that is necessary to keep the downdraft saturated. The

second term on the righthand sides is due to downdrafts and is zero above the downdraftoriginating level. Below the updraft-air originating level, the rst term on the right-hand

sides is zero and only downdrafts aect the larger-scale environment. Below the updraft-air

originating level, the convective-scale 
uxes due to updrafts are zero. Between the updraftair-originating level and the level of free convection (the LFC), Fl and R are set to zero.

Since no liquid water is assumed to be in the environment as the downdraft, the downward


ux due to updrafts as well as downdraft 
uxes in equation (5.3.2.2.33) are zero. Schubert

(1974) showed that convection will not increase the total moist static energy per unit area

in a column. In essence, only precipitation can change the dry static energy budget and

the total mass of water vapor. All variables in the 
ux terms can be determined from the

equations for the static control, except mb(). This is determined in the dynamic control,

which incorporates the closure assumption of the scheme and is described next.

Dynamic control

Arakawa-Schubert rst introduced the cloud work function, which is an integral

measure of the buoyancy force associated with a subensemble. Starting with

dwu

dt = Bu  Fr =

dwu

dz

dz

dt =

d

dt

d

dz

w2u

2 =

1

wu ddt

w2u

2

; (5:3:2:2:36)

55
where Bu is the acceleration due to buoyancy and Fr the deceleration due to friction, and

multiplying equation (5.3.2.2.36) by u(; z)wu(; z), gives

d

dt u w2u2 = uwu(Bu  Fr ): (5:3:2:2:37)

Integrating over the depth of the updraft and using mu = uwu = mbu yields

d

dt Z zT

zb u w2u

2

dz = mb()Z zT

zb uBudz  Du; (5:3:2:2:38)

where D is the updraft-scale kinetic energy dissipation. Equation (5.3.2.2.38) can be

written in the symbolic form

d

dtKEu = Au()mb()  Du(); (5:3:2:2:39)

where Au() is a measure of the eciency of kinetic energy generation inside the cloud

and is called the cloud work function. It can also be written as

Au() = Z zT

zB

g

CpT (z)

u(; z)

1 + 


(hu(; z)  h~(z))dz ; (5:3:2:2:40)

where 
 is dened as in equation (5.3.2.2.15). As with equations (5.3.2.2.36-5.3.2.2.38),

dening a kinetic energy generation inside the downdraft leads to

d

dtKEd = Ad()m0()  Dd(); (5:3:2:2:41)

where Ad, the measure of the eciency of kinetic energy generation inside the downdraft,

can be written as

Ad() = Z zsur

z0

g

CpT (z)

d(; z)

1 + 


(h~(z)  hd(; z))dz : (5:3:2:2:42)

Note that dry static energy instead of moist static energy would have to be used if

subsaturation is assumed. We can combine equation (5.3.2.2.39) and (5.3.2.2.41) and

then make use of (5.3.2.2.24) to yield

d

dtKEtot = Atot()mb()  Dtot(); (5:3:2:2:43)

56
where

Atot() = Au() + ()Ad() (5:3:2:2:44)

is the total cloud work function which was redened as a measure of the eciency of kinetic

energy generation in updrafts as well as downdrafts. Next, AS separated the change of the

cloud work function into two parts: One is due to the change in the larger-scale variables



dAtot

dt 

LS

 F (); (5:3:2:2:45)

and one is due to the modication of the environment by the clouds. Since the cumulus

feedback on the larger-scale elds is a linear function of mb, this term can be written in

the symbolic form



dAtot

dt 

CU



Z



K(; 0 )mb(0)d: (5:3:2:2:46)

Therefore

dAtot

dt = F () + Z



K(; 0 )mb(0)d; (5:3:2:2:47)

where K(; 0 ) are the kernels. The kernels are an expression for the interaction

between clouds (updrafts and downdrafts). Equation (5.3.2.2.47) is solved with a linear

programming method (Lord 1978).

In the original version of the Arakawa-Schubert scheme, the fractional entrainment

rate was the parameter which characterized the cloud. In later papers, the cloud-top

detrainment level was chosen instead. If a ne vertical resolution is assumed, the second

choice will most likely be better numerically, since no interpolation is necessary at the cloud

tops. However, in the extremely unstable environment of the mid-latitudes, it is sometimes

impossible to calculate \clouds" with cloud tops in the unstable layers. Entrainment rates

would have to be extremely large to stop cloud growth. We therefore chose the fractional

entrainment rate as the spectral parameter.

The procedure

The cloud base is a function of time and space. However, at a specic grid point the

cloud base will be the same for every member of the subensemble. We also distinguish

among an updraft-air originating level, zu, a downdraft-air originating level, z0, a cloud

57
base level, zb (the LCL), and a level of free convection, zbc (LFC). Here, zu is determined

from condition (5.3.2.2.25) and determines the thermodynamic properties of the updraft

from cloud type i. The air becomes saturated at zb; condensation will start, but no

convection can occur yet because the buoyancy is negative. In some instances this level

could be the same as the LFC. The LFC is of great importance since this is the level at

which the static control starts the calculations of individual convective elements. Since the

air that feeds the cloud originates below the LCL, compensatory subsidence is allowed to

reach the originating level of the updraft air.

For the downdraft, the originating level is also a function of time and space. If the

downdraft exists, it will always reach the surface.

For updraft and downdraft in layer k the mass budgets are dened as

eu(k; i)  du(k; i) = u(k + :5; i)  u(k  :5; i) (5:3:2:2:48a)

and

ed(k; i)  dd(k; i) = d(k + :5; i)  d(k  :5; i); (5:3:2:2:48b)

where entrainment for the updraft and downdraft is dened as

eu(k; i) = uezd u(k + :5; i) (5:3:2:2:49a)

ed(k; i) = dezd d(k  :5; i) (5:3:2:2:49b)

and detrainment is dened as

du(k; i) = udzd u(k + :5; i) (5:3:2:2:50a)

dd(k; i) = ddzd d(k  :5; i): (5:3:2:2:50b)

Combining the above three equations for the updraft and downdraft yields

u(k  :5; i) = u(k + :5; i)(1: + uezd  udzd) (5:3:2:2:51a)

for the updraft and

d(k + :5; i) = d(k  :5; i)(1: + dezd  ddzd) (5:3:2:2:51b)

58
for the downdraft. Here we dene zd = z(k + :5)  z(k  :5). The discretized form for

the downdraft moist static energy budget reads

ed(k; i)h~(k)  dd(k; i) hd(k + :5; i)  hd(k  :5; i)

2 = d(k + :5; i)hd(k + :5; i)  d(k  :5; i)hd(k  :5; i)

: (5:3:2:2:52)

Using equations (5.3.2.2.48)-(5.3.2.2.51) in equation (5.3.2.2.52) leads to

hd(k + :5; i) = hd(k  :5; i)(1:  :5ddzd) + dezd

~h(k)

1: + dezd  ddzd + :5ddzd

: (5:3:2:2:53)

The moisture budget for the downdraft is developed in several steps. First, the downdraft

water vapor mixing ratio before evaporation, but after entrainment, is calculated. This is

done using

qd(k; i) = qd(k  :5; i)(1:  :5ddzd) + dezd q~(k)

1: + dezd  ddzd + :5ddzd

: (5:3:2:2:54)

Next, equations (5.3.2.2.14) and (5.3.2.2.15) give the mixing ratio, qvd, that the updraft

or downdraft would have if saturated. Hence, the amount of moisture that is necessary to

keep the downdraft from cloud type i saturated in layer k is

qe(k; i) = [qd(k; i)  qvd(k; i)]: (5:3:2:2:55)

Next we check whether the updraft produces enough rain to sustain saturation in the

downdraft by requiring that

X

c0z(k) u(k  :5; i)ql(k  :5; i) X(i)z(k) d(k + :5; i)qe(k; i) > 0: (5:3:2:2:56)

If this is not the case, a downdraft is not allowed to exist.

Having dened the discretized versions of the equations from the static control, we

now can describe the procedure.

Using the larger-scale temperature and moisture elds (T0; q0) at time t0, and given

a functional or empirical relationship for d, de, and dd, the equations from the static

control are used to calculate ue; hu(z; i); hd(z; i); qu(z; i); qd (z; i); u(z; i), and d(z; i) for

cloud type i. These are needed to determine the total cloud work function Atot using

Atot(i) = Au(i) + Ad(i): (5:3:2:2:57)

59
The discretized versions of equations (5.3.2.2.40) and (5.3.2.2.42) that are used to

determine the cloud work functions for updrafts and downdrafts are

Au(i) =

k=Xktop

k=LF C"

g

cp T (k  :5) u(k  :5; i)

 

hu(k  :5; i)  h~(k  :5)

1 + 
(k  :5) !

 (z(k  1)  z(k))#

(5:3:2:2:58)

and

Ad(i) =

k

X=sur

k=z0

"

g

cp T (k  :5) d(k  :5; i)

 

hd(k  :5; i)  ~h(k  :5)

1 + 
(k  :5) !

 (z(k)  z(k  1))#

: (5:3:2:2:59)

The kernels of cloud type i are by denition the changes of the cloud work functions due

to another subensemble, i

0 . Thus, following Lord (1978), T0 and q0 are modied by an

arbitrary amount of mass 
ux, m0bt

0, from the i0 subensemble. This is done for every

possible subensemble and can be written in the symbolic form

T0(k; i) = T(k) + i0 (T(k))m0bt

0 ; (5:3:2:2:60)

q

0(k; i)=q(k) + i0 (q(k))m0bt0 : (5:3:2:2:61)

The  terms, which are changes per unit mb(i), are easily calculated from budget

considerations as in Lord (1978). With the downdraft terms, the moist static energy

budget of layer k and cloud type i becomes

p(k)

g

i0(h~(k; i)) = + (u(k  :5; i)  (i)d(k  :5; i))h~(k  :5)

 (u(k + :5; i)  (i)d(k + :5; i))h~(k + :5)

 (eu(k; i) + (i)ed (k; i))h~(k)

+ du(k; i) hu(k + :5; i) + hu(k  :5; i)

2

+ (i)dd (k; i) hd(k + :5; i) + hd(k  :5; i)

2

; (5:3:2:2:62)

60
where eu(k; i)anddu(k; i) are the entrainment and detrainment for the updraft, and p(k)

is dened by p(k) = p(k + :5)  p(k  :5). A simple physical interpretation of the terms

on the righthand side can be understood by looking at Fig. 5.2. The rst term is the

subsidence on top of the layer, the second is the subsidence on the bottom of the layer. This

subsidence is an environmental compensatory mass 
ux due to the updraft and downdraft

mass 
uxes inside the cloud. Note that below zu the \compensatory subsidence" may be

compensatory uplifting, since in that case only downdrafts exist. The third term represents

entrainment into the updraft and downdraft; the fourth term represents detrainment from

the edges of the updraft; the fth term represents detrainment from the edges of the

downdraft.

For the moisture budget,

p(k)

g

i0 (~q(k; i)) = + (u(k  :5; i)  (i)d(k  :5; i))~q(k  :5)

 (u(k + :5; i)  (i)d(k + :5; i))~q(k + :5)

 (eu(k; i) + (i)ed (k; i))~q(k)

+ du(k; i) qu(k + :5; i) + qu(k  :5; i)

2

+ (i)dd (k; i) qd(k + :5; i) + qd(k  :5; i)

2

: (5:3:2:2:63)

At the cloud top, downdrafts have no eects and updrafts detrain all their mass.

p(ktop)

g

i0(h~(ktop; i)) =  u(ktop + :5; i)h~(ktop + :5)

 eu(ktop; i)h~(ktop)

+ du(ktop; i) hu(ktop + :5; i) + hu(ktop; i)

2 + u(ktop; i)hu(ktop; i)

(5:3:2:2:64)

and

p(ktop)

g

i0(~q(ktop; i)) =  u(ktop + :5; i)~q(ktop + :5)

 eu(ktop; i)~q(ktop)

+ du(ktop; i) qu(ktop + :5; i) + qu(ktop; i)

2 + u(ktop; i)qu(ktop; i)

:

(5:3:2:2:65)

61
[η(k - .5, i) - ε(i)ηdo (k - .5, i)] ψ (k - .5)

[e(k,i) + ε (i)edo(k,i)]ψ(k)

[η(k + .5,i) - ε(i)hdo(k+ .5,i)]ψ(k + .5)

+ε(i)ddo(k,i) ψdo(k - 5,i) + ψdo(k + 5,i)

2

+d(k,i) ψc(k - .5,i) + ψc(k + .5,i)

2

Figure 5.2 Illustration of budget for thermodynamic variable ψ in layer k.

k - .5

k

k + .5

Cloud

Updraft

Originating

Level

Downdraft

Originating

Level

Figure 5.3 Conceptual picture of convection parameterized in Grell scheme.

62
Here p(ktop) = p(ktop+:5)p(ktop:5). Note that in the fourth term we have included

the detrainment of all the cloud mass at the cloud top. Finally, at the surface

p(ksur)

g

i0(h~(ksur; i)) =  (i)d(ksur  :5; i))h~(ksur  :5)

+ (i)d(ksur; i)hd(ksur; i)

 (i)ed (ksur  :5; i)h~(ksur  :5)

+ (i)dd (ksur; i) hd(ksur; i) + hd(ksur  :5; i)

2

(5:3:2:2:66)

and

p(ksur)

g

i0(~q(ksur; i)) =  (i)d(ksur  :5; i))~q(ksur  :5)

+ (i)d(ksur; i)qd(ksur; i)

 (i)ed(ksur  :5; i)~q(ksur  :5)

+ (i)dd(ksur; i) qd(ksur; i) + qd(ksur  :5; i)

2

; (5:3:2:2:67)

with p(ksur) = p(ksur + :5)  p(ksur  :5). Here, the rst term is the compensatory

environmental mass 
ux, the second term is the detrainment of all downdraft air at the

bottom, the third term is entrainment into the downdraft, and the fourth term is the

detrainment of air around the downdraft edges.

The new thermodynamic elds, T0

0(k; i0) and q0

0(k; i0), are then used again from the

static control to calculate new cloud properties and a new cloud work function, A0tot(i

0 ; i).

Note that T 00 and q

00 are now functions of the subensemble i0 . From the denition of the

kernel we then can calculate the kernels simply as

K(i; i0) = A0tot(i

0 ; i)  Atot(i)

mb\t\

: (5:3:2:2:68)

Next, we go back to the original elds and modify those with the large-scale advective

changes to get

T \(k) = T0 +



@T

@ t 

ADV

t (5:3:2:2:69)

and

q\(k) = q0 +



@q

@ t 

ADV

t; (5:3:2:2:70)

63
where (5.3.2.2.69) and (5.3.2.2.70) are applied over t= 30 min. The double prime

quantities are then used again by the static control, which will calculate new cloud

properties, and so new cloud work functions, Atot\(i), will be determined. Next, the

large-scale forcing (by denition the change of the cloud work function due to large-scale

eects only) is calculated using

F (i) = Atot\(i)  Atot(i)

t

: (5:3:2:2:71)

The large-scale forcing and the kernels are then both used by the dynamic control to

estimate the cloud base mass 
ux distribution function, mb, using an IMSL subroutine

to solve the linear programming problem. Finally, the feedback to the larger-scale

environment is simply given by



@T (k)

@ t 

CU =

i

0MAX X

i0=1



0

i(T (k))mb(i

0 ) (5:3:2:2:72)

and



@q(k)

@ t 

CU =

i

0MAX X

i0=1



0

i(q(k))mb(i

0 ); (5:3:2:2:73)

where the precipitation can be calculated using

P =

i

0MAX X

i0=1

k=Xktop

k=1

c0z(k)ql(k + :5; i)mu(k + :5; i)



i

0MAX X

i0=1

k=Xktop

k=1

z(k)qev (k + :5; i)md(k + :5; i)

: (5:3:2:2:74)

5.3.2.3 The Grell scheme

This is a very simple scheme that was constructed to avoid rst-order sources of

errors (Grell 1993). The very simplistic conceptual picture of how this parameterization

is envisioned to function is shown in Fig. 5.3. Clouds are pictured as two steady-state

circulations, caused by an updraft and a downdraft. There is no direct mixing between

cloudy air and environmental air, except at the top and the bottom of the circulations.

The cloud model that is used to calculate cloud properties in this scheme is formulated

64
with only a few equations. Mass 
ux is constant with height, and there is no entrainment

or detrainment along the cloud edges. We can simply write

mu(z) = mu(zb) = mb (5:3:2:3:1)

and

md(z) = md(z0) = m0 (5:3:2:3:2)

for the mass 
ux of the updraft (mu) and the downdraft (md). Here mb and m0 are simply

the mass 
uxes of the updraft and downdraft at their originating level. If it is assumed that

the conditions at originating levels are given by the environment, for any thermodynamic

variable , the budget inside the cloud simply becomes

u(z)=~(zb) + Su(z) ; (5:3:2:3:3)

and

d(z)=~(z0) + Sd(z); (5:3:2:3:4)

where  is a thermodynamic variable, the tilde denotes an environmental value, and S

stands for sources or sinks. For moist static energy

~h(z) = CpT~

(z) + gz + Lq~(z); (5:3:2:3:5)

equations (3) and (4) simply become

hu(z) = h~(zb) (5:3:2:3:6)

and

hd(z) = ~h(z0): (5:3:2:3:7)

For the moisture budget of the updraft we can make use of the approximate equations

(5.3.2.2.14) and (5.3.2.2.15) to calculate the mixing ratio inside the cloud if saturation is

assumed. Together with equations (5.3.2.3.3) and (5.3.2.3.4), this will give us Su and Sd,

the condensation and evaporation. Note also that no cloud water is assumed to exist; all

water is converted to rain.

65
Given boundary conditions, equations (5.3.2.3.1)-(5.3.2.3.7) have two unknowns, mb,

and m0. In order to leave only one unknown variable, the originating mass 
ux of the

downdraft is made a function of the updraft mass 
ux and the reevaporation of convective

condensate, as in the previous section (see equations (5.3.2.2.20)-(5.3.2.2.24)). Therefore,

m0 =

I1mb

I2 = mb: (5:3:2:3:8)

Here, 1   is the precipitation eciency. To specify boundary conditions, we assume

hu(z) = hu(zb) = MAX(h~(z)); (5:3:2:3:9)

with z  zb; and

hu(zT ) = ~h(zT ); (5:3:2:3:10)

where the asterisk denotes a saturation value. Similarly, for the downdraft,

hd(z) = hd(z0) = MIN(h~(z)): (5:3:2:3:11)

Physically, for both, updraft and downdraft, we allow for maximum buoyancy. For this

deep convection scheme, the cloud base for the updraft is not restricted to the boundary

layer, but can be anywhere in the troposphere.

Feedback to the larger-scale equations

To avoid zero-order sources of errors, the feedback must include the cooling eects

of moist convective downdrafts. Furthermore, lateral mixing should never be excessive,

especially if the cloud properties have been calculated with a steady-state cloud model.

Keeping in mind the conceptual picture in Fig. 5.3, the feedback for this scheme is entirely

determined by compensating mass 
uxes and detrainment at cloud top and bottom.

Conceptually, no averaging (such as the normally used top-hat or Reynolds averaging

methods) is necessary. This does not mean that scale-separation is not required, but for

this parameterization it is not necessary to assume that the fractional area coverage is

very small. Note, however, that any parameterization can only make sense if a clear scale

separation exists. None of the parameterized eects may be resolved by the larger-scale.

66
Assuming that the conceptual picture in Fig. 5.3 happens in only one grid box, we can

express the changes caused by the convection as

 

@~h(k)

@ t !

CU =

@hu(z)mb

@z 

@~h(z)mb

@z 

@hd(z)m0

@z +

@~h(z)m0

@z (5:3:2:3:12)

and



@q(k)

@ t 

CU =

@qu(z)mb

@z 

@q~(z)mb

@z 

@qd(z)m0

@z +

@q~(z)m0

@z : (5:3:2:3:13)

Because of the simplicity of the static control, these equations can be further simplied to

give

 

@~h(k)

@ t !

CU = mb

@~h(z)

@z (1  ) + mb( @hu(z)@z  

@hd(z)

@z ) (5:3:2:3:14)



@q(k)

@ t 

CU = mb

@q~(z)

@z (1  ) + mb( @qu(z)@z  

@qd(z)

@z ): (5:3:2:3:15)

The rainfall is dened as

R  I1mb(1  ): (5:3:2:3:16)

The second term on the righthand sides of equations (5.2.2.3.14) and (5.2.2.3.15) are due

to downdrafts and are zero above the downdraft originating level. Below the updraft-air

originating level, the rst term of the right-hand sides are zero and only downdrafts aect

the larger-scale environment. All variables in the 
ux terms can be determined from the

equations of the static control, except mb. Dynamic control

Because of the simplicity of the above equations, many closure assumptions can be

used. The most simple closure is a Kuo-type assumption, which relates the rainfall rate

to the moisture convergence. However, more applicable seems to be a stability closure.

Again we have two choices. We could assume that the clouds will remove the available

buoyant energy as in other mesoscale parameterizations, or that the clouds will stabilize

the environment as fast as the larger-scale (or also sub-grid scale) destabilizes it, or even a

mixture of both. Although both assumptions are easily implemented, we chose the closure

which depends on the rate of destabilization. In this closure the change of the available

67
buoyant energy due to convection osets the changes due to other eects (larger-scale

destabilization as well as sub-grid scale destabilization), yielding



dABE

dt 

OTH = 



dABE

dt 

CU

: (5:3:2:3:17)

Next, the change due to the convection is normalized in terms of the mass 
ux to read



dABE

dt 

CU

 mb

dABE

dt 

NCU

; (5:3:2:3:18)

where subscript NCU denotes the change of the available buoyant energy due to a cloud

normalized by the cloud-base mass 
ux. Equations (5.3.2.3.17) and (5.3.2.3.18) are used

to calculate mb.

The Procedure

This section describes in detail the procedure necessary to calculate the convective

feedback. First, we will explain the very simplistic approach to calculate a normalized

feedback, then we will describe how the closure assumption determines the mass 
ux.

Using the larger-scale temperature and moisture elds (T0; q0) at time t0,

hu(z); hd(z); qu(z); qd (z) are simply arrived at (see equations (5.3.2.3.6)-(5.3.2.3.10). The

rst calculation is the determination of the integrals I1 and I2 (calculated as residuals

using equations (5.3.2.3.8) and (5.3.2.3.9). The next step is then to estimate the convective

changes per unit mass 
ux (before knowing the actual mb's). This is done by estimating

the net change of a thermodynamic variable  in a layer k by using

p(k)

g

( ~(k)) = (1  )( ~(k  :5)  ~(k + :5)); (5:3:2:3:19)

where p(k) is dened by p(k) = p(k+:5)p(k:5). This subsidence is an environmental

compensatory mass 
ux due to the updraft and downdraft mass 
uxes inside the cloud.

Note that below zu the \compensatory subsidence" may be compensatory uplifting, since

in that case only downdrafts exist.

At the cloud top,

p(ktop)

g

( ~(ktop)) = ~(ktop  :5) + u(ktop): (5:3:2:3:20)

68
Here p(ktop) = p(ktop + :5)  p(ktop  :5). Finally, at the surface (the downdraft tops)

p(ksur)

g

( ~(ksur)) = (~(ksur  :5) + d(ksur)); (5:3:2:3:21)

with p(ksur) = p(ksur + :5)  p(ksur  :5). Here, the rst term is the compensatory

environmental mass 
ux, and the second term is the detrainment of all downdraft air at

the bottom. These normalized changes are also used in the calculation of the nal feedback

(after mb is determined), which is simply given by



@(k)

@ t 

CU = ((k))mb: (5:3:2:3:22)

To calculate the mass 
ux mb, we dene the buoyant energy which is available to a

cloud (updraft and downdraft) as

ABE =

k=Xktop

k=LF C"

g

cp T (k  :5)  

~h(kb)  ~h

(k  :5)

1 + 
(k  :5) !

 (z(k  1)  z(k))#

+

k

X=sur

k=z0

"

g

cp T (k  :5)  

h~(k0)  h~(k  :5)

1 + 
(k  :5) !

 (z(k)  z(k  1))# : (5:3:2:3:23)

where 
 is dened in equation (5.3.2.2.15). We can calculate ABE (similar to Lord 1982)

for the unchanged environment as well as for the environment which has been modied by

some arbitrary mass 
ux m0bt

0 . Hence, we can write

NA =



dABE

dt 

NCU =

ABE0  ABE

m0bt

0 : (5:3:2:3:24)

ABE are calculated using T0 and q0, while ABE0 are calculated after modication of

the thermodynamic variables by an arbitrary amount of mass 
ux, m0bt

0 , where

0(k) = (k) + ((k))m0bt

0: (5:3:2:3:25)

For a closure which depends on the rate of destabilization, we have to calculate the change

in the available buoyant energy due to large-scale or other subgrid-scale eects. We modify

the thermodynamic elds with

\(k) = 0 +



@

@ t 

LS+SUBG

t; (5:3:2:3:26)

69
where (5.3.2.3.26) is applied at every timestep t. These double prime quantities are then

used to calculate the changes in the available buoyant energy due to \non-convective"

eects. As a result, the equation for the mass 
ux becomes

mb =

ABE\  ABE

(ABE0  ABE)m0b

: (5:3:2:3:27)

5.3.3 Parameterization of shallow convection

The shallow convection scheme is constructed to be able to serve two tasks. It

parameterizes planetary boundary layer (PBL) forced shallow non-precipitating convection

as well as mid-tropospheric shallow convection caused by other sub-grid scale eects

(such as cloud top radiational cooling). The rst might not be necessary when the

parameterization is coupled to a higher order closure PBL scheme. It will transport

moisture from inside the boundary layer into the layers just above the boundary layer.

This is accomplished by emulating bubbles (forced by surface heat and moisture 
uxes

only, with strong lateral mixing) which rise without precipitation formation through the

top of the boundary layer into the free atmosphere, where they then lose their buoyancy.

Because of the strong lateral mixing, they usually do not rise more than 50-75 mb. The

physics involved in describing the second kind of shallow convection is the same, except

for the forcing.

To parameterize this type of convection we assume that a \convective element" can

be characterized by a bubble which rises through several model layers. It is assumed

to be forced by planetary boundary layer 
uxes or radiational cooling tendencies. Some

of the elements of this parameterization are based on an Arakawa-Schubert type scheme

(section 5.3.2.2) and some are based on the simple one-cloud scheme described in section

5.3.2.3. However, the clouds (shallow \convective elements\) are characterized by dierent

properties. They usually have large mixing, are non-precipitating, and have no convectivescale downdrafts. They are forced by subgrid-scale processes only. The following

description will be focused on dierences from the previously described models. Since

the sole purpose of this scheme is to represent \very" shallow convection, it is also

constructed as a one-cloud scheme. Although it implicitly uses equations (5.3.2.2.1)-

(5.3.2.2.4), considerable simplications can be made by assuming strong lateral mixing

70
(detrainment being equally as strong as entrainment). Equations (5.3.2.2.1) through

(5.3.2.2.4) then read

 = 0; (5:3:3:1)

e = d =

:2

r

; (5:3:3:2)

and

@c

@z =

:2

r

( ~  c) + Sc; (5:3:3:3)

where r in equation (5.3.3.2) is the radius of the element. The parameterization will be

sensitive to the choice of r. For this type of convection we assume r = 50m. When assuming

that no precipitation forms or evaporates, equations (5.3.3.1)-(5.3.3.3), together with initial

conditions (5.3.2.2.25) and (5.3.2.2.26), form a simple set of equations to determine the

properties of the convective element, if r is given. Without precipitation formation, Sc

in equation (5.3.3.3) is zero. For the feedback, equations (5.3.2.2.32)-(5.3.2.2.34) simply

become

F Ss(z)  [sc(z)  s(z)]mc; (5:3:3:4)

F Sq (z)  [qc(z)  q(z)]mc; (5:3:3:5)

F Sl(z)  l(z)mc = 0: (5:3:3:6)

The only unknown in these equations is the mass 
ux. It is determined in the dynamic

control, where we make use of the denition of the cloud work function (5.3.2.2.40) and

simply impose



dA(scl)

dt 

CU = 



dA(scl)

dt 

SUBG

: (5:3:3:7)

Note that since the cloud work function is independent of mass 
ux (mass 
ux is constant

with height), equation (5.3.2.2.40) for cloud-type scl simplies to

A(scl) = Z zT

zB

g

CpT (z)

1

1 + 


(hc(z)  h~(z))dz : (5:3:3:8)

Subscript CU refers to the eects due to convection, and SUBG to eects due to sub-grid

scale forcing. A(scl) becomes simply the buoyancy which is available for that particular

71
cloud scl. Therefore, physically, the change of the eciency of kinetic energy generation

due to cloud scl is directly proportional to the buoyancy generation by sub-grid scale

forcing. To arrive at a useful closure, the term on the left hand side of equation (5.3.3.7)

is normalized by the mass
ux to yield

mc

dA(scl)

dt 

NCU = 



dA(scl)

dt 

SUBG

: (5:3:3:10)

Here, the subscript NCU now stands for the change of A due to a unit mass of cloud scl.

The variables in equation (5.3.3.10) are known, except for mc. After using (5.3.3.10) to

calculate mc, we can then calculate the feedback. Note that in equation (5.3.3.2), mc is not

dependent on height, and is simply the cloud base mass 
ux. It should be noted here that

the above described parameterization will greatly benet from a high vertical resolution.

In some instances it may be of use to allow the shallow convection scheme to be called

several times in a column (stacked on top of each other), since dierent sub-grid-scale

forcing mechanisms may act at the same time in one column, but at dierent levels.

72
5.4 Planetary boundary layer parameterizations

5.4.1 Surface-Energy equation

Over land, the surface temperature Tg is computed from a surface energy budget that

is base on the \force-restore" method developed by Blackadar (Zhang and Anthes 1982).

The budget equation is

Cg

@Tg

@ t = Rn  Hm  Hs  LvEs; (5:4:1:1)

where Cg is the thermal capacity of the slab per unit area, Rn is the net radiation, Hm is

the heat 
ow into the substrate, Hs is the sensible heat 
ux into the atmosphere, Lv is the

latent heat of vaporization, and Es is the surface moisture 
ux. Blackadar (1979) shows

that the following formulation enables the amplitude and phase of the slab temperature to

be identical to the surface temperature of a real soil layer of uniform thermal conductivity

 and heat capacity per unit volume Cs, with Cg related to these parameters and the

angular velocity of the earth 
 by

Cg = :95

Cs

2
 1=2

: (5:4:1:2)

The thermal capacity, Cg , is related to a parameter called the thermal inertia, , where 

is

 = (Cs)1=2: (5:4:1:3)

From (5.4.1.2) and (5.4.1.3),

Cg = 3:293 106; (5:4:1:4)

where  (cal cm2K1s

1

2 ; 1 cal = 4:18 J ) is specied in the model as a function of land-use

characteristic (Appendix 4). The terms on the right hand side of (5.4.1.1) are described

as follows:

73
5.4.1.1 Net radiative 
ux Rn

Radiation is the driving force of the diabatic planetary boundary layer (PBL) and is

the most important component of the slab-energy budget.

Rn = Qs + Is (5:4:1:5)

where Qs and Is are the net surface shortwave and longwave irradiances.

a. Clear Sky

For clear sky, the amount of solar radiation absorbed by the slab, including multiple

re
ection of short waves, is approximated as

Qs = S0(1  A) cos ; (5:4:1:6)

where S0 is the solar constant (1395.6 W m2), A is the albedo. is the zenith angle,

and  is the short-wave transmissivity. The term cos is given by

cos = sinsin + coscoscosh0; (5:4:1:7)

where  represents the latitude of the location,  the solar declination, and h0 the local

hour angle of the sun (Sellers. 1974).

The short-wave transmissivity for multiple re
ection (Benjamin 1983) is

 =

a[s + (1  s)(1  b)]

(1  XRA)

; (5:4:1:8)

where a is the absorption transmissivity, s is the scattering transmissivity, b is the

backscattering coecient, and XR is the multiple re
ection factor

XR = ad(1  sd)bd; (5:4:1:9)

where the subscript d denotes diuse.

All the clear-air transmissivities (a; s; ad; sd) and backscattering coecients (b and

bd) are determined as a function of path length and precipitable water from a look-up table

from the Carlson and Boland (1978) radiative transfer model. Transmissivities are then

adjusted for surface pressure as follows:

 =

1+(

0  1)ps

1013:25

; (5:4:1:10)

74
where 

0 is the transmissivity from the look-up table (appendix 2) obtained by assuming

the surface pressure is 1013.25mb, and ps is the surface pressure at the location. The net

longwave radiation, Is, is equal to the sum of the outgoing (I ") and downward (I #)

longwave radiation. The outgoing longwave radiation is

I "= g SB T 4

g ; (5:4:1:11)

where g is the slab emissivity, Tg is the ground temperature, and SB the StefanBoltzmann constant. The downward longwave radiation absorbed at the surface is

I #= g aSBT 4a ; (5:4:1:12)

where Ta is the atmospheric temperature in the layer above the surface, and a, the

atmospheric longwave emissivity, is given by

a = :725 + :17log10wp; (5:4:1:13)

in which wp is the precipitable water in centimeters

b. Cloudy skies

For cloudy skies, a cloud parameterization scheme (Benjamin 1983) is used to simulate

the eects of clouds on short-wave and downward longwave radiation. Groups of sigma

levels are chosen to correspond to low-, middle-, and upper-cloud layers based upon an

assumed surface pressure of 1000mb. The clouds below 800mb are designated as low

clouds, middle clouds are those between 800mb and 450mb, and upper clouds are those

above 450mb.

The attenuation of short-wave radiation by cloud is parameterized with absorption

(ac) and scattering (sc) transmissivities. The transmissivities through the three cloud

layers are given by

ac = Y

3

i=1

[1  (1  ai)]ni (5:4:1:14)

and

sc = Y

3

i=1

[1  (1  si)]ni ; (5:4:1:15)

75
where i = 1; 2; 3 represents low, middle, and high clouds, respectively, ni is the cloud

fraction, and ai and si are given in table 5.1. The minimum short-wave absorption

transmissivity is set at 0.7, and the minimum scattering transmissivity is set at 0.44.

The cloud fraction is based on relative humidity. Cloud fraction at low and middle

levels is

n = 4:0RH  3:0; (5:4:1:16)

and in the upper atmosphere

n = 42:5RH  1:5; (5:4:1:17)

where RH is the maximum relative humidity found in the model layers within the low,

middle, or upper cloud layers. The expression for eective short-wave transmissivity under

cloudy skies is

 =

acsca[s + (1  s)(1  b)]

(1  XcA)

; (5:4:1:18)

where the multiple re
ection factor for cloudy skies (Xc) is dened as

Xc = adac(1  sdsc)bd (5:4:1:19)

in which bd, the mean backscattering coecient, is



bd =

bd(1  sd) + (1  sc)

(1  sd) + (1  sc)

: (5:4:1:20)

The cloud enhancement of long-wave radiation incident on the ground is expressed as

I #

0= I # 

1 +X

3

i=1

cini!; (5:4:1:21)

where ci are the enhancement coecients at dierent levels (table 5.2).

76
Table 5.1 Cloud absorption and scattering transmissivities.

Cloud Level Absorption

(τai)

Scattering

(τsi)

Low 0.80 0.48

Middle 0.85 0.60

High 0.98 0.80

Table 5.2 Enhancement coefficients ci on longwave radiation due to clouds.

Cloud Level Coefficient

Low 0.26

Middle 0.22

High 0.06

77
5.4.1.2 Heat Flow into the Substrate Hm

The transfer of heat due to molecular conduction is calculated from the equation

Hm = KmCg (Tg  Tm); (5:4:1:22)

where Km is the heat transfer coecient expressed as Km = 1:18
, 
 is the angular

velocity of the earth, and Tm is the temperature of the substrate, which is presently

taken to be a constant value equal to the mean surface-air temperature over the period of

simulation. If the model is used in a forecast mode rather than a research mode, Tm may

be set equal to the mean surface temperature of the previous day.

5.4.1.3 Sensible-Heat Flux Hs and Surface Moisture Flux Es

These 
uxes are computed in dierent ways, depending upon what PBL parameterization is used. Details will be described in the next sections.

5.4.2 Bulk-aerodynamic parameterization

The bulk-aerodynamic option of the PBL physics follows Deardor (1972). It is a

very inexpensive choice. The surface-heat 
uxes are given by

Hs = acpmCCu(g  a)V ; (5:4:2:1)

where a and a are density and potential temperature at the lowest model layer, C and

Cu are exchange coecients (Deardor 1972) dened as

Cu = CuN 

1  RiB

RiC 

(5:4:2:2)

and

C = CN 

1  RiB

RiC 

(5:4:2:3)

for stable conditions (0  RiB  :9Ric), and

Cu =

1

1

CuN  25exp(:26  :03 2

)

(5:4:2:4)

and

C =

1

1

CN + 1Cu  1CuN

(5:4:2:5)

78
for the unstable case (RiB  0). Here CuN and CN are the neutral values for Cu and C,

and are given by

CuN =



k1ln 

:025h

z0  + 8:41

(5:4:2:6)

and

CN =



0:74k1ln 

:025h

z0  + 7:31

; (5:4:2:7)

where Ric = 3:05, h is the depth of the lowest model layer, is dened as

 = log10(RiB)  3:5; (5:4:2:8)

and the velocity V is given by

V = (V 2a + V 2

c )1=2

: (5:4:2:9)

Va is the wind-speed at the lowest model layer, and Vc is a convective velocity, which is

important under conditions of low mean wind-speed and is dened under unstable and

neutral conditions as

Vc = 2(g  a)1=2; (5:4:2:10)

while it is zero under stable conditions.

The surface moisture 
ux is

Es = aCCuM(qvs (Tg )  qva)V ; (5:4:2:11)

where M is the moisture availability parameter which varies from 1.0 for a wet surface to

0.0 for a surface with no potential for evaporation. The moisture availability is specied as

a function of land-use category (Appendix 4). The model results are often quite sensitive

to the value used for M.

The surface momentum 
ux is given by

s = aCDV 2; (5:4:2:12)

where the drag coecient CD is dened as

CD = C0D + 3 103 

s s + 9800

: (5:4:2:13)

79
The second term in (5.4.2.13), involving the surface geopotential s, is a correction for

elevated terrain (Bleck,1977). The expression for C0D follows Deardor (1972), where

C0D = C2u (5:4:2:14)

5.4.3 Blackadar High-resolution model

A revised version of Blackadar's PBL model (Blackadar, 1976, 1979; Zhang and

Anthes, 1982) is used to forecast the vertical mixing of horizontal wind (u and v), potential

temperature (), mixing ratio (qv ), cloud water (qc), and ice (qi). The surface heat and

moisture 
uxes are computed from similarity theory. First the friction velocity, u, is

computed based on

u = MAX 

kV

ln za

z0  m

; u0!; (5:4:3:1)

where u0 is a background value (0.1ms1 over land and zero over water) and V is given

by (5.4.2.9). The surface-heat 
ux is computed from

Hs = CpmakuT; (5:4:3:2)

where

T =

a  g

ln za

z0  h

; (5:4:3:3)

where z0 is the roughness parameter, za is the height of the lowest -level, and m and

 h are nondimensional stability parameters that are a function of the bulk Richardson

number RiB, which is given by

RiB =

gza

a

va  vg

V 2

; (5:4:3:4)

where the subscript v represents virtual potential temperature. There are four cases

possible:

a. Stable case

For the stable case, RiB > Ric, where the critical Richardson number Ric is dened

as

Ric = :2: (5:4:3:5)

80
In this case,

u = u0; (5:4:3:6)

 m = h = 10ln za

z0

; (5:4:3:7)

and

Hs = M ax(250 W m2; cpmakuT): (5:4:3:8)

b. Mechanically driven turbulence

For this case 0  RiB  Ric, and we get

 m = h = 5  RiB

1:1  5RiB 

ln za

z0

: (5:4:3:9)

c. Unstable (forced convection)

Here RiB < 0 and j h=L j  1:5, where the Monin-Obukhov length, L, is dened as

L = 

cpmaau3



kgHs

(5:4:3:10)

and h is the height of the PBL. In this case, m = h = 0, and za=L = RiBln za

z0

.

d. Unstable (free convection)

Here RiB < 0 and j h=L j > 1:5. In this case

 h = 3:23 

zaL

  1:99 

zaL

2  0:474 

zaL

3

; (5:4:3:11)

and

 m = 1:86 

zaL

  1:07 

zaL

2  0:249 

zaL

3

: (5:4:3:12)

where za=L is restricted to be no less than -2.0 in this approximation. For za=L equal to

-2.0, h = 2:29, and m = 1:43.

In the general case, za=L is a function of m and (5.4.3.12) is an implicit equation

requiring an iterative solution. To save time, we approximate za=L as an explicit function

of RiB, such that

zaL = RiB ln za

z0

: (5:4:3:13)

81
The above scheme ensures continuity of m for all values of RiB. The formulation for the

surface moisture 
ux in the multi-layer case was derived from Carlson and Boland (1978),

where

Es = M aI1(qvs(Tg )  qva); (5:4:3:14)

and

I1 = ku ln 

kuza Ka +

za

zl   h1

: (5:4:3:15)

The quantity zl is the depth of the molecular layer (0.01 m over land and z0 over water)

and Ka is a background molecular diusivity equal to 2:4 105m2s1.

Over land, the roughness length z0 is specied as a function of land-use category

(Appendix 4). Over water, z0 is calculated as a function of friction velocity (Delsol et al.

, 1971) such that

z0 = 0:032u2

=g + z0c; (5:4:3:16)

where z0c is a background value of 104m.

The Blackadar scheme considers two dierent PBL regimes, the nocturnal regime and

the free-convection regime. The rst three cases (stable, mechanically driven turbulence,

and forced convection) are in the nocturnal regime, which is usually stable or at most

marginally unstable.

Nocturnal Regime

The rst-order closure approach is used to predict model variables. The ground stress

is calculated from

s = u

2

; (5:4:3:17)

where u is computed from (5.4.3.1). The components of s in the x and y directions are

sx =

u

Va

s (5:4:3:18)

and

sy =

v

Va

s; (5:4:3:19)

where Va is the wind speed at the lowest model level. For surface layer variables, the

prognostic equations are

@a

@ t =

 (H1  Hs)

(acpmz1)

; (5:4:3:20)

82
@qva

@ t =

 (E1  Es)

(az1)

; (5:4:3:21)

@ua

@ t =

(1x  sx)

(az1)

; (5:4:3:22)

@va

@ t =

(1y  sy)

(az1)

; (5:4:3:23)

and

@qca

@ t =

F1

(az1)

; (5:4:3:24)

where Hs is the surface heat 
ux computed from (5.4.3.2), ES is the surface moisture 
ux

computed from (5.4.3.14), subscript a refers to surface layer variables, subscript 1 refers to

the 
uxes at the top of the surface layer (Fig. 5.4), and z1 is the height of the lowest model

layer. The 
uxes at the full  levels are computed from K-theory, as described in section

(5.4.4). The prognostic variables above the surface layer are computed from K-theory and

an implicit diusion scheme (Richtmeyer, 1957; Zhang and Anthes, 1982).

Free-Convection Regime

During strong heating from below, large surface heat 
uxes and a super- adiabatic

layer occur in the lower troposphere. As the buoyant plumes of hot air rise under such

unstable conditions, mixing of heat, momentum, and moisture take place at each level. The

vertical mixing in this scheme is not determined by local gradients, but by the thermal

structure of the whole mixed layer. In the Blackadar PBL model, the vertical mixing

is visualized as taking place between the lowest layer and each layer in the mixed layer,

instead of between adjacent layers as in K-theory.

In the surface layer, the prognostic variables are solved by the analytic solution

+1 a = 1 a +



Fsz1

mh  2 

Fs mh  +F1

mh   



exp mh  t

z1   1 +

Fst

h

; (5:4:3:25)

where  represents any prognostic variable, Fs is the surface 
ux, F1 is the 
ux at the

top of the surface layer, h is the height of the PBL, t is the time-step, and the mixing

coecient is

m = H1

"

acpm(1  ) Z h

z1

[va  v (z

0 )] dz0#1

: (5:4:3:26)

83
Here  is the entrainment coecient (0.2) and H1 is the heat 
ux at the top of the surface

layer computed by the Priestly equation

H1 = acpmz1(va  v

1 1

2 )

3

2 

2g

27va 1

2 1

z1

h

z 1

3 1  (2z

1 1

2 ) 1

3 i 3

2

; (5:4:3:27)

where z1 is the depth of the surface layer and the subscript 1 1

2

refers to the second

prediction layer above the surface (Fig. 5.4).

For the variables above the surface layer, the prognostic equation is

@i

@ t

= m(a  i);  = ; qv ; or qc (5:4:3:28)

@i

@ t = wm (a  i);  = u; v: (5:4:3:29)

The variable w is a weighting function for reducing mixing near the top of the mixed layer,

where

w = 1 

z

h

: (5:4:3:30)

Care must be taken at the layer where the top of the mixed layer is located because the

top of the mixed layer does not necessarily coincide with a model level.

5.4.4 Vertical diusion

Above the mixed layer, K-theory is used to predict the vertical diusion of the

prognostic variables, such that

FV  = p

 @

@z Kz @@z ; (5:4:4:1)

where the eddy diusivity, Kz , is a function of the local Richardson number Ri . Specically,

Kz = Kz0 + l2 S

1

2 Ric  Ri

Ric

for Ri < Ric (5:4:4:2)

Kz = Kz0; for Ri  Ric (5:4:4:3)

where Kz0 = 1 m2 s1; l = 40 m ; andRic is a critical Richardson number which is a

function of layer thickness (m) and is dened as

Ric = :257 z

:175

: (5:4:4:4)

84
Z1

Z5

Z4 1/2

Z4

Z3

Z3 1/2

Z2

Z2 1/2

Z1 1/2

Za

N

P

F

F

F

F

F

Fs

α

α

α

α

α

θva θvg

θ

h

x

x

x

x

x

Figure 5.4 Illustration of vertical grid structure for high-resolution (Blackadar) model. The top of

the surface layer is z1; θvg and θva are the virtual potential temperatures of the ground surface and

lowest model level, respectively; P and N denote the positive and negative areas associated with a

parcel of air originating at za and rising to h, the top of the PBL.

85
According to (5.4.4.4), Ric varies from 0.58 for z = 100 m to 0.86 for z = 1000 m.

The Richardson number is

Ri =

g

S

@

@z (5:4:4:5)

and S is

S =



@u

@z 2

+



@v

@z 2

+ 109: (5:4:4:6)

5.4.5 Moist vertical diusion

There is an option with explicit moisture of including the eects of moisture on vertical

diusion. Taking into account moist-adiabatic processes in cloudy air (Durran and Klemp

1982), (5.4.4.5) is modied to

Ri = (1 + ) "

g

S

@

@z 

g2 

1+

ScpT

#

(5:4:5:1)

where

 =

L2

v qvs

cpRvT 2

(5:4:5:2)

and

 =

Lvqvs RdT

; (5:4:5:3)

and this modied value is used in (5.4.4.2) where the cloud amount exceeds 0.01 g

kg1.

86
5.5 Atmospheric radiation parameterization

The atmospheric radiation option in the model provides a longwave (infra-red) and

shortwave (visible) scheme that interact with the atmosphere, cloud and precipitation

elds, and with the surface (Dudhia 1989).

5.5.1 Longwave radiative scheme

Longwave absorption by water vapor, the primary clear-air absorber, is strongly

spectral in character, and the method employed is the commonly used broadband emissivity

method (see Stephens 1984). This involves using a precalculated emissivity function, ,

which represents the frequency-integrated absorption spectrum of water vapor, weighted

by a suitable envelope function. Rodgers (1967) gives an upward and downward emissivity

as a function of water vapor path, u, with a temperature correction term, where u includes

a pressure correction factor proportional to p0:86. The form of the tted function is

(u) = X

i=4

i=0

(ai + T bi)x

i

; (5:5:1:1)

where x = ln u and T is a u-weighted T - 250K. For u less than 10 g m2

, the form is

(u) = X

i=4

i=1

(ci + T di)y

i

; (5:5:1:2)

where y = u1=2 and ai , bi , ci and di are constants. In the tropics, e-type absorption is

an important additional component of the longwave absorption spectrum and is included

with a similar fourth-order polynomial in ln (ue) to (5.5.1.1) from Stephens and Webster

(1979), where e is the partial pressure due to water vapor. Given the emissivity functions

from (5.5.1.1-2) (u for upward 
ux and d for downward 
ux), the upward and downward


uxes at any model level are given by

Fu =

Z 1

0

B(T )du; (5:5:1:3a)

Fd =

Z 1

0

B(T )dd; (5:5:1:4a)

In (5.4.1.3a) the integration is performed downwards through the model layers. The

quantity d is calculated for each layer using the temperature (T ) of the layer and the

87
frequency-integrated Planck function B = SBT 4

, where SB is the Stefan-Boltzmann

constant. When the surface is encountered, the ground emission Fbot is multiplied by

1   and added to the integration. In (5.5.1.4a), the integration is performed upwards;

the downward longwave 
ux at the model top, Ftop, is assumed to result only from CO2

emission in the stratosphere. Thus (5.5.1.3a-4a) can be expressed as

Fu(z) = Z z

0=zsf c

z0=z B(T ) du

dz0 dz0 + Fbot[1  u(z; zsf c)]; (5:5:1:3b)

and

Fd(z) = Z z

0=ztop

z0=z B(T ) dd

dz0 dz0 + Ftop[(1  d(z; ztop)]; (5:5:1:4b)

where

(z; z1) = Z z1

z d

dz0 dz0 : (5:5:1:5)

It can be seen from the formulas that if the emissivity reaches 1 during the integration, the

remaining atmosphere makes no contribution to the 
ux. This is consistent with the idea

that an emissivity of 1 corresponds to a \black" layer with respect to longwave radiation.

Following Stephens (1978), the cloud water is assumed to have a constant absorption

coecient which is slightly dierent for upward and downward radiation. The absorption

coecients are cu = 0:130 m2 g1 and cd = 0:158 m2 g1. To combine these with water

vapor absorption, the transmissivities are multiplied since clouds are assumed to be \grey

bodies." The net emissivity is then

tot = 1  TvTc; (5:5:1:6)

with

Tv = 1  vaporand (5:5:1:7)

Tc = exp(cuc); (5:5:1:8)

where uc is the cloud water path (liquid mass per unit area).

Ice cloud is assumed to be composed of hexagonal plate-like crystals with the diametermass relation given in section (5.3.1.1). If the assumption is made that the crystals do not

re
ect longwave radiation and are suciently thick to be \black", it is possible to estimate

88
an absorption coecient as an integrated cross-sectional area. Allowing for the random

orientation of these crystals and a hemispheric integration factor of 1.66, the absorption

coecient takes a value of i = 0:0735 m2 g1, or about half that of cloud water. Since

this value agrees with some observations, it was applied in the model.

For rain and snow, the size distribution is necessary since the cross section is not

proportional to the mass of a particle. The size spectrum changes with precipitation

intensity so the absorption coecient varies with precipitation amount. The eective

absorption coecient is given by

p =

1:66

2000 

N0

3

r

1=4

m

2

g1; (5:5:1:9)

where r is the particle density. For the constants used in the explicit moisture scheme

described earlier, the absorption coecients take values of 2:34 103m2g1for snow and

0:330 103m2g1for rain. The eective water path for a layer of z meters thickness is

given by

up = (qr )3=4z 1000gm2; (5:5:1:10)

so that the transmissivity is given by

Tp = exp(pup): (5:5:1:11)

This transmissivity is multiplied with the others in (5.5.1.6) to give tot . This is known as

an overlap approximation. Rain and snow have less eect on the longwave 
ux by 2 to 3

orders of magnitude, but still are not insignicant.

Carbon dioxide is less easily treated since it cannot be assumed \grey". That is,

its absorption is concentrated in a band of infrared wavelengths. To include its eect,

an overlap method is used as discussed by Stephens (1984). In eect, the spectrum is

divided into a carbon-dioxide band and a non-carbon-dioxide region. The former requires

overlapping of the carbon dioxide transmissivity function while the latter does not. The

relative weights of these two regions is slightly temperature dependent, but they add to

give the total absorption. A pressure correction factor proportional to p1:75 is applied to

the carbon dioxide path calculation.

89
Having obtained the 
ux proles, Fu(z) and Fd(z), the radiative heating rate is

calculated from

Q_

R = cp

@T

@ t =

1



@

@z (Fd  Fu) = g

@

@p (Fd  Fu): (5:5:1:12)

In the model, the values of F are dened on the model full sigma-levels. This makes

the various integrals and derivatives easier to represent numerically.

5.5.2 Shortwave Radiative Scheme

The downward component of shortwave 
ux is evaluated taking into account 1) the

eects of solar zenith angle, which in
uences the downward component and the path length;

2) clouds, which have an albedo and absorption; 3) and clear air, where there is scattering

and water-vapor absorption. Thus,

Sd(z) = S0 

Z top

z

(dScs + dSca + dSs + dSa); (5:5:2:13)

where  is the cosine of the zenith angle and S0 is the solar constant.

As with the longwave scheme, cloud fraction in a grid box is either 0 or 1 because of

the assumed stratiform nature of the clouds. The cloud back-scattering (or albedo) and

absorption are bilinearly interpolated from tabulated functions of  and ln(w=) (where

w is the vertically integrated cloud water path) derived from Stephens' (1978) theoretical

results. The total eect of a cloud or multiple layers of cloud above a height z is found

from the above function as a percentage of the downward solar 
ux absorbed or re
ected.

Then at a height z  z, a new total percentage is calculated from the table allowing

the eect of the layer z to be estimated. However, this percentage is only applied to

S0  S(clear air); that is, the clear-air eect above z is removed.

Clear-air water vapor absorption is calculated as a function of water vapor path

allowing for solar zenith angle. The absorption function is from Lacis and Hansen (1974).

The method is a similar integration-dierence scheme to that described above for cloud.

Clear-air scattering is taken to be uniform and proportional to the atmosphere's

mass path length, again allowing for the zenith angle, with a constant giving 20 percent

scattering in one atmosphere. The heating rate is then given by

RT = RT (longwave) + 1

cp

Sabs; (5:5:2:14)

90
where Sabs is dened from the absorption part of the Sd integral given in (5.5.2.13), since

only cloud and clear-air absorption contribute to solar heating.

The solar and infrared 
uxes at the surface, calculated from the atmospheric radiative

schemes, are use in the energy budget of the land surface.

91
Appendix 1. Glossary of Symbols

a Fraction of convective cloud cover; also constant used in cloud

microphysics

ABE Available buoyant energy

AT The forcing terms of the thermodynamic equation that vary on the

time-scale of the Rossby-waves

Au; Ad; Atot Cloud Work Function for updraft, downdraft, and all of model cloud

Av The forcing terms of the v-momentum equation that vary on the timescale of the Rossby-waves

Av The forcing terms of the u-momentum equation that vary on the timescale of the Rossby-waves

A0 Parameter for heterogeneous freezing (K1)

A Antidiusive 
ux

b Backscattering coecient; also fraction of total water vapor convergence used to moisten grid column (section 5.3.2.1); also constant

(0.8) used in cloud microphysics computation

B Planck function

Bu Acceleration due to buoyancy

B0 Parameter for heterogeneous freezing (m3s1)

c0 Rainfall conversion parameter (section 5.3.2.2)

93
cm Constant used in computation of Dm

ci Coecients used in calculation of cloud eect on downward longwave

radiation (Table 5.2)

cp Specic heat at constant pressure for dry air

cpm Specic heat at constant pressure for moist air

c

 Net condensation rate averaged over grid volume

c



c Net condensation rate in cumulus cloud (section 5.3.2.1)

C Constant (2. m s1 K1=2) used in computing convective velocity

Cg Thermal capacity of slab per unit area (J m2 K1)

Cs Heat capacity per unit volume (J m3 K1)

C Surface exchange coecient for heat

Cu Surface exchange coecient for momentum; also total condensate in

updraft (section 5.3.2.2)

CD Surface drag coecient

C0D Component of surface-drag coecient

CuN Value of surface momentum exchange coecient under neutral

stability conditions

94
CN Value of surface heat exchange coecient under neutral stability

conditions

D Mass divergence (hydrostatic split-explicit scheme); also horizontal

deformation (section 5.1)

Df Diusivity of water vapor in air

Du; Dd; Dtot Updraft, downdraft, and total cloud kinetic energy dissipation

Distance between an observation and a given grid point (section 4)

Di Diameter of ice crystal (m)

D Diusion and PBL tendencies for variable 

Dm Modied distance between an observation and a given grid point

(section 4)

e Horizontal Coriolis parameter (s1)

es; esi ; esw Saturation vapor pressure, over ice, over water (cb)

E Eciency of collection of cloud by precipitation; also vertical 
ux of

water vapor

Es Flux of water vapor from surface into atmosphere

f Coriolis parameter

f1; f2 Ventilation coecients for rain or snow

95
F Larger-scale forcing (section 5.3.2.2 and 5.3.2.3); also function of

distance from lateral boundary (section 2.6.2)

FH; F L Flux from high-order and low-order advective scheme

Fbot; Ftop Longwave radiative 
ux at bottom, top of model atmosphere (W m2)

Fd; Fu Downward, upward longwave radiative 
ux (W m2)

FH Term representing contribution of horizontal diusion of a variable 

to the temporal rate of change of 

FV  Term representing contribution of vertical diusion of a variable  to

the temporal rate of change of 

Fs Flux of dry static energy (section 5.3.2.2); also Surface 
ux of heat,

moisture or momentum

Fq Flux of water vapor (section 5.3.2.2)

Fl Flux of suspended cloud liquid water (section 5.3.2.2)

Fr Deceleration

F1; F2 Amplitude factors used in computing lateral boundary conditions

(section 2.6.2)

g Acceleration of gravity (9.8 m s2)

h Moist static energy; also height of planetary boundary layer (m)

ho Local hour angle of sun

96
h; h ~

c; hu; hd; ~h Moist static energy in environment, cloud, updraft, downdraft, and

saturation value in environment

H Vertical 
ux of sensible heat (W m2)

Hm Heat 
ux into substrate (W m2)

Hs Sensible heat 
ux from surface into atmosphere (W m2)

I Function of static stability and surface friction velocity; also horizontal

grid-index in y-direction

IMAX Maximum value of grid-index in y-direction

Is Net longwave iradiance at surface (wm2)

I " Outgoing longwave radiation from surface (W m2)

I # Downward longwave radiation absorbed at surface (W m2) under

clear skies

I #0 Downward longwave radiation absorbed at surface (W m2) in

presence of clouds

I1 Normalized condensate in updraft (section 5.3.2.2)

I2 Normalized evaporate in updraft (section 5.3.2.2)

J Horizontal grid index in x-direction

JMAX Maximum value of grid index in x-direction

97
k Dimensionless x-wavenumber for upper radiative scheme; also von

Karman constant (0.4)

^k Dimensionless eective x-wavenumber for upper radiative scheme

k1 Constant used in formula for computing autoconversion of cloud drops

to rain drops

K Total horizontal wavenmumber (m1) also Kernels

Ka Background molecular diusivity (2:4 105m2s1); also thermal

conductivity of air ( J m1s1 K1)

KH Horizontal eddy diusivity (m2s1)

K0H Coecient used in fourth-order diusion (s1)

KHO Background value of horizontal eddy diusivity (m2s1)

Km Coecient of heat transfer from ground into substrate (s1)

KMAX Maximum value of index in vertical direction

Kz Coecient of vertical diusivity (m2s1)

Kz0 Background value of coecient of vertical diusivity (m2s1)

KEu;KEd;KEtot Kinetic energy for updraft, downdraft, and all of model cloud

L Hydrostatic term due to liquid water loading; also Monin-Obukhov

length

Lm Latent heat of fusion (0.35 106 J kg1)

98
Ls Latent heat of sublimation (2.85 106 J kg1

)

Lv Latent heat of condensation (2.5 106 J kg1)

l Dimensionless y-wavenumber for upper radiative scheme; also vertical

mixing length

^

l Dimensionless eective y-wavenumber for upper radiative scheme

Mi Mass of ice crystal (kg)

Mmax Maximum mass of ice crystal (kg)

M0 Initial mass of ice crystal (kg)

m Mass 
ux (updraft and downdraft) in convective parameterization

cloud (5.3.2.2); also map scale factor

m Mixing coecient used in free-convective regime of high-resolution

PBL model

mb Cloud base mass 
ux (section 5.3.2.2)

m0 Downdraft base mass 
ux (section 5.3.2.2)

mu Updraft mass 
ux in convective parameterization cloud (5.3.2.2)

md Downdraft mass 
ux in convective parameterization cloud (5.3.2.2)

M Surface moisture availability

99
Mt Vertical integral of horizontal convergence of water vapor

n Fraction of cloud

n0 Cloud microphysics parameter

nc Number concentration of ice crystals (kg1)

N Brunt-Vaisala frequency (s1)

Nc Number concentration of cloud droplets per unit volume (1010 m3)

Nh Nondimensional function for vertical prole of convective heating

Nm Nondimensional function for vertical prole of convective moistening

N0 Cloud microphysics parameter (8 106 m4for rain 2 107 m4for

snow)

p Pressure (cb)

pb Pressure (cb) at convective cloud base

ps Surface pressure (cb)

pt Pressure (cb) at top of model

pu Pressure (cb) at top of convective cloud

pLCL Pressure (cb) at lifting condensation level

100
p ps  pt (cb)

pd Dot-point p(cb)

p0 Reference-state pressure

p0 Perturbation pressure (Pa)

Pressure value representing the free atmosphere, where terrain

in
uences are small (in FDDA)

p^ Fourier transform of p

0 for upper radiative boundary condition

PCON Condensation of water vapor or evaporation of cloud drops (kg kg1

s1)

PRA Accretion of cloud drops by rain drops (kg kg1s1)

PRC Autoconversion of cloud drops to rain drops (kg kg1s1)

PRE Evaporation of rain drops (kg kg1s1)

PCI Heterogeneous freezing of cloud water (kg kg s1)

PID Deposition of vapor onto ice crystals (kg kg s1)

PI I Initiation of ice crystals (kg kg s1)

PMF Melting/freezing of cloud and precipitation due to advection (kg kg

s1)

101
PRM Melting of falling precipitation (kg kg s1

)

PSM Melting of falling snow (kg kg s1)

q; q ~ u; qd; q~ Water vapor mixing ratio in environment, updraft, downdraft, and

saturation value in environment

qc Mixing ratio of cloud water; also water vapor mixing ratio in cloud

(section 5.3.2.2)

qc0 Critical value of mixing ratio of cloud water

qr Mixing ratio of rain water

ql Suspended liquid water vapor mixing ratio inside updraft

qv Mixing ratio of water vapor

qvc Mixing ratio of water vapor in cumulus cloud

qvs Saturation mixing ratio of water vapor

Q Diabatic heating rate per unit mass (J kg1s1)

Qs Net short wave irradiance at the surface (W m2)

R Rainfall (convective-scale sink of cloud water, 5.3.2.2); also ideal gas

constant for dry air (287 J kg1 K1)

RH Relative humidity

Ri Richardson number

102
Rn Net radiation

RiB Bulk Richardson number

Ric Critical value of bulk Richardson number; also critical value of

Richardson number

Rv Gas constant for water vapor (461.5 J kg1 K1)

RT Radiative heating rate (K s1)

r Radius of convective parameterization cloud (sections 5.3.2.2)

S Supersaturation; also source or sink term (section 5.3.2.2); also square

of the vertical wind shear

Sc Schmidt number

So Solar constant (1395.6 W m2)

Su Source or sink term in updraft (section 5.3.2.2)

Sd Downward solar 
ux (W m2); also source or sink term in downdraft

(section 5.3.2.2)

Si Supersaturation over ice

s Dry static energy

t Time (s)

103
T Temperature (K)

Tc Longwave transmissivity due to cloud

Td Dewpoint temperature (K)

Tg Temperature (K) of ground

Tp Longwave transmissivity due to precipitation

Tv Virtual temperature (K); also longwave transmissivity due to vapor

T Surface friction temperature (K)

T0 Reference-state temperature (K)

T 0 Perturbation temperature (K)

u Component of wind velocity in eastward direction (m s1); also water

vapor path (g m2)

u Surface friction velocity (m s1)

uc; up Liquid water path for cloud, precipitation (g m2)

v Component of wind velocity in northward direction (m s1)

vt Mass weighted mean terminal velocity of rain drops (m s1)

V Fall speed of a precipitation particle (m s1); also modied horizontal

wind velocity in PBL

104
V Horizontal wind vector

Va Horizontal windspeed at lowest model layer

Vc Convective PBL velocity(m s1)

Vqf Divergence of vertical eddy 
ux of water vapor due to convective

clouds

w Vertical velocity (m s1); also weight function for reducing mixing

near top of mixed layer

wn Weight function for blending model tendencies and large-scale

tendencies near lateral boundaries (section 2.6.1)

wp Precipitable water (cm)

wu Vertical velocity in updraft

w^ Fourier transform of w

x Horizontal grid coordinate increasing generally eastward

X Horizontal coordinate on earth surface increasing generally eastward

Xc Multiple-re
ection factor in cloudy air

Xd Distance vector

XR Multiple-re
ection factor in clear air

105
y Horizontal grid coordinate increasing generally northward

Y Horizontal coordinate on earth surface increasing generally northward

z Height above surface (m)

za Height of lowest layer in model (m)

zb Height of updraft originating level(section 5.3.2.2) (m)

z0 Height of downdraft originating level(section 5.3.2.2); also surface

roughness length(m)

zoc Background value of surface roughness length over water (104m)

zl Depth of molecular layer

zLCL Height of lifting condensation level (m)

zT Height of updraft top (section 5.3.2.2)

 Coecient array for upper radiative boundary condition (m s1

Pa1); also any thermodynamic variable (section 5.3.2.2)

~ Any thermodynamic variable in environment

u Any thermodynamic variable in updraft

d Any thermodynamic variable in downdraft

c; p Longwave absorption coecients for cloud, precipitation (m2 g1)

106
 Parameter in sound-wave temporal dierencing; also precipitation

eciency parameter in section 5.3.2.2

 Gamma function

d Dry adiabatic lapse rate (K m1)

dp Dewpoint adiabatic lapse rate (K m1)


 Ratio of heat capacities (cp=cv) for dry air

 Solar declination

M Supersaturation or undersaturation

p Vertical grid size (Pa)

s Horizontal grid length (m)

t Time step (s)

t

0 Short time step for rain fall term (s)

x Horizontal grid length (m)

z Thickness of vertical layer (m)

 Thickness of model  levels

107
c Critical value of convective cloud depth

 Short time step for sound waves (s)

r2 Horizontal Laplacian on -surfaces

r4 Fourth order diusion operator on -surfaces

 Parameter relating updraft and downdraft mass 
ux (section 5.3.2.2);

also small value; also entrainment coecient used in high resolution

PBL-model (0.2)

a Atmospheric emissivity

g Emissivity of ground

u; d Atmospheric longwave emissivity

d Normalized mass 
ux for downdraft (section 5.3.2.2)

u Normalized mass 
ux for updraft (section 5.3.2.2)

 Potential temperature (K); also angle between y-axis and north for

full Coriolis force

a Potential temperature (K) at lowest layer in model

g Potential temperature (K) of ground surface

e Equivalent potential temperature (K)

108
es Saturation equivalent potential temperature (K)

v Virtual potential temperature (K)

 Longitude; also cloud type (section 5.3.2.2); also thermal conductivity

(J m1s1 K1); also parameter in raindrop distribution (m1)

 Dynamic viscosity of air (kg m1s1); also solar zenith angle ; also

total net fractional entrainment rate (section 5.3.2.2); also constant

in smoother (section 3.3)

u Total net fractional entrainment rate for updraft (section 5.3.2.2)

ue Gross fractional entrainment rate for updraft (section 5.3.2.2)

ud Gross fractional detrainment rate for updraft (section 5.3.2.2)

 Coecient for Asselin time lter; also for spatial smoother

 Exner function

 Density of air (kg m3)

r Particle density (kg m3)

u Density in updraft

w Density of water (kg m3)

0 Reference-state density (kg m3)

109
0 Perturbation density (kg m3

)

 Nondimensional vertical coordinate of model

0 Dummy variable of integration

_ Vertical velocity in -coordinates (s1)

_ c Vertical velocity of convective cloud in -coordinates (s1)

SB Stefan-Boltzmann constant (5:67051 108 J m2 K4s1)

 Half-period of time window of in
uence of an observation (section 4);

also short-wave transmissivity



0 Short-wave transmissivity obtained from lookup table

a Clear air absorption transmissivity

ac Cloudy air absorption transmissivity

s Clear air scattering transmissivity; also surface stress

sc Cloudy air scattering transmissivity

 Geopotential; also latitude; also scalar variable in advection equation

s Surface geopotential

110
 Symbol denoting low-order, monotonic solution to advection equation

 Diusivity of vapor in air (m2s1); also thermal inertia

	 Solar zenith angle; also function of bulk Richardson number

	m Nondimensional stability parameter for momentum

	h Nondimensional stability parameter for heat and water vapor

! Vertical velocity in pressure coordinates (cb s1)

!c Vertical velocity of convective cloud in pressure coordinates (cb s1)


 Angular velocity of earth (7:2722 105s1)

111
Appendix 2. Look-up table for transmissivities

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0

1.0 0.926 0.868 0.855 0.846 0.838 0.832 0.827 0.822 0.818 0.814 0.811

1.2 0.915 0.854 0.840 0.831 0.823 0.817 0.811 0.806 0.802 0.798 0.794

1.4 0.903 0.841 0.826 0.816 0.808 0.802 0.796 0.791 0.787 0.782 0.779

1.6 0.892 0.828 0.813 0.803 0.795 0.788 0.782 0.777 0.772 0.768 0.764

1.8 0.881 0.815 0.800 0.790 0.782 0.775 0.769 0.763 0.758 0.754 0.750

2.0 0.870 0.803 0.788 0.777 0.769 0.762 0.756 0.750 0.745 0.741 0.737

2.2 0.860 0.792 0.776 0.765 0.757 0.750 0.743 0.738 0.733 0.728 0.724

2.4 0.850 0.781 0.765 0.754 0.745 0.738 0.731 0.726 0.7210 0.716 0.712

2.6 0.839 0.770 0.753 0.742 0.733 0.726 0.720 0.714 0.709 0.704 0.700

2.8 0.830 0.759 0.743 0.731 0.722 0.715 0.709 0.703 0.698 0.693 0.689

3.0 0.820 0.748 0.732 0.721 0.712 0.704 0.698 0.692 0.687 0.682 0.678

3.2 0.810 0.738 0.722 0.710 0.701 0.694 0.687 0.681 0.676 0.671 0.667

3.4 0.801 0.728 0.712 0.700 0.691 0.683 0.677 0.671 0.666 0.661 0.656

3.6 0.791 0.719 0.702 0.690 0.681 0.674 0.667 0.661 0.656 0.651 0.646

3.8 0.782 0.709 0.692 0.681 0.671 0.664 0.657 0.651 0.646 0.641 0.636

4.0 0.773 0.700 0.683 0.671 0.662 0.654 0.648 0.642 0.636 0.631 0.627

4.2 0.764 0.691 0.674 0.662 0.653 0.645 0.638 0.632 0.627 0.622 0.618

4.4 0.756 0.682 0.665 0.653 0.644 0.636 0.629 0.623 0.618 0.613 0.608

4.6 0.747 0.673 0.656 0.644 0.635 0.627 0.621 0.615 0.609 0.604 0.600

4.8 0.738 0.665 0.647 0.636 0.626 0.619 0.612 0.606 0.600 0.596 0.591

5.0 0.730 0.656 0.639 0.627 0.618 0.610 0.603 0.597 0.592 0.587 0.582

5.2 0.722 0.648 0.631 0.619 0.610 0.602 0.595 0.589 0.584 0.5790 0.574

5.4 0.714 0.640 0.623 0.611 0.602 0.594 0.587 0.581 0.576 0.571 0.566

5.6 0.706 0.632 0.615 0.603 0.594 0.586 0.579 0.573 0.568 0.563 0.558

5.8 0.698 0.624 0.607 0.595 0.586 0.578 0.571 0.565 0.560 0.555 0.550

6.0 0.690 0.616 0.599 0.588 0.578 0.571 0.564 0.558 0.552 0.547 0.543

6.2 0.683 0.609 0.592 0.580 0.571 0.563 0.556 0.550 0.545 0.540 0.535

6.4 0.675 0.602 0.585 0.573 0.564 0.556 0.549 0.543 0.538 0.533 0.528

6.6 0.668 0.594 0.577 0.566 0.556 0.549 0.542 0.536 0.531 0.526 0.521

6.8 0.661 0.587 0.570 0.559 0.549 0.542 0.535 0.529 0.524 0.519 0.5147

7.0 0.653 0.580 0.563 0.552 0.542 0.535 0.528 0.522 0.517 0.512 0.507

7.2 0.646 0.573 0.556 0.545 0.536 0.528 0.521 0.515 0.510 0.505 0.501

7.4 0.639 0.567 0.550 0.538 0.529 0.521 0.515 0.509 0.503 0.499 0.494

7.6 0.633 0.560 0.543 0.532 0.522 0.515 0.508 0.502 0.497 0.492 0.488

7.8 0.626 0.553 0.537 0.525 0.516 0.508 0.502 0.496 0.491 0.486 0.481

8.0 0.619 0.547 0.530 0.519 0.510 0.502 0.496 0.490 0.484 0.480 0.475

8.2 0.613 0.541 0.524 0.512 0.503 0.496 0.489 0.484 0.478 0.473 0.469

8.4 0.606 0.534 0.518 0.506 0.497 0.490 0.483 0.477 0.472 0.467 0.463

8.6 0.600 0.528 0.512 0.500 0.491 0.484 0.477 0.472 0.466 0.462 0.457

8.8 0.594 0.522 0.506 0.494 0.486 0.478 0.472 0.466 0.461 0.456 0.451

9.0 0.587 0.516 0.500 0.489 0.480 0.472 0.466 0.460 0.455 0.450 0.446

9.2 0.581 0.511 0.494 0.483 0.474 0.467 0.460 0.454 0.449 0.445 0.440

9.4 0.575 0.505 0.488 0.477 0.468 0.461 0.455 0.449 0.444 0.439 0.435

9.6 0.569 0.499 0.483 0.472 0.463 0.456 0.449 0.444 0.438 0.434 0.429

9.8 0.563 0.494 0.477 0.466 0.458 0.450 0.444 0.438 0.433 0.428 0.424

10.0 0.558 0.488 0.472 0.461 0.452 0.445 0.439 0.433 0.428 0.423 0.419

precipitable water (cm)

Path length

112
Appendix 3. Map Pro jections

Map pro jections are constructed by pro jecting the surface of the earth onto a right

circular cone, cutting the cone, and 
attening it into a plane surface. Three pro jections are

available for the MM4 system { Polar stereographic, Lambert conformal, and Mercator.

Polar stereographic is preferred for high-latitude studies, Lambert conformal for middlelatitude studies, and Mercator for low-latitude studies. This appendix summarizes the

map scale factors for each pro jection and gives the equations for converting from latitude

and longitude to the x and y positions on the model grid.

Although the grid size x = y = s is constant on the model's grid, the actual

distance represented by s on the spherical earth varies with location on the grid because

the earth is curved. The map scale factor m is dened as the ratio of the distance on the

grid to the corresponding distance on the earth's surface

m =

distance on grid

actual distance on earth A:1

a. Lambert Conformal

Conformal means that the scale is equal in all directions about a point, so that shapes

of geographic features on the earth are preserved. The Lambert conformal grid is true at

latitudes 30 and 60N so that m = 1: at these latitudes. In general,

m =

sin 1

sin 

tan=2

tan 1=2

0:716

; A:2

where 1 = 30 and is the colatitude ( = 90  ).

It is sometimes necessary to compute the position (x; y) on the grid given the latitude

and longitude of a point, or vice versa. The following relations pertain to an X  Y grid

with center X = 0; Y = 0 at latitude 0 and longitude 0. Note that the relationship

between (x; y) and (X; Y ) is

x = X +

JMAX  1

2

s; A:3

y = Y +

IMAX  1

2

s; A:4

 = any longitude

113
0 = longitude of Y axis

 = any latitude

0 = latitude along 0 at which Y = 0

 = 90  

n = :716

 1 = 30

 0 = 90  0

a = 6370km

r =

a

n

sin 1 

tan =2

tan 1=2

n

; A:5

C2 =

a

n

sin 1 

tan 0=2

tan 1=2

n

; A:6

C1 = 0  90=n; A:7

0 = n( + C1); A:8

X = rcos0 ; A:9

Y = rsin0 + C2: A:10

The inverse problem to calculate latitude and longitude is done as follows:

0 = arctan 

Y  C2

X



; A:11

 =

0

n  C1; A:12

r = X

cos0 or

Y  C2

sin0 A:13

 = 2arctan "tan 1=2 

nr

asin 1 1=n#

; A:14

 = 90  : A:15

b. Polar Stereographic

114
For the polar stereographic pro jection, true at latitude 1 = 60N, the map scale

factor is

m =

1 + sin1

1 + sin A:16

The relationships between latitude and longitude and X and Y on the polar stereographic

grid are calculated as before on the Lambert conformal grid, but now n = 1.

r = amcos; A:17

C2 = asin 1 

tan 0=2

tan 1=2



; A:18

C1 = 0  90; A:19

0 =  + C1; A:20

X = rcos0 ; A:21

Y = rsin0 + C2: A:22

and for the inverse problem

0 = arctan 

Y  C2

X



; A:23

 = 0  C1; A:24

r = X

cos0 or

Y  C2

sin0 A:25

 = 2arctan tan 1=2 

r

asin 1  ; A:26

 = 90  : A:27

Note that the signs of Y  C2 and X in (A.23) must be considered to obtain the

correct quadrant for 0 .

c. Mercator

For the Mercator grid, 0(Y = 0) corresponds to the equator and the relationships

between X and Y and  and  are relatively simple

X = (acos1)(  0); A:28

115
Y = (acos1)ln 

1 + sin

cos 

; A:29a

Y = (acos1)ln[tan(45 + =2)]: A:29b

Note that (  0) in (A.28) must be expressed in radians. The latitude 1 at which the

Mercator pro jection is true is often taken to be 30.

The reverse problem, to obtain X and Y from  and , is also simple

 = 0 +

X

acos1

: A:30

To solve for , use (A.29b)

 = 90 + 2arctan exp 

Y

acos1  : A:31

116
Appendix 4. Land Use Categories

Description of land-use categories and physical parameters for summer (15 April-15

October) and winter (15 October-15 April).

Landuse

Integer

Identification

Landuse

Description

Albedo(%) Moisture

Avail. (%)

Emissivity

(% at 9 μ m)

Roughness

Length (cm)

Thermal Inertia

(cal cm-2 k-1 s-1/2)

Sum Win Sum Win Sum Win Sum Win Sum Win

1 Urban land 18 18 5 10 88 88 50 50 0.03 0.03

2 Agriculture 17 23 30 60 92 92 15 5 0.04 0.04

3 Range-grassland 19 23 15 30 92 92 12 10 0.03 0.04

4 Deciduous forest 16 17 30 60 93 93 50 50 0.04 0.05

5 Coniferous forest 12 12 30 60 95 95 50 50 0.04 0.05

6 Mixed forest and

wet land 14 14 35 70 95 95 40 40 0.05 0.06

7 Water 8 8 100 100 98 98 .01 .01 0.06 0.06

8 Marsh or wet land 14 14 50 75 95 95 20 20 0.06 0.06

9 Desert 25 25 2 5 85 85 10 10 0.02 0.02

10 Tundra 15 70 50 90 92 92 10 10 0.05 0.05

11 Permanent ice 55 70 95 95 95 95 5 5 0.05 0.05

12 Tropical or sub

tropical forest 12 12 50 50 95 95 50 50 0.05 0.05

13 Savannah 20 20 15 15 92 92 15 15 0.03 0.03

117
References

Anthes, R. A., 1972:The development of asymmetries in a three-dimensional numerical

model of the tropical cyclone. Mon. Wea. Rev., 100, 461-476.

, 1977: A cumulus parameterization scheme utilizing a one-dimensional cloud model.

Mon. Wea. Rev., 105, 270-286.

, and T. T. Warner, 1978: Development of hydrodynamic models suitable for air

pollution and other mesometeorological studies. Mon. Wea. Rev., 106, 1045-1078.

, E.-Y. Hsie, and Y.H. Kuo, 1987: Description of the Penn State/NCAR Mesoscale

Model Version 4 (MM4). NCAR/TN-282+STR, National Center for Atmospheric

Research, Boulder, CO, 66 pp.

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

large scale environment. Part I. J. Atmos. Sci., 31, 674-701.

, and V. R. Lamb, 1977: Computational design of the basic dynamical process of the

UCLA general circulation model. Methods in Computational Physics, 17, 173-265.

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

Benjamin, S. G., 1983: Some eects of surface heating and topography on the

regional severe storm environment. Ph.D. thesis, Department of Meteorology, The

Pennsylvania State University, 265 pp.

Blackadar, A. K., 1976: Modeling the nocturnal boundary layer. Preprints of Third

Symposium on Atmospheric Turbulence and Air Quality, Raleigh, NC, 19-22 October

1976, Amer. Meteor. Soc., Boston, 46-49.

, 1979: High resolution models of the planetary boundary layer. Advances in

Environmental Science and Engineering, 1, Pfain and Ziegler, Eds., Gordon and

Breach Publ. Group, Newark, 50-85.

Betts, A. K., 1974: The scientic basis and ob jectives of the U.S. subprogram for the

GATE. Bul l. Amer. Meteor. Soc., 55, 304-313.

117
Bleck, R., 1977: Numerical simulation of lee cyclogenesis in the Gulf of Genoa. Mon. Wea.

Rev., 105, 428-445.

Bougeault, P., 1983: A non-re
ective upper boundary condition for limited-height

hydrostatic models. Mon. Wea. Rev., 111, 420{429.

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

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

Carlson, T. N., and F. E. Boland, 1978: Analysis of urban-rural canopy using a surface

heat 
ux/temperature model. J. Appl. Meteor., 17, 998-1013.

Davies, H. C., and R. E. Turner, 1977: Updating prediction models by dynamical

relaxation: An examination of the technique. Quart. J. Roy. Meteor. Soc., 103, 225-245.

Deardor, J. W., 1972: Parameterization of the planetary boundary layer for use in general

circulation models. Mon. Wea. Rev., 100, 93-106.

, 1978: Ecient prediction of ground surface temperature and moisture, with inclusion

of a layer of vegetation. J. Geophys. Res., 83, 1889-1903.

Delsol, F., K. Miyakoda, and R. H. Clarke, 1971: Parameterized processes in the surface

boundary layer of an atmospheric circulation model. Quart. J. Roy. Meteor. Soc., 97, 181-208.

Dudhia, J., 1989: Numerical study of convection observed during the winter monsoon

experiment using a mesoscale two-dimensional model. J. Atmos. Sci., 46, 3077-3107.

, 1993: A nonhydrostatic version of the Penn State / NCAR mesoscale model: Validation

tests and simulation of an Atlantic cyclone and cold front. Mon. Wea. Rev., 121,

1493{1513.

Durran, D.R., and J.B. Klemp, 1982: On the eects of moisture on the Brunt-Vaisala

frequency. J. Atmos. Sci., 39, 2152{2158.

Fritsch, J. M., and C. F. Chappel, 1980: Numerical prediction of convectively driven

118
mesoscale pressure systems. Part I: Convective parameterization. J. Atmos. Sci., 37, 1722-1733.

Garratt, J.R., J.C. Wyngaard and R. J. Francey, 1982: Winds in the atmospheric boundary

layer - prediction and observation. J. Atmos. Sci., 39, 1307-1316.

Grell, G. A., Y.-H. Kuo, and R. Pasch, 1991: Semi-prognostic tests of cumulus

parameterization schemes in the middle latitudes. Mon. Wea. Rev., 119, 5-31.

, 1993: Prognostic evaluation of assumptions used by cumulus parameterizations.Mon.

Wea. Rev., 121, 764-787.

Haagenson, P.L., J. Dudhia, G. A. Grell, and D. R. Stauer, 1994: The Penn State /

NCAR mesoscale model (MM5) source code documentation. NCAR Technical Note,

NCAR/TN-392+STR, 200 pp. [Available from NCAR Information Services, P.O. Box

3000, Boulder, CO 80307.]

Hsie, E.-Y. and R.A. Anthes 1984: Simulations of frontogenesis in a moist atmosphere using

alternative parameterizations of condensation and precipitation. J. Atmos. Sci., 41, 2701-2716.

Hoke, J.E. and R.A. Anthes, 1976: The initialization of numerical models by a dynamical

initialization technique. Mon. Wea. Rev., 104, 1551-1556.

Houze, R. A., C.-P. Cheng, C. A. Leary, and J. F. Gamache, 1979: Diagnosis of cloud

mass and heat 
uxes from radar and synoptic data. J. Atmos. Sci., 37, 754-773.

Kessler, E., 1969: On the distribution and continuity of water substance in atmospheric

circulation. Meteor. Monogr., 10, No.32, Amer. Meteor. Soc., 84 pp.

Kistler, R.E., 1974: A study of data assimilation techniques in an autobarotropic primitive

equation channel model. M.S. thesis, The Pennsylvania State University, 84 pp.

[Available from the Dept. of Meteorology, The Pennsylvania State University, University Park, PA 16802]

Klemp, J.B., and D.R. Durran, 1983: An upper boundary condition permitting internal

119
gravity wave radiation in numerical mesoscale models. Mon. Wea. Rev., 111, 430{

444.

Klemp, J.B., and R.B.Wilhelmson, 1978: Simulations of three-dimensional convective

storm dynamics. J. Atmos. Sci., 35, 1070{1096.

Kuo, H. L., 1965: On formation and intensication of tropical cyclones through latent heat

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

, 1974: Further studies of the parameterization of the eect of cumulus convection on

large scale 
ow. J. Atmos. Sci., 31, 1232-1240.

Lacis, A. A., and J. E. Hansen, 1974: A parameterization for the absorption of solar

radiation in the earth's atmosphere. J. Atmos. Sci., 31, 118-133.

Lord, S., 1978: Development and observational verication of a cumulus cloud

parameterization. Ph. D. dissertation, University of California, Los Angeles, 359

pp.

, 1982: Interaction of a cumulus cloud ensemble with the large scale environment. Part

III: Semi-prognostic test of the Arakawa-Schubert parameterization. J. Atmos. Sci., 39, 88-103.

Madala, R. V., 1981: Finite-Dierence Techniques for Vectorized Fluid Dynamics

Calculations. Edited by D. L. Book, Springer Verlag, New York.

Marshall, J. S., and W. M. Palmer, 1948: The distribution of raindrops with size. J.

Meteor., 5, 154-166.

Montieth, J. L., 1961: An empirical method for estimating long wave radiation exchanges

in the British Isles. Quart. J. Roy. Meteor. Soc., 87, 171-179.

Perkey, D. J., and C. W. Kreitzberg, 1976: A time-dependent lateral boundary scheme for

limited area primitive equation models. Mon. Wea. Rev., 104, 744-755.

Richtmeyer, R. D., 1957: Dierence methods for Initial-Value Problems. Interscience, New

York, 283 pp.

120
Rogers, C. D., 1967: The use of emissivity in atmosphere radiation calculation. Quart. J.

Roy. Metero. Soc., 93, 43-52.

Schubert, W. H., 1974: Cumulus parameterization theory in terms of feedback and control.

Atmos. Sci. Pap. No. 226, Colorado State University, 19 pp.

Sekhon, R.S., and R.C. Srivastava, 1970: Snow size spectra and re
ectivity. J. Atmos.

Sci., 27, 299{307.

Shapiro, R., 1970: Smoothing, ltering and boundary eects. Rev. of Geophys. Space

Phy.,8,359-387.

Skamarock, W.C., and J.B. Klemp, 1992: The stability of time-split numerical methods

for the hydrostatic and nonhydrostatic elastic equations. Mon. Wea. Rev. 120,

2109{2127.

Smagorinsky, J., S. Manabe, and J. L. Holloway, Jr., 1965: Numerical results from a

nine-level general circulation model of the atmosphere, Mon. Wea. Rev., 93, 727-768.

Smolarkiewicz, P. K., and G. A. Grell, 1992: A class of monotone interpolation schemes.

J. Comp. Phys., 101, 431-440.

Stauer, D.R. and N.L. Seaman, 1990: Use of four-dimensional data assimilation in a

limited-area mesoscale model. Part I: Experiments with synoptic-scale data. Mon.

Wea. Rev., 118, 1250-1277.

, N.L. Seaman and F.S. Binkowski, 1991: Use of four- dimensional data assimilation

in a limited-area mesoscale model. Part II: Eects of data assimilation within the

planetary boundary layer. Mon. Wea. Rev., 119, 734-754.

, and N.L. Seaman, 1994: On multi-scale four-dimensional data assimilation. J. Appl.

Meteor., 33, 416-434.

Stephens, G. L., 1978: Radiation proles in extended water clouds. Part II:

Parameterization schemes. J. Atmos. Sci., 35, 2123-2132.

, 1984: The parameterization of radiation for numerical weather prediction and climate

121
models. Mon. Wea. Rev., 112, 826-867.

, and P. J. Webster, 1979: Sensitivity of radiative forcing to variable cloud and moisture.

J. Atmos. Sci., 36, 1542-1556.

Zhang, D.-L., and R.A. Anthes, 1982: A high-resolution model of the planetary boundary

layer{sensitivity tests and comparisons with SESAME-79 data. J. Appl. Meteor., 21, 1594-1609.

, H.-R. Chang, N. L. Seaman, T. T. Warner, and J. M. Fritsch, 1986: A two-way

interactive nesting procedure with variable terrain resolution. Mon. Wea. Rev., 114, 1330-1339.

122