NCAR/TN–468+STR

NCAR TECHNICAL NOTE

June 2005

A Description of the

Advanced Research WRF Version 2

William C. Skamarock

Joseph B. Klemp

Jimy Dudhia

David O. Gill

Dale M. Barker

Wei Wang

Jordan G. Powers

Mesoscale and Microscale Meteorology Division

National Center for Atmospheric Research

Boulder, Colorado, USA
NCAR TECHNICAL NOTES

The Technical Note series provides an outlet for a variety of NCAR manuscripts that contribute

in specialized ways to the body of scientific knowledge but which are not suitable for journal,

monograph, or book publication. Reports in this series are issued by the NCAR Scientific Divisions; copies may be obtained on request from the Publications Office of NCAR. Designation

symbols for the series include:

EDD: Engineering, Design, or Development Reports

Equipment descriptions, test results, instrumentation,

and operating and maintenance manuals.

IA: Instructional Aids

Instruction manuals, bibliographies, film supplements,

and other research or instructional aids.

PPR: Program Progress Reports

Field program reports, interim and working reports,

survey reports, and plans for experiments.

PROC: Proceedings

Documentation of symposia, colloquia, conferences, workshops,

and lectures. (Distribution may be limited to attendees.)

STR: Scientific and Technical Reports

Data compilations, theoretical and numerical

investigations, and experimental results.

The National Center for Atmospheric Research (NCAR) is operated by the University Corporation for Atmospheric Research (UCAR) and is sponsored by the National Science Foundation.

Any opinions, findings, conclusions, or recommendations expressed in this publication are those

of the authors and do not necessarily reflect the views of the National Science Foundation.
Contents

Acknowledgments ix

1 Introduction 1

1.1 Advanced Research WRF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Major Features of the ARW System . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Governing Equations 5

2.1 Vertical Coordinate and Variables . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2 Flux-Form Euler Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.3 Inclusion of Moisture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.4 Map Projections, Coriolis and Curvature Terms . . . . . . . . . . . . . . . . . . 7

2.5 Perturbation Form of the Governing Equations . . . . . . . . . . . . . . . . . . . 9

3 Model Discretization 11

3.1 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3.1.1 Runge-Kutta Time Integration Scheme . . . . . . . . . . . . . . . . . . . 11

3.1.2 Acoustic Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3.1.3 Full Time-Split Integration Sequence . . . . . . . . . . . . . . . . . . . . 14

3.1.4 Diabatic Forcing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.1.5 Hydrostatic Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.2 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.2.1 Acoustic Step Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.2.2 Coriolis and Curvature Terms . . . . . . . . . . . . . . . . . . . . . . . . 18

3.2.3 Advection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.3 Stability Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.3.1 RK3 Time Step Constraint . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.3.2 Acoustic Time Step Constraint . . . . . . . . . . . . . . . . . . . . . . . 22

4 Turbulent Mixing and Model Filters 23

4.1 Explicit Spatial Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

4.1.1 Horizontal and Vertical Diffusion on Coordinate Surfaces . . . . . . . . . 23

4.1.2 Horizontal and Vertical Diffusion in Physical Space . . . . . . . . . . . . 24

4.1.3 Computation of the Eddy Viscosities . . . . . . . . . . . . . . . . . . . . 26

4.1.4 TKE equation for the 1.5 Order Turbulence Closure . . . . . . . . . . . . 27

4.2 Filters for the Time-split RK3 scheme . . . . . . . . . . . . . . . . . . . . . . . . 29

4.2.1 Divergence Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

i
4.2.2 External Mode Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.2.3 Semi-Implicit Acoustic Step Off-centering . . . . . . . . . . . . . . . . . . 29

4.3 Other Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4.3.1 Gravity Wave Absorbing Layer . . . . . . . . . . . . . . . . . . . . . . . 30

4.3.2 Rayleigh Damping Layer . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4.3.3 Vertical-Velocity Damping . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5 Initial Conditions 33

5.1 Initialization for Idealized Conditions . . . . . . . . . . . . . . . . . . . . . . . . 33

5.2 Initialization for Real-Data Conditions . . . . . . . . . . . . . . . . . . . . . . . 35

5.2.1 Use of the Standard Initialization by the ARW . . . . . . . . . . . . . . . 35

5.2.2 Reference and Perturbation State . . . . . . . . . . . . . . . . . . . . . . 36

5.2.3 Generating Lateral Boundary Data . . . . . . . . . . . . . . . . . . . . . 38

5.2.4 Masking of Surface Fields . . . . . . . . . . . . . . . . . . . . . . . . . . 38

6 Lateral Boundary Conditions 39

6.1 Periodic Lateral Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 39

6.2 Open Lateral Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 39

6.3 Symmetric Lateral Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . 40

6.4 Specified Lateral Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 40

7 Nesting 43

7.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

7.2 Nesting and Staggering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

7.3 Nested Lateral Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . 47

7.4 Steps to Generate a Nest Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

8 Physics 51

8.1 Microphysics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

8.1.1 Kessler scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

8.1.2 Purdue Lin scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

8.1.3 WRF Single-Moment 3-class (WSM3) scheme . . . . . . . . . . . . . . . 53

8.1.4 WSM5 scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

8.1.5 WSM6 scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

8.1.6 Eta Grid-scale Cloud and Precipitation (2001) scheme . . . . . . . . . . . 53

8.1.7 Thompson et al. scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

8.2 Cumulus parameterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

8.2.1 Kain-Fritsch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

8.2.2 Betts-Miller-Janjic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

8.2.3 Grell-Devenyi ensemble . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

8.3 Surface Layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

8.3.1 Similarity theory (MM5) . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

8.3.2 Similarity theory (Eta) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

8.4 Land-Surface Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

8.4.1 5-layer thermal diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

8.4.2 Noah LSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

ii
8.4.3 Rapid Update Cycle (RUC) Model LSM . . . . . . . . . . . . . . . . . . 58

8.5 Planetary Boundary Layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

8.5.1 Medium Range Forecast Model (MRF) PBL . . . . . . . . . . . . . . . . 59

8.5.2 Yonsei University (YSU) PBL . . . . . . . . . . . . . . . . . . . . . . . . 59

8.5.3 Mellor-Yamada-Janjic (MYJ) PBL . . . . . . . . . . . . . . . . . . . . . 59

8.6 Atmospheric Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

8.6.1 Rapid Radiative Transfer Model (RRTM) Longwave . . . . . . . . . . . . 61

8.6.2 Eta Geophysical Fluid Dynamics Laboratory (GFDL) Longwave . . . . . 61

8.6.3 Eta Geophysical Fluid Dynamics Laboratory (GFDL) Shortwave . . . . . 61

8.6.4 MM5 (Dudhia) Shortwave . . . . . . . . . . . . . . . . . . . . . . . . . . 61

8.6.5 Goddard Shortwave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

8.7 Physics Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

9 Variational Data Assimilation 65

9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

9.2 Improvements to the WRF-Var Algorithm . . . . . . . . . . . . . . . . . . . . . 67

9.2.1 Improved vertical interpolation . . . . . . . . . . . . . . . . . . . . . . . 67

9.2.2 Improved minimization and “outer loop” . . . . . . . . . . . . . . . . . . 67

9.2.3 Flexible choice of control variables . . . . . . . . . . . . . . . . . . . . . . 67

9.2.4 First Guess at Appropriate Time (FGAT) . . . . . . . . . . . . . . . . . 69

9.2.5 Radar Data Assimilation . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

9.2.6 Unified Regional/Global 3D-Var Assimilation . . . . . . . . . . . . . . . 69

9.3 Background Error Covariances . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

9.3.1 Removal of time-mean . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

9.3.2 Multivariate Covariances: Regression coefficients and unbalanced variables 72

9.3.3 Vertical Covariances: Eigenvectors/eigenvalues and control variable projections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

9.3.4 Horizontal Covariances: Recursive filter lengthscale (regional), or power

spectra (global) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

Appendices

A Physical Constants 75

B List of Symbols 77

C Acronyms 81

References 83

iii
iv
List of Figures

1.1 WRF system components. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.1 ARW η coordinate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3.1 Time step integration sequence. Here n represents the number of acoustic time

steps for a given substep of the RK3 integration, and ns is the ratio of the RK3

time step to the acoustic time step for the second and third RK3 substeps. . . . 14

3.2 Horizontal and vertical grids of the ARW . . . . . . . . . . . . . . . . . . . . . . 17

5.1 Schematic showing the data flow and program components in the SI, and how

the SI feeds initial data to the ARW. Letters in the rectangular boxes indicate

program names. GRID GEN: defines the model domain and creates static files

of terrestrial data. GRIB PREP: decodes GriB data. HINTERP: interpolates

meteorological data to the model domain. VINTERP: vertically interpolates

data to model coordinate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

6.1 Specified and relaxation zones for a grid with a single specified row and column,

and four rows and columns for the relaxation zone. These are typical values used

for a specified lateral boundary condition for a real-data case. . . . . . . . . . . 41

7.1 1-way and 2-way nesting options in the ARW. . . . . . . . . . . . . . . . . . . . 44

7.2 Various nest configurations for multiple grids. (a) Telescoping nests. (b) Nests at

the same level with respect to a parent grid. (c) Overlapping grids: not allowed

(d) Inner-most grid has more than one parent grid: not allowed . . . . . . . . . 45

7.3 Arakawa-C grid staggering for a portion of a parent domain and an imbedded

nest domain with a 3:1 grid size ratio. The solid lines denote coarse grid cell

boundaries, and the dashed lines are the boundaries for each fine grid cell. The

horizontal components of velocity (“U” and “V”) are defined along the normal cell

face, and the thermodynamic variables (“θ”) are defined at the center of the grid

cell (each square). The bold typeface variables along the interface between the

coarse and the fine grid define the locations where the specified lateral boundaries

for the nest are in effect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

7.4 Similar to Fig. 7.3, but with a 2:1 grid-distance ratio. . . . . . . . . . . . . . . 48

7.5 Nest grid integration sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

9.1 Sketch showing the relationship between datasets (circles), and algorithms (rectangles) of the ARW system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

v
9.2 Sketch of the role of Stage 0 converters in transforming model-specific data (e.g.,

ARW, KMA global model, etc.) to standard perturbation fields and relevant

metadata (e.g., latitude, height, land/sea, etc.). . . . . . . . . . . . . . . . . . . 72

vi
List of Tables

3.1 Maximum stable Courant numbers for one-dimensional linear advection. From

Wicker and Skamarock (2002). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

8.1 Microphysics Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

8.2 Cumulus Parameterization Options . . . . . . . . . . . . . . . . . . . . . . . . . 55

8.3 Land Surface Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

8.4 Planetary Boundary Layer Options . . . . . . . . . . . . . . . . . . . . . . . . . 59

8.5 Radiation Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

8.6 Physics Interactions. Columns correspond to model physical processes: radiation

(Rad), microphysics (MP), cumulus parameterization (CP), planetary boundary

layer/vertical diffusion (PBL), and surface physics (Sfc). Rows corresponds to

model variables where i and o indicate whether a variable is input or output

(updated) by a physical process. . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

9.1 The definitions of the various stages of the control variable transform given by

(9.2) for the unified global/regional WRF-Var system. Indices (i, j, k) refer to

grid-point space, index m to vertical mode, and l, n to global spectral mode.

The variables are: u, v: velocity components; T: temperature; q: mixing ratio;

ps: surface pressure; ψ: streamfunction; χ: velocity potential; r: relative humidity. The subscript u indicates an unbalanced field. The acronym RF stands for

recursive filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

vii
viii
Acknowledgments

Many people beyond this document’s author list have contributed to the development of

dynamics and physics. Regarding the dynamics solver development, we would like to thank

Louis Wicker for his assistance with the Runge-Kutta integration scheme. George Bryan and

Jason Knievel contributed considerably in the turbulent mixing algorithms and other model

filters.

Regarding the physics, we would like to thank Tom Black, John Brown, Kenneth Campana,

Fei Chen, Shu-Hua Chen, Ming-Da Chou, Mike Ek, Brad Ferrier, Georg Grell, Bill Hall, Songyou

Hong, Zavisa Janjic, Jack Kain, Jeong-Ock Lim, Ken Mitchell, Eli Mlawer, Tanya Smirnova,

Wei-Kuo Tao, Greg Thompson, and others.

The development of WRF-Var represents an international team effort. We would like to

acknowledge the following people for their contributions to the WRF-Var system: Yong-Run

Guo, Wei Huang, Qingnong Xiao, Syed Rizvi, Francois Vandenberghe, Mike McAtee, Roy Peck,

Wan-Shu Wu, Dezso Devenyi, Jianfeng Gu, Mi-Seon Lee, Ki-Han Youn, Eunha Lim, Hyun-Cheol

Shin, Shu-Hua Chen, and Hui-Chuan Lin.

We would like express our appreciation to John Michalakes for developing the WRF Software

Framework. We would also like thank Jacques Middlecoff, Tom Henderson, Al Bourgeois, Matt

Bettencourt, Todd Hutchinson, Kent Yang, and Dan Schaffer for their contributions to the

software development.

Finally, we would like to thank Hans Huang and Louisa Nance for their careful review of

this document. We also thank Cindy Bruy`ere for assistance with the graphics.

The Weather Research and Forecasting model effort is supported by the National Science

Foundation, the Department of Defense, the National Oceanic and Atmospheric Administration,

and the Federal Aviation Administration.

ix
x
Chapter 1

Introduction

The development of the Weather Research and Forecasting (WRF) modeling system is a multiagency effort intended to provide a next-generation mesoscale forecast model and data assimilation system that will advance both the understanding and prediction of mesoscale weather

and accelerate the transfer of research advances into operations. The model is being developed as a collaborative effort among the NCAR Mesoscale and Microscale Meteorology (MMM)

Division, the National Oceanic and Atmospheric Administration’s (NOAA) National Centers

for Environmental Prediction (NCEP) and Forecast System Laboratory (FSL), the Department

of Defense’s Air Force Weather Agency (AFWA) and Naval Research Laboratory (NRL), the

Center for Analysis and Prediction of Storms (CAPS) at the University of Oklahoma, and the

Federal Aviation Administration (FAA), along with the participation of a number of university

scientists.

The WRF model is designed to be a flexible, state-of-the-art, portable code that is efficient

in a massively parallel computing environment. A modular single-source code is maintained

that can be configured for both research and operations. It offers numerous physics options,

thus tapping into the experience of the broad modeling community. Advanced data assimilation

systems are being developed and tested in tandem with the model. WRF is maintained and

supported as a community model to facilitate wide use, particularly for research and teaching,

in the university community. It is suitable for use in a broad spectrum of applications across

scales ranging from meters to thousands of kilometers. Such applications include research and

operational numerical weather prediction (NWP), data assimilation and parameterized-physics

research, downscaling climate simulations, driving air quality models, atmosphere-ocean coupling, and idealized simulations (e.g boundary-layer eddies, convection, baroclinic waves). With

WRF as a common tool in the university and operational centers, closer ties will be promoted

between these communities, and research advances will have a direct path to operations. These

hallmarks make the WRF modeling system unique in the history of NWP in the United States.

The principal components of the WRF system are depicted in Figure 1.1. The WRF Software

Framework (WSF) provides the infrastructure that accommodates multiple dynamics solvers,

physics packages that plug into the solvers through a standard physics interface, programs

for initialization, and the WRF variational data assimilation (WRF-Var) system. As of this

writing there are two dynamics solvers in the WSF: the Advanced Research WRF (ARW) solver

(originally referred to as the Eulerian mass or “em” solver) developed primarily at NCAR,

and the NMM (Nonhydrostatic Mesoscale Model) solver developed at NCEP, which will be

documented and supported to the community by the Developmental Testbed Center (DTC).

1
Obs Data,

Analyses,

Forecast

Post Processors,

Verification

WRF Software Framework

Standard Physics Interface

Physics Packages

Dynamics Solvers

ARW Solver

NMM Solver

...

Initialization

WRF-Var

Figure 1.1: WRF system components.

While there are multiple solvers, and while not all physics are available to both solvers, the

WSF is common to all components.

1.1 Advanced Research WRF

The ARW system consists of the ARW dynamics solver together with other components of the

WRF system needed to produce a simulation. Thus, it also encompasses physics schemes, initialization routines, and a data assimilation package. The ARW shares the WSF, the framework

common to all WRF modeling system components, including the NMM solver. Similarly, the

physics packages are available to both the ARW and NMM solvers. In this vein, it should be

understood that the association of a component of the WRF system with the ARW does not

preclude it from being a component of any other WRF configuration. The following section

highlights the major features of the ARW system, Version 2, and reflects elements of WRF

Version 2, which was first released in May 2004.

This technical note will focus on the scientific and algorithmic approaches in the ARW.

Discussed are the ARW solver, available physics options, initialization capabilities, boundary

conditions, and grid-nesting techniques. The WSF provides the software infrastructure for all

WRF configurations and is documented separately (Michalakes et al., 1999, 2004). The WRFVar program, a component of the broader WRF system, has been adapted from MM5 3DVAR

(Barker et al., 2004) and is encompassed within the ARW. As a separate document detailing the

broader WRF-Var system will be forthcoming, this technical note will focus on a summary of

the changes and updates implemented to adapt this data assimilation capability from the MM5

to WRF. For those seeking information on running the ARW modeling system, details on its

use are contained in the ARW User’s Guide (Wang et al., 2004).

2
1.2 Major Features of the ARW System

ARW Solver

• Equations: Fully compressible, Euler nonhydrostatic with a run-time hydrostatic option

available. Conservative for scalar variables.

• Prognostic Variables: Velocity components u and v in Cartesian coordinate, vertical velocity

w, perturbation potential temperature, perturbation geopotential, and perturbation surface pressure of dry air. Optionally, turbulent kinetic energy and any number of scalars

such as water vapor mixing ratio, rain/snow mixing ratio, and cloud water/ice mixing

ratio.

• Vertical Coordinate: Terrain-following hydrostatic-pressure, with vertical grid stretching permitted. Top of the model is a constant pressure surface.

• Horizontal Grid: Arakawa C-grid staggering.

• Time Integration: Time-split integration using a 3rd order Runge-Kutta scheme with smaller

time step for acoustic and gravity-wave modes.

• Spatial Discretization: 2nd to 6th order advection options in horizontal and vertical.

• Turbulent Mixing and Model Filters: Sub-grid scale turbulence formulation in both coordinate and physical space. Divergence damping, external-mode filtering, vertically implicit

acoustic step off-centering. Explicit filter option also available.

• Initial Conditions: Three dimensional for real-data, and one-, two- and three-dimensional

using idealized data. A number of test cases are provided.

• Lateral Boundary Conditions: Periodic, open, symmetric, and specified options available.

• Top Boundary Conditions: Gravity wave absorbing (diffusion or Rayleigh damping). w = 0

top boundary condition at constant pressure level.

• Bottom Boundary Conditions: Physical or free-slip.

• Earth’s Rotation: Full Coriolis terms included.

• Mapping to Sphere: Three map projections are supported for real-data simulation: polar

stereographic, Lambert-conformal, and Mercator. Curvature terms included.

• Nesting: One-way, two-way, and moving nests.

Model Physics

• Microphysics: Bulk schemes ranging from simplified physics suitable for mesoscale modeling

to sophisticated mixed-phase physics suitable for cloud-resolving modeling.

• Cumulus parameterizations: Adjustment and mass-flux schemes for mesoscale modeling including NWP.

• Surface physics: Multi-layer land surface models ranging from a simple thermal model to full

vegetation and soil moisture models, including snow cover and sea ice.

• Planetary boundary layer physics: Turbulent kinetic energy prediction or non-local K schemes.

• Atmospheric radiation physics: Longwave and shortwave schemes with multiple spectral

bands and a simple shortwave scheme. Cloud effects and surface fluxes are included.

3
WRF-Var System

• Incremental formulation of the model-space cost function.

• Quasi-Newton or conjugate gradient minimization algorithms.

• Analysis increments on un-staggered Arakawa-A grid.

• Representation of the horizontal component of background error B via recursive filters (regional) or power spectra (global). The vertical component is applied through projection

onto climatologically-averaged eigenvectors of vertical error. Horizontal/vertical errors are

non-separable (horizontal scales vary with vertical eigenvector).

• Background cost function (Jb) preconditioning via a control variable transform U defined as

B = UUT.

• Flexible choice of background error model and control variables.

• Climatological background error covariances estimated via either the NMC-method of averaged forecast differences or suitably averaged ensemble perturbations.

• Unified 3D-Var (4D-Var under development), global and regional, multi-model capability.

WRF Software Framework

• Highly modular, single-source code for maintainability.

• Portable across a range of available computing platforms.

• Support for multiple dynamics solvers and physics modules.

• Separation of scientific codes from parallelization and other architecture-specific codes.

• Input/Output Application Program Interface (API) enabling various external packages to be

installed with WRF, hence allowing WRF to easily support various data formats.

• Efficient execution on a range of computing platforms (distributed and shared memory, vector

and scalar types).

• Use of Earth System Modeling Framework (ESMF) timing package.

• Model coupling API enabling WRF to be coupled with other models such as ocean, and land

models.

4
Chapter 2

Governing Equations

The ARW dynamics solver integrates the compressible, nonhydrostatic Euler equations. The

equations are cast in flux form using variables that have conservation properties, following the

philosophy of Ooyama (1990). The equations are formulated using a terrain-following mass

vertical coordinate (Laprise, 1992). In this chapter we define the vertical coordinate and present

the flux form equations in Cartesian space, we extend the equations to include the effects of

moisture in the atmosphere, and we further augment the equations to include projections to the

sphere.

2.1 Vertical Coordinate and Variables η

1.0

0.8

0.6

0.4

0.2

0

Pht = constant

Phs

Figure 2.1: ARW η coordinate.

The ARW equations are formulated using a

terrain-following hydrostatic-pressure vertical coordinate denoted by η and defined as

η = (ph −pht)/μ where μ = phs −pht. (2.1)

ph is the hydrostatic component of the pressure,

and phs and pht refer to values along the surface

and top boundaries, respectively. The coordinate

definition (2.1), proposed by Laprise (1992), is

the traditional σ coordinate used in many hydrostatic atmospheric models. η varies from a value

of 1 at the surface to 0 at the upper boundary of

the model domain (Fig. 2.1). This vertical coordinate is also called a mass vertical coordinate.

Since μ(x, y) represents the mass per unit area

within the column in the model domain at (x, y),

the appropriate flux form variables are

V = μv = (U, V, W), Ω = μη, ̇ Θ = μθ. (2.2)

v = (u, v, w) are the covariant velocities in the

two horizontal and vertical directions, respectively, while ω =  ̇η is the contravariant ‘vertical’

5
velocity. θ is the potential temperature. Also appearing in the governing equations of the ARW

are the non-conserved variables φ = gz (the geopotential), p (pressure), and α = 1/ρ (the inverse

density).

2.2 Flux-Form Euler Equations

Using the variables defined above, the flux-form Euler equations can be written as

∂tU + (∇ · Vu) − ∂x(pφη) + ∂η(pφx) = FU (2.3)

∂tV + (∇ · Vv) − ∂y(pφη) + ∂η(pφy) = FV (2.4)

∂tW + (∇ · Vw) − g(∂ηp − μ) = FW (2.5)

∂tΘ + (∇ · Vθ) = FΘ (2.6)

∂tμ + (∇ · V) = 0 (2.7)

∂tφ + μ

−1

[(V · ∇φ) − gW] = 0 (2.8)

along with the diagnostic relation for the inverse density

∂ηφ = −αμ, (2.9)

and the equation of state

p = p0(Rdθ/p0α)

γ

. (2.10)

In (2.3) – (2.10), the subscripts x, y and η denote differentiation,

∇ · Va = ∂x(U a) + ∂y(V a) + ∂η(Ωa),

and

V · ∇a = U ∂xa + V ∂ya + Ω∂ηa,

where a represents a generic variable. γ = cp/cv = 1.4 is the ratio of the heat capacities for dry

air, Rd is the gas constant for dry air, and p0 is a reference pressure (typically 105 Pascals). The

right-hand-side (RHS) terms FU , FV , FW , and FΘ represent forcing terms arising from model

physics, turbulent mixing, spherical projections, and the earth’s rotation.

The prognostic equations (2.3) – (2.8) are cast in conservative form except for (2.8) which

is the material derivative of the definition of the geopotential. (2.8) could be cast in flux form

but we find no advantage in doing so since μφ is not a conserved quantity. We could also

use a prognostic pressure equation in place of (2.8) (see Laprise, 1992), but pressure is not a

conserved variable and we could not use a pressure equation and the conservation equation for Θ

(2.6) because they are linearly dependent. Additionally, prognostic pressure equations have the

disadvantage of possessing a mass divergence term multiplied by a large coefficient (proportional

to the sound speed) which makes spatial and temporal discretization problematic. It should be

noted that the relation for the hydrostatic balance (2.9) does not represent a constraint on the

solution, rather it is a diagnostic relation that formally is part of the coordinate definition. In the

hydrostatic counterpart to the nonhydrostatic equations, (2.9) replaces the vertical momentum

equation (2.5) and it becomes a constraint on the solution.

6
2.3 Inclusion of Moisture

In formulating the moist Euler equations, we retain the coupling of dry air mass to the prognostic

variables, and we retain the conservation equation for dry air (2.7), as opposed to coupling the

variables to the full (moist) air mass and hence introducing source terms in the mass conservation

equation (2.7). Additionally, we define the coordinate with respect to the dry-air mass. Based

on these principles, the vertical coordinate can be written as

η = (pdh − pdht)/μd (2.11)

where μd represents the mass of the dry air in the column and pdh and pdht represent the

hydrostatic pressure of the dry atmosphere and the hydrostatic pressure at the top of the dry

atmosphere. The coupled variables are defined as

V = μdv, Ω = μdη, ̇ Θ = μdθ. (2.12)

With these definitions, the moist Euler equations can be written as

∂tU + (∇ · Vu)η + μdα∂xp + (α/αd)∂ηp∂xφ = FU (2.13)

∂tV + (∇ · Vv)η + μdα∂yp + (α/αd)∂ηp∂yφ = FV (2.14)

∂tW + (∇ · Vw)η − g[(α/αd)∂ηp − μd] = FW (2.15)

∂tΘ + (∇ · Vθ)η = FΘ (2.16)

∂tμd + (∇ · V)η = 0 (2.17)

∂tφ + μ

−1

d

[(V · ∇φ)η − gW] = 0 (2.18)

∂tQm + (V · ∇qm)η = FQm (2.19)

with the diagnostic equation for dry inverse density

∂ηφ = −αdμd (2.20)

and the diagnostic relation for the full pressure (vapor plus dry air)

p = p0(Rdθm/p0αd)

γ

(2.21)

In these equations, αd is the inverse density of the dry air (1/ρd) and α is the inverse density

taking into account the full parcel density α = αd(1 + qv + qc + qr + qi + ...)

−1 where q∗ are

the mixing ratios (mass per mass of dry air) for water vapor, cloud, rain, ice, etc. Additionally,

θm = θ(1 + (Rv/Rd)qv) ≈ θ(1 + 1.61qv), and Qm = μdqm; qm = qv, qc, qi, ... .

2.4 Map Projections, Coriolis and Curvature Terms

The ARW solver currently supports three projections to the sphere— the Lambert conformal,

polar stereographic, and Mercator projections. These projections are described in Haltiner and

Williams (1980). These projections, and the ARW implementation of the map factors, assume

that the map factor transformations for x are y are identical at a given point; that is, the

7
transformation is isotropic. Anisotropic transformations, such as a latitude-longitude grid, can

be accommodated by defining separate map factors for the x and y transformations.

In the ARW’s computational space, ∆x and ∆y are constants. Orthogonal projections to

the sphere require that the physical distances between grid points in the projection vary with

position in the grid. To transform the governing equations, a map scale factor m is defined as

the ratio of the distance in computational space to the corresponding distance on the earth’s

surface:

m =

(∆x, ∆y)

distance on the earth. (2.22)

The ARW solver includes the map-scale factors in the governing equations by redefining the

momentum variables as

U = μdu/m, V = μdv/m, W = μdw/m, Ω = μdη/m.  ̇

Using these redefined momentum variables, the governing equations, including map factors and

rotational terms, can be written as

∂tU + m[∂x(Uu) + ∂y(V u)] + ∂η(Ωu) + μdα∂xp + (α/αd)∂ηp∂xφ = FU (2.23)

∂tV + m[∂x(Uv) + ∂y(V v)] + ∂η(Ωv) + μdα∂yp + (α/αd)∂ηp∂yφ = FV (2.24)

∂tW + m[∂x(Uw) + ∂y(V w)] + ∂η(Ωw) − m−1g[(α/αd)∂ηp − μd] = FW (2.25)

∂tΘ + m2[∂x(Uθ) + ∂y(V θ)] + m∂η(Ωθ) = FΘ (2.26)

∂tμd + m2[Ux + Vy] + m∂η(Ω) = 0 (2.27)

∂tφ + μ

−1

d

[m2(Uφx + V φy) + mΩφη − gW] = 0 (2.28)

∂tQm + m2[∂x(Uqm) + ∂y(V qm)] + m∂η(Ωqm) = FQm, (2.29)

and, for completeness, the diagnostic relation for the dry inverse density

∂ηφ = −αdμd, (2.30)

and the diagnostic equation for full pressure (vapor plus dry air)

p = p0(Rdθm/p0αd)

γ

. (2.31)

The right-hand-side terms of the momentum equations (2.23) – (2.25) contain the Coriolis

and curvature terms along with mixing terms and physical forcings. The Coriolis and curvature

terms can be written as follows:

FUcor = +f + u

∂m

∂y − v

∂m

∂x 

V − eW cos αr −

uW

re

(2.32)

FVcor = −



f + u

∂m

∂y − v

∂m

∂x 

U + eW sin αr −

vW

re

(2.33)

FWcor = +e(U cos αr − V sin αr) + 

uU + vV

re



, (2.34)

where αr is the local rotation angle between the y-axis and the meridians, ψ is the latitude,

f = 2Ωe sin ψ, e = 2Ωe cos ψ, Ωe is the angular rotation rate of the earth, and re is the radius of

8
the earth. In this formulation we have approximated the radial distance from the center of the

earth as the mean earth radius re, and we have not taken into account the change in horizontal

grid distance as a function of the radius. The terms containing m are the horizontal curvature

terms, those containing re relate to vertical (earth-surface) curvature, and those with e and f

are the Coriolis force.

In idealized cases, the map scale factor m = 1, f is often taken to be constant, and e = 0.

2.5 Perturbation Form of the Governing Equations

Before constructing the discrete solver, it is advantageous to recast the governing equations

using perturbation variables so as to reduce truncation errors in the horizontal pressure gradient

calculations in the discrete solver, in addition to reducing machine rounding errors in the vertical

pressure gradient and buoyancy calculations. For this purpose, new variables are defined as

perturbations from a hydrostatically-balanced reference state, and we define reference state

variables (denoted by overbars) that are a function of height only and that satisfy the governing

equations for an atmosphere at rest. That is, the reference state is in hydrostatic balance and

is strictly only a function of z. In this manner, p =  ̄p(z) + p

0

, φ = φ ̄(z) + φ

0

, α =  ̄α(z) + α

0

,

and μd =  ̄μd(x, y) + μ

0

d

. Because the η coordinate surfaces are generally not horizontal, the

reference profiles  ̄p, φ ̄, and  ̄α are functions of (x, y, η). The hydrostatically balanced portion

of the pressure gradients in the reference sounding can be removed without approximation to

the equations using these perturbation variables. The momentum equations (2.23) – (2.25) are

written as

∂tU + m[∂x(Uu) + ∂y(V u)] + ∂η(Ωu) + (μdα∂xp

0 + μdα0

∂xp ̄)

+(α/αd)(μd∂xφ

0 + ∂ηp0

∂xφ − μ

0

d∂xφ) = FU (2.35)

∂tV + m[∂x(Uv) + ∂y(V v)] + ∂η(Ωv) + (μdα∂yp

0 + μdα0

∂yp ̄)

+(α/αd)(μd∂yφ

0 + ∂ηp0

∂yφ − μ

0

d∂yφ) = FV (2.36)

∂tW + m[∂x(Uw) + ∂y(V w)] + ∂η(Ωw)

−m−1g(α/αd)[∂ηp

0 − μ ̄d(qv + qc + qr)] + m−1μ0

d

g = FW , (2.37)

and the mass conservation equation (2.27) and geopotential equation (2.28) become

∂tμ

0

d + m2

[∂xU + ∂yV ] + m∂ηΩ = 0 (2.38)

∂tφ

0 + μ

−1

d

[m2(Uφx + V φy) + mΩφη − gW] = 0. (2.39)

Remaining unchanged are the conservation equations for the potential temperature and scalars

∂tΘ + m2[∂x(Uθ) + ∂y(V θ)] + m∂η(Ωθ) = FΘ (2.40)

∂tQm + m2[∂x(Uqm) + ∂y(V qm)] + m∂η(Ωqm) = FQm. (2.41)

In the perturbation system the hydrostatic relation (2.30) becomes

∂ηφ

0 = −μ ̄dα0

d − αdμ

0

d

. (2.42)

Equations (2.35) – (2.41), together with the equation of state (2.21), represent the equations

solved in the ARW. The RHS terms in these equations include the Coriolis terms (2.32) –

9
(2.34), mixing terms (described in Chapter 4), and parameterized physics (described in Chapter

8). Also note that the equation of state (2.21) cannot be written in perturbation form because

of the exponent in the expression. For small perturbation simulations, accuracy for perturbation

variables can be maintained by linearizing (2.21) for the perturbation variables.

10
Chapter 3

Model Discretization

3.1 Temporal Discretization

The ARW solver uses a time-split integration scheme. Generally speaking, slow or low-frequency

(meteorologically significant) modes are integrated using a third-order Runge-Kutta (RK3) time

integration scheme, while the high-frequency acoustic modes are integrated over smaller time

steps to maintain numerical stability. The horizontally propagating acoustic modes (including

the external mode present in the mass-coordinate equations using a constant-pressure upper

boundary condition) are integrated using a forward-backward time integration scheme, and vertically propagating acoustic modes and buoyancy oscillations are integrated using a vertically

implicit scheme (using the acoustic time step). The time-split integration is similar to that first

developed by Klemp and Wilhelmson (1978) and analyzed by Skamarock and Klemp (1992).

The time-split RK3 scheme is described in general terms in Wicker and Skamarock (2002).

The primary differences between the descriptions found in the references and the ARW implementation are associated with our use of the mass vertical coordinate and a flux-form set of

equations, along with our use of perturbation variables for the acoustic component of the timesplit integration. The acoustic-mode integration is cast in the form of a correction to the RK3

integration.

3.1.1 Runge-Kutta Time Integration Scheme

The RK3 scheme, described in Wicker and Skamarock (2002), integrates a set of ordinary

differential equations using a predictor-corrector formulation. Defining the prognostic variables

in the ARW solver as Φ = (U, V, W, Θ, φ0, μ0, Qm) and the model equations as Φt = R(Φ), the

RK3 integration takes the form of 3 steps to advance a solution Φ(t) to Φ(t + ∆t):

Φ

∗ = Φt +

∆t

3

R(Φt) (3.1)

Φ

∗∗ = Φt +

∆t

2

R(Φ∗) (3.2)

Φ

t+∆t = Φt + ∆tR(Φ∗∗) (3.3)

where ∆t is the time step for the low-frequency modes (the model time step). In (3.1) – (3.3),

superscripts denote time levels. This scheme is not a true Runge-Kutta scheme per se because,

11
while it is third-order accurate for linear equations, it is only second-order accurate for nonlinear

equations. With respect to the ARW equations, the time derivatives Φt are the partial time

derivatives (the leftmost terms) in equations (2.35) – (2.41), and R(Φ) are the remaining terms

in (2.35) – (2.41).

3.1.2 Acoustic Integration

The high-frequency but meteorologically insignificant acoustic modes would severely limit the

RK3 time step ∆t in (3.1) – (3.3). To circumvent this time step limitation we use the approach

described in Wicker and Skamarock (2002). Additionally, to increase the accuracy of the splitting, we integrate a perturbation form of the governing equations using smaller acoustic time

steps within the RK3 large-time-step sequence. To form the perturbation equations for the RK3

time-split acoustic integration, we define small time step variables that are deviations from the

most recent RK3 predictor (denoted by the superscript t

∗ and representing either Φt

, Φ∗, or Φ∗∗

in (3.1) – (3.3)):

V00 = V − Vt

∗

, Ω

00 = Ω − Ωt

∗

, Θ

00 = Θ − Θt

∗

,

φ

00 = φ0 − φ0t

∗

, α00

d = α

0

d − α

0

d

t

∗

, μ00

d = μ

0

d − μ

0t

∗

d

.

The hydrostatic relation (i.e., the vertical coordinate definition) becomes

α

00

d = −

1

μ

t

∗

d



∂ηφ

00 + αt

∗

d μ

00

d



. (3.4)

Additionally, we also introduce a version of the equation of state that is linearized about t

∗

,

p

00 =

c

2

s

α

t

∗

d



Θ00

Θt

∗ −

α

00

d

α

t

∗

d

−

μ

00

d

μ

t

∗

d



, (3.5)

where c

2

s = γpt

∗

α

t

∗

d

is the square of the sound speed. The linearized state equation (3.5) and

the vertical coordinate definition (3.4) are used to cast the vertical pressure gradient in (2.37)

in terms of the model’s prognostic variables. By combining (3.5) and (3.4), the vertical pressure

gradient can be expressed as

∂ηp

00 = ∂η(C∂ηφ

00) + ∂η



c

2

s

α

t

∗

d

Θ00

Θt

∗



, (3.6)

where C = c

2

s

/μt

∗

α

t

∗ 2

. This linearization about the most recent large time step should be highly

accurate over the time interval of the several small time steps.

These variables along with (3.6) are substituted into the prognostic equations (2.35) – (2.41)

12
and lead to the acoustic time-step equations:

δτU

00 + μt

∗

α

t

∗

∂xp

00τ + (μt

∗

∂xp ̄)α

00τ + (α/αd)[μt

∗

∂xφ

00τ + (∂xφt

∗

)(∂ηp

00 − μ00)τ

] = RU

t

∗

(3.7)

δτV

00 + μt

∗

α

t

∗

∂yp

00τ + (μt

∗

∂yp ̄)α

00τ + (α/αd)[μt

∗

∂yφ

00τ + (∂yφ

t

∗

)(∂ηp

00 − μ00)τ

] = RV

t

∗

(3.8)

δτμ

00

d + m2

[∂xU

00 + ∂yV00]τ+∆τ + m∂ηΩ00τ+∆τ = Rμ

t

∗

(3.9)

δτΘ

00 + m2

[∂x(U

00θt

∗

) + ∂y(V

00θt

∗

)]τ+∆τ + m∂η(Ω00τ+∆τθ

t

∗

) = RΘ

t

∗

(3.10)

δτW00 − m−1g



(α/αd)

t

∗

∂η(C∂ηφ00) + ∂η



c

2

s

αt

∗

Θ00

Θt

∗



− μ

00

d

τ

= RW

t

∗

(3.11)

δτφ

00 +

1

μ

t

∗

d

[mΩ

00τ+∆τφη − gW00τ

] = Rφ

t

∗

. (3.12)

The RHS terms in (3.7) – (3.12) are fixed for the acoustic steps that comprise the time integration

of each RK3 sub-step (i.e., (3.1) – (3.3)), and are given by

R

t

∗

U = − m[∂x(Uu) + ∂y(V u)] + ∂η(Ωu) − (μdα∂xp

0 − μdα0

∂xp ̄)

− (α/αd)(μd∂xφ

0 − ∂ηp

0

∂xφ + μ

0

d∂xφ) + FU (3.13)

R

t

∗

V = − m[∂x(Uv) + ∂y(V v)] + ∂η(Ωv) − (μdα∂yp

0 − μdα0

∂yp ̄)

− (α/αd)(μd∂yφ

0 − ∂ηp0

∂yφ + μ

0

d∂yφ) + FV (3.14)

R

t

∗

μd = − m2

[∂xU + ∂yV ] + m∂ηΩ (3.15)

R

t

∗

Θ = − m2

[∂x(Uθ) + ∂y(V θ)] − m∂η(Ωθ) + FΘ (3.16)

R

t

∗

W = − m[∂x(Uw) + ∂y(V w)] + ∂η(Ωw)

+ m−1g(α/αd)[∂ηp

0 +  ̄μd(qv + qc + qr)] − m−1μ

0

d

g + FW (3.17)

R

t

∗

φ = − μ

−1

d

[m2(Uφx + V φy) + mΩφη − gW], (3.18)

where all variables in (3.13) – (3.18) are evaluated at time t

∗

(i.e., using Φt, Φ∗, or Φ∗∗ for the

appropriate RK3 sub-step in (3.1) – (3.3)). Equations (3.7) – (3.12) utilize the discrete acoustic

time-step operator

δτa =

a

τ+∆τ − aτ

∆τ

,

where ∆τ is the acoustic time step, and an acoustic time-step averaging operator

a

τ =

1 + β

2

a

τ+∆τ +

1 − β

2

a

τ

, (3.19)

where β is a user-specified parameter (see Section 4.2.3).

The integration over the acoustic time steps proceeds as follows. Beginning with the small

time-step variables at time τ , (3.7) and (3.8) are stepped forward to obtain U

00τ+∆τ

and V

00τ+∆τ

.

Both μ

00τ+∆τ

and Ω00τ+∆τare then calculated from (3.9). This is accomplished by first integrating

(3.9) vertically from the surface to the material surface at the top of the domain, which removes

the ∂ηΩ

00 term such that

δτμd = m2

Z 0

1

[∂xU

00 + ∂yV00]τ+∆τ

dη. (3.20)

13
After computing μ

00

d

τ+∆τ

from (3.20), Ω00τ+∆τis recovered by using (3.9) to integrate the ∂ηΩ

00

term vertically using the lower boundary condition Ω00 = 0 at the surface. Equation (3.10)

is then stepped forward to calculate Θ00τ+∆τ. Equations (3.11) and (3.12) are combined to

form a vertically implicit equation that is solved for W00τ+∆τsubject to the boundary condition

W00 = V00 · ∇h at the surface (z = h(x, y)) and p

0 = 0 along the model top. φ00τ+∆τ

is then

obtained from (3.12), and p

00τ+∆τ

and α

00

d

τ+∆τ

are recovered from (3.5) and (3.4).

3.1.3 Full Time-Split Integration Sequence

The time-split RK3 integration technique is summarized below. It consists of two primary

loops— an outer loop for the large-time-step Runge-Kutta integration, and an inner loop for

the acoustic mode integration.

Begin Time Step

Begin RK3 Loop: Steps 1, 2, and 3

(1) If RK3 step 1, compute and store FΦ

(i.e., physics tendencies for RK3 step, including mixing).

(2) Compute Rt

∗

Φ , (3.13)–(3.18)

Begin Acoustic Step Loop: Steps 1 → n,

RK3 step 1, n = 1, ∆τ = ∆t/3;

RK3 step 2, n = ns/2, ∆τ = ∆t/ns;

RK3 step 3, n = ns, ∆τ = ∆t/ns.

(3) Advance horizontal momentum, (3.7) and (3.8)

(4) Advance μd (3.9) and compute Ω00τ+∆τ

then advance Θ (3.10)

(5) Advance W and φ (3.11) and (3.12)

(6) Diagnose p

00 and α00 using (3.5) and (3.4)

End Acoustic Step Loop

(7) Scalar transport: Advance scalars (2.41)

over RK3 substep (3.1), (3.2) or (3.3)

(using mass fluxes U, V and Ω time-averaged over the acoustic steps).

(8) Compute p

0 and α0 using updated prognostic variables in (2.31) and (2.42)

End RK3 Loop

(9) Compute non-RK3 physics (currently microphysics), advance variables.

End Time Step

Figure 3.1: Time step integration sequence. Here n represents the number of acoustic time

steps for a given substep of the RK3 integration, and ns is the ratio of the RK3 time step to

the acoustic time step for the second and third RK3 substeps.

In the RK3 scheme, physics can be integrated within the RK3 time integration (using a time

14
forward step, i.e., step (1) in Fig. 3.1, or the RK3 time integration if higher temporal accuracy

is desired, i.e., in step (2)— implying a physics evaluation every RK3 substep) or external to it

using additive timesplitting, i.e., step (9).

Within the acoustic integration, the acoustic time step ∆τ is specified by the user through

the choice of ns (see Section 3.3.2). Within the first RK3 substep, however, a single acoustic

time step is used to advance the solution regardless of ns. Within the full RK3-acoustic timesplit

integration, this modified acoustic time step does not impose any additional stability constraints

(see Wicker and Skamarock, 2002).

The major costs in the model arise from the evaluation of the right hand side terms Rt

∗

in

(3.7) – (3.12). The efficiency of the RK3 timesplit scheme arises from the fact that the RK3

time step ∆t is much larger than the acoustic time step ∆τ , hence the most costly evaluations

are only performed in the less-frequent RK3 steps.

3.1.4 Diabatic Forcing

Within the RK3 integration sequence outlined in Fig. 3.1, the RHS term Rt

∗

Θ in the thermodynamic equation (3.10) contains contributions from the diabatic physics tendencies that are

computed in step (1) at the beginning of the first RK3 step. This diabatic forcing is integrated

within the acoustic steps (specifically, in step 4 in the time integration sequence shown in Fig.

3.1). Additional diabatic contributions are integrated in an additive-time-split manner in step

(9) after the RK3 update is complete. Thus, the diabatic forcing computed in step (9) (the

microphysics in the current release of the ARW) does not appear in Rt

∗

Θ from (3.10) in the

acoustic integration. We have discovered that this time splitting can excite acoustic waves and

can give rise to noise in the solutions for some applications. Note that the non-RK3 physics are

integrated in step (9) because balances produced in the physics are required at the end of the

time step (e.g., the saturation adjustment in the microphysics). So while moving these non-RK3

physics into step (1) would eliminate the noise, the balances produced by these physics would

be altered.

We have found that the excitation of the acoustic modes can be circumvented while leaving

the non-RK3 physics in step (9) by using the following procedure that is implemented in the

ARW. In step (1) of the integration procedure (Fig. 3.1), an estimate of the diabatic forcing

of Θ arising from the non-RK3 physics in step (9) is included in the diabatic forcing term

Rt

∗

Θ in (3.10) (which is advanced in step 4). This estimated diabatic forcing is then removed

from the updated Θ after the RK3 integration is complete and before the evaluation of the

non-RK3 physics in step (9). We use the diabatic forcing from the previous time step as the

estimated forcing; hence this procedure results in few additional computations outside of saving

the diabatic forcing between time steps.

3.1.5 Hydrostatic Option

A hydrostatic option is available in the ARW solver. The time-split RK3 integration technique

summarized in Fig. 3.1 is retained, including the acoustic step loop. Steps (5) and (6) in the

acoustic-step loop, where W and φ are advanced and p

00 and α00 are diagnosed, are replaced by

15
(1), the diagnosis of the hydrostatic pressure using the definition of the vertical coordinate

δηph =

αd

α

μd =



1 +Xqm



μd,

followed by (2), the diagnosis of αd using the equation of state (2.31) and the prognosed θ, and

(3), the diagnosis of the geopotential using the hydrostatic equation

δηφ

0 = −μ ̄dα0

d − μ

0

dαd.

The vertical velocity w can be diagnosed from the geopotential equation, but it is not needed in

the solution procedure. The acoustic step loop advances gravity waves, including the external

mode, when the hydrostatic option is used; there are no horizontally propagating acoustic modes

in this hydrostatic system.

3.2 Spatial Discretization

The spatial discretization in the ARW solver uses a C grid staggering for the variables as shown

in Fig. 3.2. That is, normal velocities are staggered one-half grid length from the thermodynamic

variables. The variable indices, (i, j) for the horizontal plane and (i, k) for the vertical plane,

indicate variable locations where (x, y, η) = (i∆x, j∆y, k∆η). We will denote the points where

θ is located as being mass points, and likewise we will denote locations where u, v, and w are

defined as u points, v points, and w points, respectively. Not shown in Fig. 3.2 are the column

mass μ, defined at the (i, j) points (mass points) on the discrete grid, the geopotential φ that

is defined at the w points, and the moisture variables qm are defined at the mass points. The

diagnostic variables used in the model, the pressure p and inverse density α, are computed at

mass points. The grid lengths ∆x and ∆y are constants in the model formulation; changes in

the physical grid lengths associated with the various projections to the sphere are accounted

for using the map factors introduced in Section 2.4. The vertical grid length ∆η is not a fixed

constant; it is specified in the initialization. The user is free to specify the η values of the model

levels subject to the constraint that η = 1 at the surface, η = 0 at the model top, and η decreases

monotonically between the surface and model top. Using these grid and variable definitions, we

can define the spatial discretization for the ARW solver.

3.2.1 Acoustic Step Equations

We begin by defining the column-mass-coupled variables relative to the uncoupled variables.

The vertical velocity is staggered only in k, so it can be coupled directly to the column mass

with no averaging or interpolation. The horizontal velocities are horizontally staggered relative

to the column mass such that the continuous variables are represented discretely as

U =

μdu

m

→

μd

xu

mx

, V =

μdv

m

→

μd

y

v

my

,

where the discrete operator a

x denotes a linear interpolation operator. The grid lengths ∆x and

∆y are constant, hence in this case the operator reduces to a

x = (ai+1/2 + ai−1/2)/2.

16
x

y

θi,j+1 θi+1,j+1

u θi,j θi+1,j i-1/2,j ui+1/2,j

ui-1/2,j+1 ui+1/2,j+1 ui+3/2,j+1

ui+3/2,j

vi,j-1/2 vi+1,j-1/2

vi+1,j+1/2

vi+1,j+3/2

vi,j+1/2

vi,j+3/2

horizontal grid vertical grid

{ Δy

Δx

{

η

ui-1/2,k+1 θi,k+1 ui+1/2,k+1 θi+1,k+1 ui+3/2

ui-1/2,k θi,k ui+1/2,k θi+1,k ui+3/2,

wi,k-1/2 wi+1,k-1/2

wi+1,k+3/2

wi,k+1/2 wi+1,k+1/2

wi,k+3/2

Δηk

Δx

{

Δηk+1/2

{

{

Figure 3.2: Horizontal and vertical grids of the ARW

Using these definitions, we can write the spatially discrete acoustic step equations

(3.7) – (3.12) as

δτU

00 + μt

∗ x

αt

∗ x

δxp

00τ + (μt

∗ x

δxp ̄)α00τ

x

+(α/αd)

x

[μt

∗ x

δxφ00τ

η

+ (δxφt

∗ η

)(δηp

00x

η

− μ00x)

τ

] = R

t

∗

U

(3.21)

δτV

00 + μt

∗ y

αt

∗ y

δyp

00τ + (μt

∗ y

δyp ̄)α00τ

y

+(α/αd)

y

[μt

∗ y

δyφ00τ

η

+ (δyφt

∗ η

)(δηp

00y

η

− μ00y)

τ

] = R

t

∗

V

(3.22)

δτμ

00

d + m2

[δxU

00 + δyV00]τ+∆τ + mδηΩ00τ+∆τ = Rt

∗

μ

(3.23)

δτΘ

00 + m2

[δx(U

00θ

t

∗ x

) + δy(V

00θ

t

∗ y

)]τ+∆τ + mδη(Ω00τ+∆τθ

t

∗ η

) = RΘ

t

∗

(3.24)

δτW00 − m−1g



(α/αd)

t

∗

η

δη(Cδηφ00) + δη



c

2

s

αt

∗

Θ00

Θt

∗



− μ

00

d

τ

= RW

t

∗

(3.25)

δτφ

00 +

1

μ

t

∗

d

[mΩ

00τ+∆τ

δηφ − gW00τ

] = Rφ

t

∗

, (3.26)

where the discrete operator

δxa = ∆x

−1

(ai+1/2 − ai−1/2) (3.27)

with the operators δy and δη similarily defined. Additionally, the operator a

η

is a vertical

interpolation operator. Using the notation given for the vertically stretched grid depicted in

Fig. 3.2, it is defined as

a

η

|k+1/2 =

1

2



∆ηk

∆ηk+1/2

ak+1 +

∆ηk+1

∆ηk+1/2

ak



. (3.28)

17
The operator vertically interpolates variables on mass levels k to the w levels (k +

1

2

). It should

be noted that the vertical grid is defined such that vertical interpolation from w levels to mass

levels reduces to a

η

k = (ak+1/2 + ak−1/2)/2 (see Fig. 3.2).

The RHS terms in the discrete acoustic step equations for momentum (3.21), (3.22) and

(3.25) are discretized as

R

t

∗

U = − (μd

xαx

δxp

0 − μdxα0

x

δxp ̄) − (α/αd)

x

(μd

x

δxφ0

η

− δηp

0

x

η

δxφ

η

+ μ

0

d

x

δxφ

η

)

+ FUcor + advection + mixing + physics (3.29)

R

t

∗

V = − (μd

yαy

δyp

0 − μdyα0

y

δyp ̄) − (α/αd)

y

(μd

y

δyφ0

η

− δηp

0

y

η

δyφ

η

+ μ

0

d

y

δyφ

η

)

+ FVcor + advection + mixing + physics (3.30)

R

t

∗

W = m−1

g(α/αd)

η

[δηp

0 +  ̄μdqmη

] − m−1μ

0

d

g

+ FWcor + advection + mixing + buoyancy + physics. (3.31)

3.2.2 Coriolis and Curvature Terms

The terms FUcor, FVcor, and FWcor in (3.29) – (3.31) represent Coriolis and curvature effects

in the equations. These terms in continuous form are given in (2.32) – (2.34). Their spatial

discretization is

FUcor = +f

x

+ u

x

δym − v

y

δxm

x

V

xy

− e

xW

xη

cos αr

x −

uW

xη

re

FVcor = −



f

y

+ u

x

δym − v

y

δxm

y

U

xy

+ e

yW

yη sin αr

y

−

vW

yη

re

FWcor = +e(U

xη

cos αr − V

yη sin αr) + 

u

xηU

xη

+ v

yηV

yη

re



.

Here the operators ()xy

= ()x

y

, and likewise for ()xη and ()yη.

3.2.3 Advection

The advection terms in the ARW solver are in the form of a flux divergence and are a subset of

the RHS terms in equations (3.13) – (3.18):

R

t

∗

Uadv = − m[∂x(Uu) + ∂y(V u)] + ∂η(Ωu) (3.32)

R

t

∗

Vadv = − m[∂x(Uv) + ∂y(V v)] + ∂η(Ωv) (3.33)

R

t

∗

μadv = − m2

[Ux + Vy] + mΩη (3.34)

R

t

∗

Θadv = − m2

[∂x(Uθ) + ∂y(V θ)] − m∂η(Ωθ) (3.35)

R

t

∗

Wadv = − m[∂x(Uw) + ∂y(V w)] + ∂η(Ωw) (3.36)

R

t

∗

φadv = − μ

−1

d

[m2

(Uφx + V φy) + mΩφη]. (3.37)

For the mass conservation equation, the flux divergence is discretized using a 2nd-order centered

approximation:

R

t

∗

μadv = −m2

[δxU + δyV ]

t

∗

+ mδηΩ

t

∗

. (3.38)

18
In the current version of the ARW, the advection of vector quantities (momentum) and scalars

is performed using the RK3 time integration as outlined in Fig. 3.1. The spatial discretization

used in this approach is outlined in the next section. For many applications it is desirable to

use positive definite or monotonic advection schemes for scalar transport. In the next major

release of the ARW we will be including a forward-in-time scheme for scalar transport that has

positive definite and monotonic options. We describe that scheme in the section following the

description of the RK3 advection.

RK3 Advection

2

nd through 6th order accurate spatial discretizations of the flux divergence are available in the

ARW for momentum, scalars and geopotential using the RK3 time-integration scheme (scalar

advection option 1, step 7 in the time-split integration sequence in Fig. 3.1). The discrete

operators can be illustrated by considering the flux divergence equation for a scalar q in its

discrete form:

R

t

∗

qadv = −m2

[δx(Uq

xadv ) + δy(V qyadv )] − mδη(Ωqηadv ). (3.39)

As in the pressure gradient discretization, the discrete operator is defined as

δx(Uq

xadv ) = ∆x−1



(Uq

xadv )i+1/2 − (Uqxadv )i−1/2



. (3.40)

The different order advection schemes correspond to different definitions for the operator q

xadv

.

The even order operators (2nd, 4th, and 6th) are

2

nd order: (qxadv )i−1/2 =

1

2

(qi + qi−1)

4

th order: (q

xadv )i−1/2 =

7

12

(qi + qi−1) −

1

12

(qi+1 + qi−2)

6

th order: (qxadv )i−1/2 =

37

60

(qi + qi−1) −

2

15

(qi+1 + qi−2) + 1

60

(qi+2 + qi−3),

and the odd order operators (3rd and 5th) are

3

rd order: (qxadv )i−1/2 = (qxadv )4

th

i−1/2

+ sign(U)

1

12



(qi+1 − qi−2) − 3(qi − qi−1)



5

th order: (q

xadv )i−1/2 = (qxadv )6

th

i−1/2

− sign(U)

1

60



(qi+2 − qi−3) − 5(qi+1 − qi−2) + 10(qi − qi−1)



.

The even-order advection operators are spatially centered and thus contain no implicit diffusion outside of the diffusion inherent in the RK3 time integration. The odd-order schemes are

upwind-biased, and the spatial discretization is inherently diffusive. As is evident in their formulation, the odd-order schemes are comprised of the next higher (even) order centered scheme

plus an upwind term that, for a constant transport mass flux, is a diffusion term of that next

higher (even) order with a hyper-viscosity proportional to the Courant number (Cr). Further

details concerning RK3 advection can be found in Wicker and Skamarock (2002)

19
Forward-In-Time Scalar Advection

A forward-in-time scalar advection scheme having monotonic and positive definite options will

be available in the next major release of the ARW. This option entails bypassing step (7) in the

time-split integration sequence in Fig. 3.1, and adding a single advection evaluation before step

(9), after the end of the RK3 loop in the integration sequence. The new advection algorithm is

patterned after Easter (1993) and can be written as follows

(μdq)

∗ = (μdq)t + m2Fx(qt

) (3.41)

μ

∗

d = μ

t

d + m2Fx(I) (3.42)

q

∗ = (μdq)∗

/(μd)

∗

(3.43)

(μdq)

∗∗ = (μdq)∗ + m2Fy(q∗

) (3.44)

μ

∗∗

d = μ

∗

d + m2Fy(I) (3.45)

q

∗∗ = (μdq)∗∗/(μd)∗∗ (3.46)

(μdq)

t+∆t = (μdq)∗∗ + mFη(q∗∗) (3.47)

μ

t+∆t

d = μ

∗∗

d + mFη(I) (3.48)

q

t+∆t = (μdq)t+∆t

/μt+∆t

d

. (3.49)

In (3.41) - (3.49) the operator

Fx(q) = −∆t∆x

−1

((U

t

q)i+1/2 − (U

t

q)i−1/2), (3.50)

with similar definitions for Fy and Fη, and I is a vector with all values equal to 1. The discrete

mass continuity equation ((3.42)+(3.45)+(3.48)) can be written as

δτμd = −m2[δxU

t

+ δyV

t

] − mδηΩ

t

. (3.51)

It can easily be seen that the scheme (3.41) - (3.49) collapses to (3.51) for q = I, and hence it is

consistent. The mass fluxes U

t

, V

t

, and Ω

t

, represent time-averaged values where the averaging

is performed over the final RK3 small-time-step cycle. Hence, the mass conservation equation

(3.51) produces the mass conservation equation integrated within the RK3 scheme. Equation

(3.51) is re-integrated within the time-split transport scheme only because a consistent column

mass μd is needed on the sub-steps to retrieve q from the prognostic variable μdq.

Any forward-in-time flux operator that is stable for Cr ≤ 1 can be used to evaluate the

operators Fx(q), Fy(q), or Fη(q). A scheme based on the Piecewise Parabolic Method (PPM)

has been implemented in the ARW. It is the non-monotonized PPM advection described in

Carpenter et al. (1990), and it provides the fluxes (U

t

q)i± 1

2

used in (3.50) in the flux divergence

operators. Defining qi as the control-volume average mixing ratio for cell i, zone edge values

q

adv can be defined as

q

adv

i+ 1

2

=



7(qi+1 + qi) − (qi+2 + qi−1)



/12. (3.52)

The flux through the i +

1

2

face can be written as

(U

t

q)i+ 1

2

= U

t

i+ 1

2



q

adv

i+ 1

2

− Cr(q

adv

i+ 1

2

− qi) − Cr(1 − Cr)(q

adv

i− 1

2

− 2qi + q

adv

i+ 1

2

)



(3.53)

20
for Cr > 0, and

(U

t

q)i+ 1

2

= U

t

i+ 1

2



q

adv

i+ 1

2

− Cr(q

adv

i+ 1

2

− qi+1) + Cr(1 + Cr)(q

adv

i+ 1

2

− 2qi+1 + q

adv

i+ 3

2

)



(3.54)

for Cr < 0, where Cr = ui+ 1

2

∆t/∆x. In addition, to maintain second-order accuracy for the

time-split forward-in-time advection scheme, the sequence (3.41) – (3.49) is reversed every other

time step; that is, the flux divergence in η is computed first followed by y and then x. Finally,

the scheme has also been augmented for Courant numbers greater than one using the 1D “semiLagrangian” flux calculation described in Lin and Rood (1996). More information about this

scheme and a description of the limiters can be found in Skamarock (2005).

3.3 Stability Constraints

There are two time steps that a user must specify when running the ARW: the model time step

(the time step used by the RK3 scheme, see Section 3.1.1) and the acoustic time step (used

in the acoustic sub-steps of the time-split integration procedure, see Section 3.1.2). Both are

limited by Courant numbers. In the following sections we describe how to choose time steps for

applications.

3.3.1 RK3 Time Step Constraint

The RK3 time step is limited by the advective Courant number u∆t/∆x and the user’s choice

of advection schemes— users can choose 2nd through 6th order discretizations for the advection

terms. The time-step limitations for 1D advection in the RK3 scheme using these advection

schemes is given in Wicker and Skamarock (2002), and is reproduced here.

Spatial order Time Scheme 3rd 4th 5th 6th

Leapfrog Unstable 0.72 Unstable 0.62

RK2 0.88 Unstable 0.30 Unstable

RK3 1.61 1.26 1.42 1.08

Table 3.1: Maximum stable Courant numbers for one-dimensional linear advection. From Wicker

and Skamarock (2002).

As is indicated in the table, the maximum stable Courant numbers for advection in the RK3

scheme are almost a factor of two greater than those for the leapfrog time-integration scheme.

For advection in three spatial dimensions, the maximum stable Courant number is 1/

√

3 times

the Courant numbers given in Table 3.1. For stability, the time step used in the ARW should

produce a maximum Courant number less than that given by theory. Thus, for 3D applications,

the time step should satisfy the following equation:

∆tmax <

Crtheory

√

3

·

∆x

umax

, (3.55)

21
where Crtheory is the Courant number taken from the RK3 entry in Table 3.1 and umax is the

maximum velocity expected in the simulation. For example in real-data applications, where jet

stream winds may reach as high as 100 ms−1, the maximum time step would be approximately

80 s on a ∆x = 10 km grid using 5th order advection. Given additional constraint from the time

splitting, and to provide a safety buffer, we usually choose a time step that is approximately

25% less than that given by (3.55). This time step is typically a factor of two greater than that

used in leapfrog-based models. For those users familiar with the MM5 model, the rule of thumb

for choosing a time step is that the time step, in seconds, should be approximately 3 times

the horizontal grid distance, in kilometers. For the ARW, the time step (in seconds) should be

approximately 6 times the grid distance (in kilometers).

3.3.2 Acoustic Time Step Constraint

The forward-backward time integration scheme used in the ARW’s 2D explicit acoustic step

integration allows a maximum Courant number Crmax = cs∆τ/∆x < 1/

√

2, where cs is the

speed of sound. We typically use a more conservative estimate for this by replacing the limiting

value 1/

√

2 with 1/2. Thus, the acoustic time step used in the model is

∆τ < 2 ·

∆x

cs

. (3.56)

For example, on a ∆x = 10 km grid, using a sound speed cs = 300 ms−1, the acoustic time

step given in (3.56) is approximately 17 s. In the ARW, the ratio of the RK3 time step to the

acoustic time step must be an even integer. For our example using a ∆x = 10 km grid in a

real-data simulation, we would specify the RK3 time step ∆t = 60s (i.e., 25% less than the 80 s

step given by (3.55), and an acoustic time step ∆τ = 15 s (i.e., 1/4 of the RK3 step, rounding

down the ∆τ = 17 s step given by (3.56)). Note that it is the ratio of the RK3 time step to the

acoustic time step that is the required input in the ARW.

22
Chapter 4

Turbulent Mixing and Model Filters

A number of formulations for turbulent mixing and filtering are available in the ARW solver.

Some of these filters are used for numerical reasons. For example, divergence damping is used to

filter acoustic modes from the solution. Other filters are meant to represent sub-grid turbulence

processes that cannot be resolved on the chosen grid. These filters remove energy from the

solution and are formulated in part on turbulence theory and observations, or represent energy

sink terms in some approximation to the Euler equation. In this section, we begin by outlining

the formulation and discretization of turbulent mixing processes in the ARW solver commonly

associated with sub-gridscale turbulence as parameterized in cloud-scale models— the secondorder horizontal and vertical mixing. In large-scale models and most NWP models, vertical

mixing is parameterized within the planetary boundary layer (PBL) physics. Vertical mixing

parameterized within the PBL physics is described later in Chapter 8. Here we note that, when

a PBL parameterization is used, all other vertical mixing is disabled. Following the outline of

turbulent mixing parameterizations in this chapter, other numerical filters available in the ARW

solver are described.

4.1 Explicit Spatial Diffusion

The ARW solver has two formulations for spatial dissipation— diffusion along coordinate surfaces and diffusion in physical (x, y, z) space. In the following sections we present the diffusion

operators for these two formulations, followed by the four separate formulations that can be used

to compute the eddy viscosities. We conclude with a description of the prognostic turbulent

kinetic energy (TKE) equation used in one set of these formulations.

4.1.1 Horizontal and Vertical Diffusion on Coordinate Surfaces

For any model variable, horizontal and vertical second order spatial filtering on model coordinate

surfaces is considered part of the RHS terms in the continuous equations (2.35) – (2.41) and can

be expressed as follows for a model variable a:

∂t(μda) = ... + μd



m∂x(mKh∂xa) + m∂y(mKh∂ya)



+ g

2

(μdα)

−1

∂η(Kvα

−1

∂ηa). (4.1)

23
For the horizontal and vertical momentum equations, (4.1) is spatially discretized as

∂tU = ... + μd

xmx



δx(mKhδxu) + δy(mxyKh

xyδyu)



+ g

2

(μd

xαx

)

−1

δη(Kv(α

xη)−1

δηu)

∂tV = ... + μd

ymy



δx(mxyKh

xyδxv) + δy(mKhδyv)



+ g

2

(μd

yα

y

)

−1

δη(Kv(α

yη)−1

δηv)

∂tW = ... + μdm



δx(mxKh

xηδxw) + δy(myKhyηδyw)



+ g

2

(μdα

η

)

−1

δη(Kvα

−1

δηw).

The spatial discretization for a scalar q, defined at the mass points, is

∂t(μdq) = ... + μdm



δx(mxP

−1

r Kh

x

δxq) + δy(myP

−1

r Kh

y

δyq)



+ g

2

(μdα)

−1

δη(Kvα

−1

δηq).

In the current ARW formulation for mixing on coordinate surfaces, the horizontal eddy viscosity

Kh is allowed to vary in space, whereas the vertical eddy viscosity does not vary in space; hence

there is no need for any spatial averaging of Kv. Additionally, note that the horizontal eddy

viscosity Kh is multiplied by the inverse turbulent Prandtl number P

−1

r

for horizontal scalar

mixing.

4.1.2 Horizontal and Vertical Diffusion in Physical Space

Coordinate Metrics

We use the geometric height coordinate in this physical space formulation. The coordinate

metrics are computed using the prognostic geopotential in the ARW solver. At the beginning

of each Runge-Kutta time step, the coordinate metrics are evaluated as part of the overall

algorithm. The definitions of the metrics are

zx = g

−1

δxφ and zy = g

−1

δyφ.

These metric terms are defined on w levels, and (zx, zy) are horizontally coincident with (u,v)

points. Additionally, the vertical diffusion terms are evaluated directly in terms of the geometric

height, avoiding the need for metric terms in the vertical.

Continuous Equations

The continuous equations for evaluating diffusion in physical space, using the velocity stress

tensor, are as follows for horizontal and vertical momentum:

∂tU = ... − m



∂xτ11 + ∂yτ12 − ∂z(zxτ11 + zyτ12)



− ∂zτ13 (4.2)

∂tV = ... − m



∂xτ12 + ∂yτ22 − ∂z(zxτ12 + zyτ22)



− ∂zτ23 (4.3)

∂tW = ... − m



∂xτ13 + ∂yτ23 − ∂z(zxτ13 + zyτ23)



− ∂zτ33. (4.4)

The stress tensor τ can be written as follows:

τ11 = −μdKhD11

τ22 = −μdKhD22

τ33 = −μdKvD33

τ12 = −μdKhD12

τ13 = −μdKhD13

τ23 = −μdKhD23.

24
Symmetry sets the remaining tensor values; τ21 = τ12, τ31 = τ13, and τ32 = τ23. The stress tensor

τ is calculated from the deformation tensor D. The continuous deformation tensor is defined as

D11 = 2 m2



∂x(m−1u) − zx∂z(m−1u)



D22 = 2 m2



∂y(m−1v) − zy∂z(m−1v)



D33 = 2 ∂zw

D12 = m2



∂y(m−1u) − zy∂z(m−1u) + ∂x(m−1v) − zx∂z(m−1v)



D13 = m2



∂x(m−1w) − zx∂z(m−1w)



+∂z(u)

D23 = m2



∂y(m−1w) − zy∂z(m−1w)



+∂z(v).

The deformation tensor is symmetric, hence D21 = D12, D31 = D13, and D32 = D23.

The diffusion formulation for scalars is

∂t(μdq) = ... +



m



∂x − ∂zzx

μdmK(∂x − zx∂z)

+

m



∂y − ∂zzy

μdmK(∂y − zy∂z)

+ ∂zμdK∂z



q. (4.5)

Spatial Discretization

Using the definition of the stress tensor, the spatial discretization of the ARW physical-space

diffusion operators for the horizontal and vertical momentum equations (4.2) - (4.4) are

∂tU = ... − mx



δxτ11 + δyτ12 − δz(zxτ11

xη + zyxyτ12yη)



− δzτ13

∂tV = ... − my



δyτ22 + δxτ12 − δz(zyτ22

yη + zxxyτ12xη)



− δzτ23

∂tW = ... − m



δxτ13 + δyτ23 − δz(zx

xητ13xη + zyyητ23yη)



− δzτ33.

The discrete forms of the stress tensor and deformation tensor are

τ11 = −μdKhD11

τ22 = −μdKhD22

τ33 = −μdKvD33

τ12 = −μd

xyKh

xyD12

τ13 = −μd

xKh

xηD13

τ23 = −μd

yKh

yηD23,

and

D11 = 2 m2



δx(m−1

x

u) − zx

xηδz(m−1

x

u)



D22 = 2 m2



δy(m−1

y

v) − zy

xηδz(m−1

y

v)



D33 = 2 δzw

D12 = (mxy)

2



δy(m−1

x

u) − zy

xηδz(m−1

x

u)

yη

+ δx(m−1

y

v) − zx

yηδz(m−1

y

v)

xη

D13 = m2



δx(m−1w) − zxδz(m−1w)

xη

+δzu

D23 = m2



δy(m−1w) − zyδz(m−1w)

yη

+δzv.

25
The spatial discretization for the scalar diffusion (4.5) is

∂t(μdq) = ... + m



δx



μd

xH1(q)



− μdδz



zx

xH1(q)

xη

+ m



δy



μd

yH2(q)



− μdδz



zy

yH2(q)

yη

+ μdδz



Kv

η

δzq



,

where

H1(q) = mxKh

x

δxq − zxδz(q

xη)



,

H2(q) = myKh

y

δyq − zyδz(q

yη)



.

4.1.3 Computation of the Eddy Viscosities

There are four options for determining the eddy viscosities Kh and Kv in the ARW solver.

External specification of Kh and Kv

Constant values for Kh and Kv can be input in the ARW namelist.

Kh determined from the horizontal deformation

The horizontal eddy viscosity Kh can be determined from the horizontal deformation using a

Smagorinsky first-order closure approach. In this formulation, the eddy viscosity is defined and

discretized as

Kh = C

2

s

l

2



0.25(D11 − D22)

2 + D2

12

xy 1

2

.

The deformation tensor components have been defined in the previous section. The length scale

l = (∆x∆y)

1/2 and Cs is a constant with a typical value Cs = 0.25. For scalar mixing, the

eddy viscosity is divided by the turbulent Prandtl number Pr that typically has a value of 1/3

(Deardorff, 1972). This option is most often used with a planetary boundary layer scheme that

independently handles the vertical mixing.

3D Smagorinsky Closure

The horizontal and vertical eddy viscosities can be determined using a 3D Smagorinsky turbulence closure. This closure specifies the eddy viscosities as

Kh,v = C

2

s

l

2

h,v max

0., D2 −



P

−1

r N

2

1/2



, (4.6)

where

D =

1

2



D

2

11 + D

2

22 + D

2

33

+



D12

xy2

+



D13

xη2

+



D23

yη2

,

and N is the Brunt-V ̈ais ̈al ̈a frequency; the computation of N, including moisture effects, is

outlined in Section 4.1.4.

26
If ∆x is less than the user-specified critical length scale lcr, then the length scale used for

calculating both Kh and Kv in (4.6) is lh,v = (∆x∆y∆z)

1/3

(and Kh = Kv = K). If ∆x is

greater than an critical length scale lcr, then the length scale lh =

√

∆x∆y in the calculation

of the horizontal eddy viscosity Kh using (4.6), and lv = ∆z for the calculation of the vertical

eddy viscosity Kv using (4.6).

Additionally, the eddy viscosities for scalar mixing are divided by the turbulent Prandtl

number Pr = 1/3.

Prognostic TKE Closure

For the predicted turbulent kinetic energy option (TKE; see section 4.1.4), the eddy viscosities

are computed using

Kh,v = Cklh,v √e,

where e is the turbulent kinetic energy (a prognostic variable in this scheme), Ck is a constant

(typically 0.15 < Ck < 0.25), and l is a length scale.

If the grid spacing ∆x is less than the critical length scale lcr, then

lh,v = min(∆x∆y∆z)

1/3

, 0.76√

e/Nfor N

2 > 0,

lh,v = (∆x∆y∆z)

1/3

for N

2 ≤ 0

(see section 4.1.4 for the calculation of N2

). Both the horizontal and vertical eddy viscosities

are multiplied by an inverse turbulent Prandtl number P

−1

r = 1 + 2l/(∆x∆y∆z)

1/3

for scalar

mixing. In this case (lcr > ∆x) the horizontal and vertical eddy viscosities are equivalent.

If the grid spacing ∆x is greater than the critical length scale lcr, then lh =

√

∆x∆y for the

calculation of Kh. For calculating Kv,

lv = min∆z, 0.76√

e/Nfor N

2 > 0,

lv = ∆z for N

2 ≤ 0.

The eddy viscosity used for mixing scalars is divided by a turbulent Prandtl number Pr. The

Prandtl number is 1/3 for the horizontal eddy viscosity Kh, and P

−1

r = 1+2l/∆z for the vertical

eddy viscosity Kv.

4.1.4 TKE equation for the 1.5 Order Turbulence Closure

The prognostic equation governing the evolution of the turbulent kinetic energy e is

∂t(μde) + (∇ · Ve)η = μd( shear production + buoyancy + dissipation ). (4.7)

The time integration and the transport terms in (4.7) are integrated in the same manner as for

other scalars (as described in Chapter 3). The right-hand side source and sink terms for e are

given as follows.

Shear Production

The shear production term in (4.7) can be written as

shear production = KhD

2

11 + KhD

2

22 + KvD

2

33 + KhD2

12

xy

+ KvD2

13

xη

+ KvD2

23

yη

.

27
Buoyancy

The buoyancy term in the TKE equation (4.7) is written as

buoyancy = −KvN

2

,

where the Brunt-V ̈ais ̈al ̈a frequency N is computed using either the formula for a moist saturated

or unsaturated environment:

N

2 = g



A

∂θe

∂z −

∂qw

∂z 

if qv ≥ qs or qc ≥ 0.01 g/Kg;

N

2 = g



1

θ

∂θ

∂z + 1.61

∂qv

∂z −

∂qw

∂z 

if qv < qs or qc < 0.01 g/Kg.

The coefficient A is defined as

A = θ

−1

1 + 1.61Lqv

RdT

1 + L2qv

CpRvT2

,

where qw represents the total water (vapor + all liquid species + all ice species), L is the latent

heat of condensation and  is the molecular weight of water over the molecular weight of dry

air. θe is the equivalent potential temperature and is defined as

θe = θ



1 +

Lqvs

CpT



,

where qvs is the saturation vapor mixing ratio.

Dissipation

If ∆x is less than the critical length scale lcr, the dissipation term in (4.7) is

dissipation = −

Ce3/2

l

,

where

C = 1.9 Ck +

(0.93 − 1.9 Ck)l

∆s

,

∆s = (∆x∆y∆z)

1/3

, and

l = min(∆x∆y∆z)

1/3

, 0.76√

e/N.

If ∆x is greater than the critical length scale lcr, the dissipation term in (4.7) is

dissipation = −

2

√

2

15

e

3/2

l

,

where

l =

kz

1 + kz/l0

,

l0 = min

αb

R zi

0

√

ez dz

R zi

0

√

e dz

, 80

,

αb = 0.2, and k = 0.4 is the von Karman constant.

28
4.2 Filters for the Time-split RK3 scheme

Three filters are used in the ARW time-split RK3 scheme apart from those in the model physics:

three-dimensional divergence damping (an acoustic model filter); an external-mode filter that

damps vertically-integrated horizontal divergence; and off-centering of the vertically implicit

integration of the vertical momentum equation and geopoential equation. Each of these is

described in the following sections.

4.2.1 Divergence Damping

The damping of the full mass divergence is a filter for acoustic modes in the ARW solver. This

3D mass divergence damping is described in Skamarock and Klemp (1992). The filtering is

accomplished by forward weighting the pressure update in the acoustic step loop described in

Section 3.1.3, step (6). The linearized equation of state (3.5) is used to diagnose the pressure at

the new time τ after the U

00

, V

00

, μ

00

d

, and Θ00 have been advanced. Divergence damping consists

of using a modified pressure in the computation of the pressure gradient terms in the horizontal

momentum equations in the acoustic steps (Equations (3.7) and (3.8)). Denoting the updated

value as p

τ

, the modified pressure p

∗τ used in (3.7) and (3.8) can be written as

p

∗τ = pτ + γd



p

τ − pτ−∆τ



, (4.8)

where γd is the damping coefficient. This modification is equivalent to inserting a horizontal

diffusion term into the equation for the 3D mass divergence, hence the name divergence damping.

A divergence damping coefficient γd = 0.1 is typically used in the ARW applications, independent

of time step or grid size.

4.2.2 External Mode Filter

The external modes in the solution are damped by filtering the vertically-integrated horizontal

divergence. This damping is accomplished by adding a term to the horizontal momentum

equations. The additional term added to (3.7) and (3.8) are

δτU

00 = ... − γeδx(δτ−∆τμ00

d

) (4.9)

and

δτV

00 = ... − γeδy(δτ−∆τμ00

d

). (4.10)

The quantity δτ−∆τμ is the vertically-integrated mass divergence (i.e., (3.20)) from the previous

acoustic step (that is computed using the time τ values of U and V ), and γe is the external

mode damping coefficient. An external mode damping coefficient γe = 0.01 is typically used in

the ARW applications, independent of time step or grid size.

4.2.3 Semi-Implicit Acoustic Step Off-centering

Forward-in-time weighting of the vertically-implicit acoustic-time-step terms damps instabilities

associated vertically-propagating sound waves. The forward weighting also damps instabilities

associated with sloping mode levels and horizontally propagating sound waves (see Durran and

29
Klemp, 1983; Dudhia, 1995). The off-centering is accomplished by using a positive (non-zero)

coefficient β (3.19) in the acoustic time-step vertical momentum equation (3.11) and geopotential

equation (3.12). An off-centering coefficient β = 0.1 is typically used in the ARW applications,

independent of time step or grid size.

4.3 Other Damping

4.3.1 Gravity Wave Absorbing Layer

A gravity-wave absorbing layer is available in the ARW solver. The absorbing layer increases the

second-order horizontal and vertical eddy viscosities in the absorbing layer using the following

formulation:

Kdh =

∆x

2

∆t

γg cos

ztop − z

zd

π

2



,

and

Kdv =

∆z

2

∆t

γg cos

ztop − z

zd

π

2



.

Here γg is a user-specified damping coefficient, ztop is the height of the model top for a particular

grid column, zd is the depth of the damping layer (from the model top), and Kdh and Kdv are

the horizontal and vertical eddy viscosities for the wave absorbing layer. If other spatial filters

are being used, then the eddy viscosities that are used for the second-order horizontal and

vertical eddy viscosities are the maximum of (Kh, Kdh) and (Kv, Kdv). The effect of this filter

on gravity waves is discussed in Klemp and Lilly (1978), where guidance on the choice of the

damping coefficeint γg can also be found.

4.3.2 Rayleigh Damping Layer

A Rayleigh damping layer is also available in the ARW solver. This scheme applies a tendency

to u, v, w, and θ to gradually relax the variable back to a predetermined reference state value,

∂u

∂t = τ (z) (u − u),

∂v

∂t = τ (z) (v − v),

∂w

∂t = τ (z)w,

∂θ

∂t = τ (z)



θ − θ



.

Overbars indicate the reference state fields, which are functions of z only and are typically

defined as the initial horizontally homogeneous fields in idealized simulations. The reference

state vertical velocity is assumed to be zero. The variable τ defines the vertical structure of the

damping layer, and has a form similar to the Rayleigh damper in Durran and Klemp (1983),

τ (z) = (

−γr sin2

h

π

2



1 −

ztop−z

zd

i for z ≥ (ztop − zd);

0 otherwise,

30
where γr is a user specified damping coeficient, ztop is the height of the model top for a particular

grid column, and zd is the depth of the damping layer (from the model top).

Because the model surfaces change height with time in the ARW solver, the reference state

values at each grid point need to be recalculated at every time step. Thus, a linear interpolation

scheme is used to calculate updated reference state values based on the height of the model

levels at each time step.

The effect of this filter on gravity waves is discussed in Klemp and Lilly (1978), where

guidance on the choice of the damping coefficeint γr can also be found.

4.3.3 Vertical-Velocity Damping

This is also called w-damping. In semi-operational or operational NWP applications, the model

robustness can be improved by detecting locations where the vertical motion approaches the

limiting Courant number for stability, and applying a Rayleigh damping term in the vertical

momentum equation to stabilize the motion. This term is non-physical and should only be used

in the situation where many, or long, model runs are being done, and there is no option for a

re-run with a shorter time-step if a failure occurs due to an excessively strong updraft. This

might be the case, for example, in an operational setting where real-time forecasts have to be

produced on time to be useful. However, if this term activates frequently, consideration should

be given to reducing the model time-step.

The term is calculated from

Cr =









Ωdt

μdη









.

If Cr > Crβ, then

∂tW = ... − μd sign(W)γw(Cr − Crβ),

where γw is the damping coefficient (typically 0.3 ms−2), and Crβ is the activation Courant

number (typically 1.0). The ARW outputs the location of this damping when it is active.

31
32
Chapter 5

Initial Conditions

The ARW may be run with initial conditions that are defined analytically for idealized simulations, or it may be run using interpolated data from either a large-scale analysis or forecast

for real-data cases. Both 2D and 3D tests cases for idealized simulations are provided. Several

sample cases for real-data simulations are provided, which rely on pre-processing from an external package that converts the large-scale GriB data into a format suitable for ingest by the

ARW’s real-data processor.

The programs that generate the specific initial conditions for the selected idealized or realdata case function similarly. They provide the ARW with:

• input data that is on the correct horizontal and vertical staggering;

• hydrostatically balanced reference state and perturbation fields; and

• metadata specifying such information as the date, grid physical characteristics, and projection details.

For neither the idealized nor the real-data cases are the initial conditions enhanced with observations. However, output from the ARW system initial condition programs is suitable as input

to the WRF variational assimilation package (see Chapter 9).

5.1 Initialization for Idealized Conditions

The ARW comes with a number of test cases using idealized environments, including mountain waves (em hill2d x), squall lines (em squall2d x, em squall2d y), supercell thunderstorms

(em quarter ss), gravity currents (em grav2d x), and baroclinic waves (em b wave). A brief description of these test cases can be found in the README test cases file provided in the ARW

release. The test cases include examples of atmospheric flows at fine scales (e.g., the gravity

current example has a grid-spacing of 100 meters and a time step of 1 second) and examples of

flow at large scales (e.g., the baroclinic wave test case uses a grid-spacing of 100 km and a time

step of 600 s), in addition to the traditional mesoscale and cloudscale model simulations. The

test suite allows an ARW user to easily reproduce these known solutions. The test suite is also

the starting point for constructing idealized flow simulations by modifying initializations that

closely resemble a desired initialization.

All of these tests use as input a 1D sounding specified as a function of geometric height z

(except for the baroclinic wave case that uses a 2D sounding specified in [y, z]), and, with the

33
exception of the baroclinic wave test case, the sounding files are in text format that can be

directly edited by the user. The 1D specification of the sounding in these test files requires the

surface values of pressure, potential temperature, and water vapor mixing ratio, followed by the

potential temperature, vapor mixing ratio, and horizontal wind components at some heights

above the surface. The initialization programs for each case assume that this moist sounding

represents an atmosphere in hydrostatic balance.

Two sets of thermodynamic fields are needed for the model— the reference state and the

perturbation state (see Chapter 2 for further discussion of the equations). The reference state

used in the idealized initializations is computed using the input sounding from which the moisture is discarded (because the reference state is dry). The perturbation state is computed using

the full moist input sounding. The procedure for computing the hydrostatically-balanced ARW

reference and perturbation state variables from the input sounding is as follows. First, density

and both a dry and full hydrostatic pressure are computed from the input sounding at the input

sounding levels z. This is accomplished by integrating the hydrostatic equation vertically up the

column using the surface pressure and potential temperature as the lower boundary condition.

The hydrostatic equation is

δzp = −ρ

z

g(1 + (Rd/Rv)qv

z

), (5.1)

where ρ

z

is a two point average between input sounding levels, and δzp is the difference of

the pressure between input sounding levels divided by the height difference. Additionally, the

equation of state is needed to close the system:

αd =

1

ρd

=

Rdθ

po



1 +

Rd

Rv

qv

 p

po

−

cv

cp

, (5.2)

where qv and θ are given in the input sounding. (5.1) and (5.2) are a coupled set of nonlinear

equations for p and ρ in the vertical integration, and they are solved by iteration. The dry

pressure on input sounding levels is computed by integrating the hydrostatic relation down from

the top, excluding the vapor component.

Having computed the full pressure (dry plus vapor) and dry air pressure on the input sounding levels, the pressure at the model top (pdht) is computed by linear interpolation in height

(or possibly extrapolation) given the height of the model top (an input variable). The column

mass μd is computed by interpolating the dry air pressure to the surface and subtracting it from

pdht. Given the column mass, the dry-air pressure at each η level is known from the coordinate

definition (2.1), repeated here

η = (pdh − pdht)/μd where μd = pdhs − pdht,

and the pressures pdhs and pdht refer to the dry atmosphere. The potential temperature from

the input sounding is interpolated to each of the model pressure levels, and the equation of state

(5.2) is used to compute the inverse density αd. Finally, the ARW’s hydrostatic relation (2.9),

in discrete form

δηφ = −αdμd

is used to compute the geopotential. This procedure is used to compute the reference state

(assuming a dry atmosphere) and the full state (using the full moist sounding). The perturbation

variables are computed as the difference between the reference and full state. It should also be

34
noted that in the nonhydrostatic model integration, the inverse density αd is diagnosed from

the geopotential using this equation of state, and the pressure is diagnosed from the equation of

state using the inverse density αd and the prognostic potential temperature θ. Thus, the ARW’s

prognostic variables μd, θ, and φ are in exact hydrostatic balance for the model equations (to

machine roundoff).

5.2 Initialization for Real-Data Conditions

The initial conditions for the real-data cases are pre-processed through a separate package called

the Standard Initialization (SI, see Fig. 5.1). The output from the SI is passed to the realdata pre-processor in the ARW— program real— which generates initial and lateral boundary

conditions. This section is primarily about the steps taken to build the initial and the lateral

boundary conditions for a real-data case. Even though the SI is outside of the ARW system, a

brief description is appropriate to see how the raw meteorological and static terrestrial data are

brought into the model for real-data cases.

5.2.1 Use of the Standard Initialization by the ARW

The SI is a set of programs that takes terrestrial and meteorological data (typically in GriB

format) and transforms them for input to the ARW pre-processor program for real-data cases

(real). Figure 5.1 shows the flow of data into and out of the SI system. The first step for the

SI is to define a physical grid (including the projection type, location on the globe, number of

grid points, nest locations, and grid distances) and to interpolate static fields to the prescribed

domain. Independent of the domain configuration, an external analysis or forecast is processed

by the SI’s GriB decoder, which diagnoses required fields and reformats the GriB data into an

internal binary format. With a specified domain, the SI horizontally interpolates the meteorological data onto the projected domain(s), and vertically interpolates the data to the ARW’s

η coordinate. The output data from the SI supplies a complete 3-dimensional snapshot of the

atmosphere on the selected model grid’s horizontal and vertical staggering at the selected time

slices, which is sent to the ARW pre-processor program for real-data cases.

The input to the ARW real-data processor from the SI contains 3-dimensional fields of potential temperature (K), mixing ratio (kg/kg), and the horizontal components of momentum (m/s,

already rotated to the model projection). The 2-dimensional static terrestrial fields include:

albedo, Coriolis parameters, terrain elevation, vegetation/land-use type, land/water mask, map

scale factors, map rotation angle, soil texture category, vegetation greenness fraction, annual

mean temperature, and latitude/longitude. The 2-dimensional time-dependent fields from the

external model, after processing by the SI, include: μd (Pa), layers of soil temperature (K) and

soil moisture (kg/kg, either total moisture, or binned into total and liquid content), snow depth

(m), skin temperature (K), and fractional sea ice. All of the fields in the final output from the

SI are on the correct horizontal and vertical staggering for the ARW. The input data from the

SI is assumed to be hydrostatically balanced.

35
GRID_GEN

GRIB_PREP

HINTERP

VINTERP

STATIC

DATA

GRIB

DATA

SI

SYSTEM

GRID_GEN

GRIB_PREP

HINTERP

VINTERP

STATIC

DATA

GRIB

DATA

SI

SYSTEM

REAL

WRF

REAL-DATA

ARW

SYSTEM

REAL

WRF

Figure 5.1: Schematic showing the data flow and program components in the SI, and how the

SI feeds initial data to the ARW. Letters in the rectangular boxes indicate program names.

GRID GEN: defines the model domain and creates static files of terrestrial data. GRIB PREP:

decodes GriB data. HINTERP: interpolates meteorological data to the model domain. VINTERP: vertically interpolates data to model coordinate.

5.2.2 Reference and Perturbation State

Identical to the idealized initializations, there is a partitioning of some of the meteorological

data into reference and perturbation fields. For real-data cases, the reference state is defined by

terrain elevation and the following three constants:

• p0 (105 Pa) reference sea level pressure;

• T0 (usually 270 to 300 K) reference sea level temperature; and

• A (50 K) temperature difference between the pressure levels of p0 and p0/e.

Using these parameters, the dry reference state surface pressure is

pdhs = p0 exp 

−T0

A

+

s

T0

A

2

−

2φsfc

A Rd

!

. (5.3)

From (5.3), the 3-dimensional reference pressure is computed as a function of the vertical coordinate η levels and the model top pdht (input provided by SI for real-data cases):

pd = η (pdhs − pdht) + pdht. (5.4)

With (5.4), the reference temperature is a straight line on a skew-T plot, defined as

T = T0 + A lnpd

p0

.

36
From the reference temperature and pressure, the reference potential temperature is then defined

as

θd =



T0 + A lnpd

p0

p0

pd

 Rd

Cp

, (5.5)

and the reciprocal of the reference density using (5.4) and (5.5) is given by

αd =

1

ρd

=

Rd θd

p0



pd

p0

−

Cv

Cp

. (5.6)

The base state difference of the dry surface pressure from (5.3) and the model top is given as

μd = pdhs − pdht. (5.7)

From (5.6) and (5.7), the reference state geopotential defined from the hydrostatic relation is

δηφ = −αd μd.

One of the total fields provided to the real-data cases by the SI is μd. The perturbation field

given the reference value (5.7) is

μ

0

d = μd − μd

. (5.8)

Starting with the reference state fields (5.4, 5.6, and 5.7) and the dry surface pressure perturbation (5.8), the perturbation fields for pressure and inverse density are diagnosed. The pressure

perturbation includes moisture and is diagnosed from the hydrostatic equation

δηp

0 = μ0

d



1 + qv

η



+ qv

η μd

,

which is integrated down from at the model top (where p

0 = 0) to recover p0

. The total dry

inverse density is given as

αd =

Rd

p0

θ



1 +

Rv

Rd

qv

 p

0

d + pd

p0

−

Cv

Cp

,

which defines the perturbation field for inverse density

α

0

d = αd − αd.

The perturbation geopotential is diagnosed from the hydrostatic relation

δηφ

0 = −



μdα

0

d + μ

0

dαd



by upward integration using the terrain elevation as the lower boundary condition. In future

versions of the real-data pre-processor, p

0 will be re-diagnosed consistent with the method used

in the model (2.21) as a final step. No modifications to the original μd, u, v, qv, or θ from the

SI are performed. The vertical component of velocity is initialized to zero.

37
5.2.3 Generating Lateral Boundary Data

This section deals with generating the separate lateral boundary condition file used exclusively

for the real-data cases. For information on which lateral boundary options are available for both

the idealized and real-data cases, see Chapter (6).

The specified lateral boundary condition for the coarse grid for real-data cases is supplied

by an external file that is generated by program real. This file contains records for the fields u,

v, θ, qv, φ

0

, and μ

0

d

that are used by the ARW to constrain the lateral boundaries (other fields

are in the boundary file, but are not used). The lateral boundary file holds one less time period

than was processed by the SI. Each of these variables has both a valid value at the initial time

of the lateral boundary time and a tendency term to get to the next boundary time period. For

example, assuming a 3-hourly availability of data from the SI, the first time period of the lateral

boundary file for u would contain data for both coupled u (map scale factor and μd interpolated

to the variable’s staggering) at the 0 h time

U0h =

μd

xu

mx









0h

,

and a tendency value defined as

Ut =

U3h − U0h

3h

,

which would take a grid point from the initial value to the value at the next large-scale time

during 3 simulation hours. The horizontal momentum fields are coupled both with μd and the

inverse map factor. The other 3-dimensional fields (θ, qv, and φ

0

) are coupled only with μd. The

μ

0

d

lateral boundary field is not coupled.

Each lateral boundary field is defined along the four sides of the rectangular grid (loosely

referred to as the north, south, east, and west sides). The boundary values and tendencies for

vertical velocity and the non-vapor moisture species are included in the external lateral boundary

file, but act as place-holders for the nested boundary data for the fine grids. The width of the

lateral boundary along each of the four sides is user selectable at run-time.

5.2.4 Masking of Surface Fields

Some of the meteorological and static fields are “masked”. A masked field is one in which the

values are typically defined only over water (e.g., sea surface temperature) or defined only over

land (e.g., soil temperature). The need to match all of the masked fields consistently to each

other requires additional steps for the real-data cases due to the masked data’s presumed use in

various physics packages in the soil, at the surface, and in the boundary layer. If the land/water

mask for a location is flagged as a water point, then the vegetation and soil categories must also

recognize the location as the special water flag for each of their respective categorical indices.

The values for the soil temperature and soil moisture come from the SI on the native levels

originally defined for those variables in the large-scale model. The SI does no vertical interpolation for the soil data. While it is typical to try to match the ARW soil scheme with the

incoming data, that is not a requirement. Pre-processor real will vertically interpolate (linear in

depth below the ground) from the incoming levels to the requested soil layers to be used within

the model.

38
Chapter 6

Lateral Boundary Conditions

Several lateral boundary condition options exist for the ARW that are suitable for idealized

flows, and a specified lateral boundary condition for real-data simulations is available. These

choices are handled via a run-time option in the Fortran namelist file. The coarsest grid of

any single simulation is eligible for any of the lateral boundary selections. For example, realdata cases could use combinations of periodic, symmetric, or open lateral boundary conditions

instead of the more traditional time-dependent conditions provided by an external boundary

file. However, use of the specified time-dependent lateral boundary conditions for one of the

idealized simulations is not possible because an external boundary file is not generated. The

ARW supports rectangular horizontal grid refinement with integer ratios of the parent and child

grid distances and time steps.

For nesting, all fine domains use the nest time-dependent lateral boundary condition where

the outer row and column of the fine grid is specified from the parent domain, described in

Section 7.3.

6.1 Periodic Lateral Boundary Conditions

Periodic lateral boundary conditions in the ARW can be specified as periodic in x (west-east),

y (south-north), or doubly periodic in (x, y). The periodic boundary conditions constrain the

solutions to be periodic; that is, a generic model state variable ψ will follow the relation

ψ(x + nLx, y + mLy) = ψ(x, y)

for all integer (n, m). The periodicity lengths (Lx, Ly) are

[(dimension of the domain in x) - 1]∆x and [(dimension of the domain in y) - 1]∆y.

6.2 Open Lateral Boundary Conditions

Open lateral boundary conditions, also referred to as gravity-wave radiating boundary conditions, can be specified for the west, east, north, or south boundary, or any combination thereof.

The gravity wave radiation conditions follow the approach of Klemp and Lilly (1978) and Klemp

and Wilhelmson (1978).

There are a number of changes in the base numerical algorithm described in Chapter 3 that

accompany the imposition of these conditions. First, for the normal horizontal velocities along

39
a boundary on which the condition is specified, the momentum equation for the horizontal

velocity, (3.7) or (3.8), is replaced by

δτU

00 = −U∗

δxu,

where U

∗ = min(U − cb, 0) at the x = 0 (western) boundary, U∗ = max(U + cb, 0) at the x = L

(eastern) boundary, and likewise for the south and north boundaries for the V momentum. The

horizontal difference operator δx is evaluated in a one-sided manner using the difference between

the boundary value and the value one grid-point into the grid from the boundary. cb is the phase

speed of the gravity waves that are to be radiated; it is specified as a model constant (for more

details see Klemp and Lilly, 1978; Klemp and Wilhelmson, 1978).

For scalars and non-normal momentum variables, the boundary-perpendicular flux divergence term is replaced with

δx(Uψ) = U

∗

δxψ + ψ δxU,

where U

∗ = min(U, 0) at the x = 0 + ∆x/2 (western) scalar boundary, U∗ = max(U, 0) at the

x = L− ∆x/2 (eastern) boundary, and likewise for the south and north boundaries using V . As

was the case for the momentum equations, the horizontal difference operator δx is evaluated in a

one-sided manner using the difference between the boundary value and the value one grid-point

into the grid from the boundary.

6.3 Symmetric Lateral Boundary Conditions

Symmetry lateral boundary conditions can be specified for the west, east, north, or south boundary, or any combination thereof. The symmetry boundaries are located on the normal-velocity

planes at the lateral edges of the grids. The normal velocities are zero at these boundaries, and

on either side of the boundary the normal velocity satisfies the relation

U⊥(xb − x) = −U⊥(xb + x),

where xb is the location of the symmetry boundary. All other variables satisfy the relation

ψ(xb − x) = ψ(xb + x).

6.4 Specified Lateral Boundary Conditions

Primarily for real-data cases, the specified boundary condition is also referred to as a relaxation,

or nudging, boundary condition. There are two uses of the specified boundaries in the ARW: for

the outer-most coarse grid or for the time-dependent boundaries supplied to a nested grid. The

specified lateral boundary conditions for the nest are automatically selected for all of the fine

grids, even if the coarse grid is using combinations of the symmetry, periodic, or open options.

If the specified lateral boundary condition is selected for the coarse grid, then all four grid sides

(west, east, north, and south) use specified lateral conditions.

The coarse grid specified lateral boundary is comprised of both a specified and a relaxation

zone as shown in Fig. 6.1). For the coarse grid, the specified zone is determined entirely by

temporal interpolation from an external forecast or analysis (supplied by the SI). The width

40
Relaxation

Rows

South

West East

North

N=1

N=2

N=3

N=4

N=5

N=5

N=4

N=3

N=2

N=1

1 2 3 4 5 5 4 3 2 1

Relaxation

Columns

Specified

Columns

Specified

Rows

Real-Data Lateral Boundary Condition: Location of Specified and Relaxation Zones

Figure 6.1: Specified and relaxation zones for a grid with a single specified row and column,

and four rows and columns for the relaxation zone. These are typical values used for a specified

lateral boundary condition for a real-data case.

of the specified zone is run-time configurable, but is typically set to 1 (i.e., the last row and

column along the outer edge of the coarse grid). The second region of the lateral boundary for

the coarse grid is the relaxation zone. The relaxation zone is where the model is nudged or

relaxed towards the large-scale forecast (e.g., rows and columns 2 through 5 in Fig. 6.1). The

size of the relaxation zone is a run-time option.

The specified lateral boundary condition for the coarse grid requires an external file, generated during the same pre-processing as the initial condition file. Let ψ be any prognostic value

having a lateral boundary entry, after Davies and Turner (1977),

∂tψ





n

= F1(ψLS − ψ) − F2∆

2

(ψLS − ψ), (6.1)

where n is the number of grid points in from the outer row or column along the boundary

(SpecZone+ 1 ≤ n ≤ SpecZone+RelaxZone−1; see Fig. 6.1) and ψLS is the large-scale value

obtained by spatial and temporal interpolation from the analyses. ∆2is a 5-point horizontal

smoother applied along η-surfaces. The weighting function coefficients F1 and F2 in (6.1) are

given by

F1 =

1

10∆t

SpecZone + RelaxZone − n

RelaxZone − 1

,

F2 =

1

50∆t

SpecZone + RelaxZone − n

RelaxZone − 1

,

where n extends only through the relaxation zone (SpecZone+1 ≤ n ≤ SpecZone+RelaxZone−

1). F1 and F2 are linear ramping functions with a maximum at the first relaxation row or column

nearest the coarse grid boundary (just inside of the specified zone).

41
On the coarse grid, the specified boundary condition applies to the horizontal wind components, potential temperature, φ

0

, μ

0

d

, and water vapor. The lateral boundary file contains

enough information to update the boundary zone values through the entire simulation period.

The momentum fields are coupled with μd and the inverse map factor (both at the specific

variable’s horizontal staggering location), and the other 3-dimensional fields are coupled with

μd. The μ

0

d field is not coupled in the lateral boundary file.

Vertical velocity has a zero gradient boundary condition applied in the specified zone (usually

the outer-most row and column). Microphysical variables, except vapor, and all other scalars

have flow-dependent boundary conditions applied in the specified zone. This boundary condition

specifies zero on inflow and zero-gradient on outflow. Since these boundary conditions require

only information from the interior of the grid, these variables are not in the specified boundary

condition file.

42
Chapter 7

Nesting

The ARW supports horizontal nesting that allows resolution to be focused over a region of

interest by introducing an additional grid (or grids) into the simulation. In the current implementation, only horizontal refinement is available: there is no vertical nesting option. The

nested grids are rectangular and are aligned with the parent (coarser) grid within which they are

nested. Additionally, the nested grids allow any integer spatial (∆xcoarse/∆xf ine) and temporal

refinements of the parent grid. This nesting implementation is in many ways similar to the

implementations in other mesoscale and cloudscale models (e.g. MM5, ARPS, COAMPS). The

major improvement in the ARW’s nesting infrastruture compared with techniques used in other

models is the ability to compute nested simulations efficiently on parallel distributed-memory

computer systems, which includes support for moving nested grids. The WRF Software Framework, described in Michalakes et al. (2004), makes these advances possible. In this chapter we

describe the various nesting options available in the ARW and the numerical coupling between

the grids.

7.1 Overview

1-Way and 2-Way Grid Nesting

Nested grid simulations can be produced using either 1-way nesting or 2-way nesting as outlined

in Fig. 7.1. The 1-way and 2-way nesting options refer to how a coarse grid and the fine grid

interact. In both the 1-way and 2-way simulation modes, the fine grid boundary conditions

(i.e., the lateral boundaries) are interpolated from the coarse grid forecast. In a 1-way nest, this

is the only information exchange between the grids (from coarse grid to fine grid). Hence, the

name 1-way nesting. In the 2-way nest integration, the fine grid solution replaces the coarse grid

solution for coarse grid points that lie inside the fine grid. This information exchange between

the grids is now in both directions (coarse-to-fine and fine-to-coarse). Hence, the name 2-way

nesting.

The 1-way nest option may be run in one of two different methods. One option is to

produce the nested simulation as two separate ARW simulations as described in the leftmost

box in Fig. 7.1. In this mode, the coarse grid is integrated first. Output from the coarse

grid integration is then processed to provide boundary conditions for the nested run (usually

at a much lower temporal frequency than the coarse grid time step), and this is followed by

the complete time integration of fine (nested) grid. Hence, this 1-way option is equivalent to

43
Concurrent ARW

simulation with two

domains

Both CG and FG

simulations run within

the same WRF

2) FG LBC from CG

at each coarse time

step

3) CG integrates one

time step, then the

FG integrates up to

the same time

4) Feedback FG to CG

2-Way ARW Simulation

Two Consecutive

ARW simulations

(using ndown)

1) Run coarse grid

(CG) simulation

2) Process CG for

initial condition (IC)

for fine grid (FG)

3) Process CG for

lateral boundary

condition (LBC) for

FG

4) Run FG simulation

Concurrent ARW

simulation with two

domains

1) Both CG and FG

simulations run within

the same WRF

2) FG LBC from CG

at each coarse time

step

3) CG integrates one

time step, then the

FG integrates up to

the same time

1-Way AWRRWF Simulation

Figure 7.1: 1-way and 2-way nesting options in the ARW.

running two separate simulations with a processing step in between. Also with separate grid

simulations, an intermediate re-analysis (such as via 3D-Var, see Section 9) can be used.

The second 1-way option (lockstep with no feedback), depicted in the middle box in Fig. 7.1,

is run as a traditional simulation with two (or more) grids integrating concurrently, except with

the feedback runtime option shut off. This option provides lateral boundary conditions to the

fine grid at each coarse grid time step, which is an advantage of the concurrent 1-way method

(no feedback).

Fine Grid Initialization Options

The ARW supports several strategies to refine a coarse-grid simulation with the introduction of

a nested grid. When using 1-way and 2-way nesting, several options for initializing the fine grid

are provided.

• All of the fine grid variables can be interpolated from the coarse grid.

• All of the fine grid variables can be input from an external file which has high-resolution

information for both the meteorological and the terrestrial fields.

• The fine grid can have some of the variables initialized with a high-resolution external

data set, while other variables are interpolated from the coarse grid.

• For a moving nest, an external orography file may be used to update the fine grid terrain

elevation. This option is not generally available in this release.

These fine grid initialization settings are user specified at run-time, and the ARW allows nested

grids to instantiate and cease during any time that the fine grid’s parent is still integrating. (The

44
(c) 1

2

3

(c) 1

2

3

(b) 1

2 3

(b) 1

2 3

(a) 1

2 3 4

(a) 1

2 3 4

(d) 1

2 4 3

(d) 1

2 4 3

x x

Figure 7.2: Various nest configurations for multiple grids. (a) Telescoping nests. (b) Nests at

the same level with respect to a parent grid. (c) Overlapping grids: not allowed (d) Inner-most

grid has more than one parent grid: not allowed

system is currently constrained to starting nests at the beginning of the coarse grid simulations

if runs require input of nest-resolution terrain or other lower boundary data. This limitation

will be addressed in the near future.)

Possible Grid Configurations

A simulation involves one outer grid and may contain multiple inner nested grids. In the ARW,

each nested region is entirely contained within a single coarser grid, referred to as the parent

grid. The finer, nested grids are referred to as child grids. Using this terminology, children

are also parents when multiple levels of nesting are used. The fine grids may be telescoped to

any depth (i.e., a parent grid may contain one or more child grids, each of which in turn may

successively contain one or more child grids; Fig. 7.2a), and several fine grids may share the

same parent at the same level of nesting (Fig. 7.2b). Any valid fine grid may either be a static

domain or it may be a moving nest with prescribed incremental shifts. The ARW does not

permit overlapping grids, where a coarse grid point is contained within more than a single child

grid (i.e., both of which are at the same nest level with respect to the parent; Fig. 7.2c). In

addition, no grid can have more than a single parent (Fig. 7.2d).

For 2-way nested grid simulations, the ratio of the parent horizontal grid distance to the

child horizontal grid distance (the spatial refinement ratio) must be an integer. This is also

true for the time steps (the temporal refinement ratio). The model does allow the time step

refinement ratio to differ from the spatial refinement ratio. Also, nested grids on the same level

(i.e., children who have the same parent) may have different spatial and temporal refinement

45
U U U U U U U

U U U U U U U

U U U U U U U

U U U U U U U

U U U U U U U

U U U U U U U

V V V V V V

V V V V V V

V V V V V V

V V V V V V

V V V V V V

V V V V V V

V V V V V V

q q q q q q

q q q q q q

q q q q q q

q q q q q q

q q q q q q

q q q q q q

U U U U U

U U

U U

U U U U U

V V V V

V V

V V

V V

V V V V

q q q q

q q

q q

q q q q

Figure 7.3: Arakawa-C grid staggering for a portion of a parent domain and an imbedded nest

domain with a 3:1 grid size ratio. The solid lines denote coarse grid cell boundaries, and the

dashed lines are the boundaries for each fine grid cell. The horizontal components of velocity

(“U” and “V”) are defined along the normal cell face, and the thermodynamic variables (“θ”)

are defined at the center of the grid cell (each square). The bold typeface variables along the

interface between the coarse and the fine grid define the locations where the specified lateral

boundaries for the nest are in effect.

ratios.

7.2 Nesting and Staggering

The ARW uses an Arakawa-C grid staggering. As shown in Fig. 7.3, the u and v components of horizontal velocity are normal to the respective faces of the grid cell, and the

mass/thermodynamic/scalar variables are located in the center of the cell.

The variable staggering has an additional column of u in the x-direction and an additional row of v in the y-direction because the normal velocity points define the grid boundaries. The horizontal momentum components reflect an average across each cell-face, while

each mass/thermodynamic/scalar variable is the representative mean value throughout the cell.

Feedback is handled to preserve these mean values: the mass/thermodynamic/scalar fields are

fed back with an average from within the entire coarse grid point (Fig. 7.3), and the horizontal

momentum variables are averaged along their respective normal coarse grid cell faces.

The horizontal interpolation (to instantiate a grid and to provide time-dependent lateral

46
boundaries) does not conserve mass. The feedback mechanism, for most of the unmasked fields,

uses cell averages (for mass/thermodynamic/scalar quantities) and cell-face averages for the

horizontal momentum fields.

The staggering defines the way that the fine grid is situated on top of the coarse grid. For all

odd ratios there is a coincident point for each variable: a location that has the coarse grid and

the fine grid at the same physical point. The location of this point depends on the variable. In

each of the coarse-grid cells with an odd ratio, the middle fine-grid cell is the coincident point

with the coarse grid for all of the mass-staggered fields (Fig. 7.3). For the horizontal momentum

variables the normal velocity has coincident points along the grid boundaries for odd ratios.

For fields that are averaged back to the coarse grid in the feedback, the mean of the nine

mass/thermodynamic/scalar (for example, due to the 3:1 grid-distance ratio in the example

shown in Fig. 7.3) fine grid points is fed back to the coarse grid. These fields include most

3D and 2D arrays. For the horizontal momentum fields averaged back to the coarse grid in

the feedback, the mean of three (for example, due to the 3:1 grid-distance ratio in the example

shown in Fig. 7.3) fine grid points is fed back to the coarse grid from along the coincident cell

face. The fields that are masked due to the land/sea category are fed back directly from the

coincident points for odd ratios. Only masked fields use the feedback method where a single

point from the fine grid is assigned to the coarse grid.

A difference between the odd and even grid-distance ratios is in the feedback from the fine

grid to the coarse grid. No coincident points exist for the single point feedback mechanisms

for even grid distance ratios (such as used for the land/sea masked 2D fields). For a 2:1 even

grid distance ratio, Figure 7.4 shows that each coarse grid point has four fine grid cells that

are equally close, and therefore four equally eligible grid points for use as the single fine-grid

point that feeds back to the coarse grid. The single-point feedback is arbitrarily chosen as the

south-west corner of the four neighboring points. This arbitrary assignment to masked fields

implies that even grid distance ratios are more suited for idealized simulations where masked

fields are less important.

7.3 Nested Lateral Boundary Conditions

For the fine grid with 2-way nesting or 1-way nesting (using a concurrent ARW simulation,

see Fig. 7.1), the boundary conditions are specified by the parent grid at every coarse-grid

time step. The nest lateral boundary condition behaves similarly to the specified boundary

condition for real-data cases (see Section 6.4), but the relaxation zone is not active. Prognostic

variables are entirely specified in the outer row and column of the fine grid through spatial and

temporal interpolation from the coarse grid (the coarse grid is stepped forward in time prior to

advancement of any child grid of that parent).

7.4 Steps to Generate a Nest Grid

Only the concurrent 1-way nest option or the 2-way nest option are considered in this section.

The 1-way nest option (using two consecutive ARW simulations, see Fig. 7.1) is functionally

similar to two separate, single-grid simulations and does not fit the following description. For a

multiple grid simulation within a single model run, there are some additional infrastructure steps

47
U U U U U

U U U U U

U U U U U

U U U U U

V V V V

V V V V

V V V V

V V V V

V V V V

q q q q

q q q q

q q q q

q q q q

U U U U U

U U

U U

U U U U U

V V V V

V V

V V

V V

V V V V

q q q q

q q

q q

q q q q

Figure 7.4: Similar to Fig. 7.3, but with a 2:1 grid-distance ratio.

that are required (briefly described in Fig. 7.5). While the following text details a simulation

with a single coarse-grid and a single fine-grid, this implies no lack of generality when handling

multiple grid levels or multiple grids on the same level.

Nest Instantiation

The fine grid is instantiated as a child of a parent grid at the requested start time. This

initialization is within the integration step for the parent grid, so no child grid can begin if

the parent is not active. To fill in the correct meteorological fields, an initialization routine is

called to horizontally interpolate the coarse-grid data to the fine grid locations using a monotone

interpolation scheme (described in Smolarkiewicz and Grell, 1990) for most fields (i.e., the same

scheme employed for generating the fine grid lateral boundary conditions) and a simple linear

interpolation or averaging scheme for masked or categorical fields. For fields that are masked

with the land/sea background (such as land only fields (e.g., snow), or water only fields (e.g.,

sea ice)), the interpolator needs to know what field defines the template for the masking (such

as the land use category). Part of the automatic code generation handles calling each field with

its associated interpolator.

48
Integrate Parent Grid One Time Step

If Nest Grid Start Time

(1) Horizontally Interpolate Parent to Child Grid

(2) Optionally Input High-Resolution Child Data

(3) Compute Child Reference State

(4) Feedback Child Initial Data to Parent Grid

(5) Re-Compute Parent Reference State

End If Nest Grid Start Time

Solve Time Step for Parent Grid (see Fig. 3.1)

While Existing Nest Grids to Integrate

(1) Lateral Forcing from Parent Grid to Child

(2) Integrate Child Grid to Current Time of Parent Grid

(3) Feedback Child Grid Information to Parent Grid

End While Existing Nest Grids to Integrate

End Grid Integrate

Figure 7.5: Nest grid integration sequence.

Fine Grid Input

After the horizontal interpolation is completed, a few orographic-based variables are saved so

that they may be used to blend the lateral boundaries along the coarse-grid/fine-grid interface.

The terrain elevation, μd

, and the reference geopotential (φ) are stored for later use. The fields

selected as input from the fine grid input file (for the concurrent 1-way and 2-way forecast

methods shown in Fig. 7.1) are ingested, and they overwrite the arrays that were horizontally

interpolated from the coarse grid. No quality control for data consistency is performed for the

fine grid input. All such masked checks are completed by the ARW real-data pre-processor real.

Interface Blended Orography

When the fine grid data has been input, the previously-saved orographic-based fields are blended

across the four outer rows and columns of the fine grid. The blending is a simple linear weighting

between the interpolated coarse-grid values (the saved data) and the fine grid values from the

input file. The weighting scheme is given as:

• row/column 1: 100% interpolated coarse grid, 0% fine grid,

• row/column 2: 75% interpolated coarse grid, 25% fine grid,

• row/column 3: 50% interpolated coarse grid, 50% fine grid,

• row/column 4: 25% interpolated coarse grid, 75% fine grid, and

• row/column 5: 0% interpolated coarse grid, 100% fine grid,

where the row or column nearest the outer edge takes precedence in ambiguous corner zones.

The blended arrays are required to compute the reference state for the fine grid. The first row

and column (100% interpolated from the coarse grid) ensures that the reference state for the

49
coarse grid and fine grid is consistent along the fine grid boundary interface. The blending along

the inner rows and columns ramps the coarse grid reference state to the fine grid reference state

for a smooth transition between the grids.

Feedback

So that the coarse grid and the fine grid are consistent at coincident points, the fine grid values

are fed back to the coarse grid. There are two available options for feedback: either the mean of

all fine grid cells contained within each coarse grid cell is fed back (or cell faces in the case of the

horizontal momentum fields), or a single-point feedback is selected for the masked or categorical

fields.

Subsequent to the feedback step, the coarse grid may be optionally smoothed in the area of

the fine grid. Two smoothers are available: a 5-point 1-2-1 smoother and a smoother-desmoother

with a similar stencil size. Both the feedback and the smoothers are run one row and column

in from the interface row and column of the coarse grid that provides the lateral boundary

conditions to the fine grid.

Reference State

The initial feedback when the nest is instantiated ensures that the coarse grid is consistent

with the fine grid, particularly with regards to elevation and the reference state fields inside the

blended region, and for such terrestrial features as coasts, lakes, and islands. The adjustment of

the elevation in the coarse grid forces a base state recalculation. The fine-grid needs an initial

base state calculation, and after the terrain feedback, the coarse grid is also in need of a base

state recalculation.

Note that with the horizontal interpolation of the coarse grid to the fine grid and the feedback

of the fine grid to the coarse grid, the coarse grid base state is recomputed even without a separate

fine-grid initial data file, since the coarse grid topography is adjusted.

With the completed base state computations, the routines return back to the integration

step for the coarse and fine grids. The fine grid data is now properly initialized for integration

and can be advanced forward a time step.

Integration

The integration by grid is recursive. At the end of each grid’s time step, a check is made to

determine if a child grid exists for that parent and if the current time is bracketed by the child’s

start/end time. This is shown in Fig. 7.5. The integration process for the nest (step 2 under

the while loop) is recursively calling the top step in the overall sequence as a coarse grid itself.

This is a “depth first” traversal of the tree of grids. If a child grid does exist, that child grid is

integrated up through the current time of the parent grid.

50
Chapter 8

Physics

This chapter outlines the physics options available in the ARW. The WRF physics options

fall into several categories, each containing several options. The physics categories are (1)

microphysics, (2) cumulus parameterization, (3) planetary boundary layer (PBL), (4) landsurface model, and (5) radiation. Diffusion, which may also be considered part of the physics,

is described in Chapter 4.

The physics section is insulated from the rest of the dynamics solver by the use of physics

drivers. These are between solver-dependent routines: a pre-physics preparation and postphysics modifications of the tendencies. The physics preparation involves filling arrays with

physics-required variables that include the temperature, pressure, heights, layer thicknesses,

and other state variables in MKS units at half-level grid points and on full levels. The velocities

are also de-staggered so that the physics part is independent of the dynamical solver’s velocity

staggering. Physics packages compute tendencies for the velocity components (un-staggered),

potential temperature, and moisture fields. The solver-dependent post-physics step will restagger these tendencies as necessary, couple tendencies with coordinate metrics, and convert to

variables or units appropriate to the dynamics solver.

In the first Runge-Kutta step, prior to the acoustic steps (see Fig. 3.1, step(1)), tendencies

are computed for radiation, surface, PBL, and cumulus physics. These tendencies are then held

fixed through the Runge-Kutta steps. Microphysics is computed after the last Runge-Kutta

step (see Fig. 3.1, step (9)) in order to maintain proper saturation conditions at the end of the

time-step.

The initialization of the physics is called prior to the first model step. This initialization

may include reading in data files for physics tables or calculating look-up tables of functions.

Each physics module includes an initialization routine for this purpose. Often physics packages

will have many of their own constants that should also be included in their own module, while

common physical constants are passed in from the physics drivers.

8.1 Microphysics

Microphysics includes explicitly resolved water vapor, cloud, and precipitation processes. The

model is general enough to accommodate any number of mass mixing-ratio variables, and other

moments such as number concentrations. Four-dimensional arrays with three spatial indices and

one species index are use to carry such scalars. Memory, i.e., the size of the fourth dimension

51
Table 8.1: Microphysics Options

Scheme Number of Ice-Phase Mixed-Phase

Variables Processes Processes

Kessler 3 N N

Purdue Lin 6 Y Y

WSM3 3 Y N

WSM5 5 Y N

WSM6 6 Y Y

Eta GCP 2 Y Y

Thompson 7 Y Y

in these arrays, is allocated depending on the needs of the scheme chosen, and advection of the

species also applies to all those required by the microphysics option. In the current version of

the ARW, microphysics is carried out at the end of the time-step as an adjustment process, and

so does not provide tendencies. The rationale for this is that condensation adjustment should

be at the end of the time-step to guarantee that the final saturation balance is accurate for the

updated temperature and moisture. However, it is also important to have the latent heating

forcing for potential temperature during the dynamical sub-steps, and this is done by saving the

microphysical heating as an approximation for the next time-step as described in Section 3.1.4.

Currently, the sedimentation process is accounted for in the microphysics, and a smaller

time step is allowed to calculate the vertical flux of precipitation to prevent instability. The

saturation adjustment is also included inside the microphysics. In the future, however, it might

be separated into an individual subroutine to enable the remaining microphysics to be called

less frequently than the model’s advection step for efficiency.

Table 8.1 shows a summary of the options indicating the number of moisture variables, and

whether ice-phase and mixed-phase processes are included. Mixed-phase processes are those

that result from the interaction of ice and water particles, such as riming that produces graupel

or hail. As a general rule, for grid sizes less than 10 km, where updrafts may be resolved,

mixed-phase schemes should be used, particularly in convective or icing situations. For coarser

grids the added expense of these schemes is not worth it because riming is not likely to be well

resolved.

8.1.1 Kessler scheme

This scheme (Kessler, 1969), which was taken from the COMMAS model (Wicker and Wilhelmson, 1995), is a simple warm cloud scheme that includes water vapor, cloud water, and rain. The

microphysical processes included are: the production, fall, and evaporation of rain; the accretion

and autoconversion of cloud water; and the production of cloud water from condensation.

52
8.1.2 Purdue Lin scheme

Six classes of hydrometeors are included: water vapor, cloud water, rain, cloud ice, snow, and

graupel. All parameterization production terms are based on Lin et al. (1983) and Rutledge

and Hobbs (1984) with some modifications, including saturation adjustment following Tao et al.

(1989) and ice sedimentation. This is a relatively sophisticated microphysics scheme in WRF,

and it is more suitable for use in research studies. The scheme is taken from the Purdue cloud

model, and the details can be found in Chen and Sun (2002).

8.1.3 WRF Single-Moment 3-class (WSM3) scheme

This scheme follows Hong et al. (2004) including ice sedimentation and other new ice-phase

parameterizations revised from the older NCEP3 scheme (Hong et al., 1998) that was in WRF

Version 1. A major difference from other schemes is that a diagnostic relation is used for

ice number concentration that is based on ice mass content rather than temperature. Three

categories of hydrometers are included: vapor, cloud water/ice, and rain/snow. As with Dudhia

(1989), this is a so-called simple-ice scheme wherein the cloud ice and cloud water are counted

as the same category. They are distinguished by temperature: namely, cloud ice can only exist

when the temperature is less than or equal to the freezing point; otherwise, cloud water can

exist. The same condition is applied to rain and snow. Though the ice phase is included, it is

considered efficient enough for using in operational models.

8.1.4 WSM5 scheme

This scheme is similar to the WSM3 simple ice scheme. However, vapor, rain, snow, cloud ice,

and cloud water are held in five different arrays. Thus, it allows supercooled water to exist, and

a gradual melting of snow as it falls below the melting layer. Details can be found in Hong et

al. (2004). It replaces WRF Version 1’s NCEP5 scheme (Hong et al., 1998).

8.1.5 WSM6 scheme

The six-class scheme extends the WSM5 scheme to include graupel and its associated processes.

Many of these processes are parameterized similarly to Lin et al. (1983), but there are differences

for the accretion calculation and in some other parameters. The freezing/melting processes are

computed during the fall-term sub-steps to increase accuracy in the vertical heating profile of

these processes. The order of the processes is also optimized to decrease the sensitivity of the

scheme to the time step of the model. As with WSM3 and WSM5, saturation adjustment follows

Dudhia (1989) and Hong et al. (1998) in separately treating ice and water saturation processes,

rather than a combined saturation such as the Purdue Lin (above) and Goddard (Tao et al.,

1989) schemes.

8.1.6 Eta Grid-scale Cloud and Precipitation (2001) scheme

This is also known as EGCP01 or the Eta Ferrier scheme. The scheme predicts changes in

water vapor and condensate in the forms of cloud water, rain, cloud ice, and precipitation ice

(snow/graupel/sleet). The individual hydrometeor fields are combined into total condensate,

53
and it is the water vapor and total condensate that are advected in the model. Local storage arrays retain first-guess information that extract contributions of cloud water, rain, cloud ice, and

precipitation ice of variable density in the form of snow, graupel, or sleet. The density of precipitation ice is estimated from a local array that stores information on the total growth of ice by

vapor deposition and accretion of liquid water. Sedimentation is treated by partitioning the timeaveraged flux of precipitation into a grid box between local storage in the box and fall out through

the bottom of the box. This approach, together with modifications in the treatment of rapid

microphysical processes, permits large time steps to be used with stable results. The mean size

of precipitation ice is assumed to be a function of temperature following the observational results

of Ryan (1996). Mixed-phase processes are now considered at temperatures warmer than -30◦C

(previously -10◦C), whereas ice saturation is assumed for cloudy conditions at colder temperatures. Further description of the scheme can be found in Sec. 3.1 of the November 2001 Technical

Procedures Bulletin (TPB) at http://www.emc.ncep.noaa.gov/mmb/mmbpll/eta12tpb/ and on

the COMET page at http://meted.ucar.edu/nwp/pcu2/etapcp1.htm.

8.1.7 Thompson et al. scheme

The Thompson et al. (2004) microphysical parameterization scheme includes improvements to

the earlier bulk scheme of Reisner et al. (1998) and has been extensively tested and compared

with both idealized case studies and documented real case studies of mid-latitude wintertime

observations. The scheme includes six classes of moisture species plus number concentration

for ice as prognostic variables. The scheme was designed to improve the prediction of freezing

drizzle events for aircraft safety warnings. Generally microphysical parameterizations have had

problems of overpredicting the amount of snow and graupel fields and under predicting the ice

in outflow regions and often not accurately predicting freezing drizzle. Key improvements are

the following:

• Primary ice nucleation as in Cooper (1986), replaces the Fletcher (1962) curve.

• Auto-conversion as in Walko et al. (1995), replaces the Kessler (1969) scheme.

• A generalized gamma distribution for graupel replaces the exponential distribution.

• The associated intercept parameter depends on mixing ratio instead of remaining constant.

• Riming growth of snow must exceed depositional growth of snow by a factor of 3 before

rimmed snow transfers into the graupel category.

• The intercept parameter of the snow size distribution depends on temperature.

• The intercept parameter for the rain size distribution depends on rain mixing ratio, thereby

simulating the fall velocity of drizzle drops as well as raindrops.

8.2 Cumulus parameterization

These schemes are responsible for the sub-grid-scale effects of convective and/or shallow clouds.

The schemes are intended to represent vertical fluxes due to unresolved updrafts and downdrafts

and compensating motion outside the clouds. They operate only on individual columns where

the scheme is triggered and provide vertical heating and moistening profiles. Some schemes

additionally provide cloud and precipitation field tendencies in the column, and future schemes

54
Table 8.2: Cumulus Parameterization Options

Scheme Cloud Type of scheme Closure

Detrainment

Kain-Fritsch Y Mass flux CAPE removal

Betts-Miller-Janjic N Adjustment Sounding adjustment

Grell-Devenyi Y Mass flux Various

may provide momentum tendencies due to convective transport of momentum. The schemes all

provide the convective component of surface rainfall.

Cumulus parameterizations are theoretically only valid for coarser grid sizes, (e.g., greater

than 10 km), where they are necessary to properly release latent heat on a realistic time scale

in the convective columns. While the assumptions about the convective eddies being entirely

sub-grid-scale break down for finer grid sizes, sometimes these schemes have been found to be

helpful in triggering convection in 5–10 km grid applications. Generally, they should not be

used when the model can resolve the convective eddies itself (e.g., ≤ 5 km grid).

Table 8.2 summarizes the basic characteristics of the available cumulus parameterization

options in the ARW.

8.2.1 Kain-Fritsch

The modified version of the Kain-Fritsch scheme (KF-Eta) is based on Kain and Fritsch (1990)

and Kain and Fritsch (1993), but has been modified based on testing within the Eta model. As

with the original KF scheme, it utilizes a simple cloud model with moist updrafts and downdrafts,

including the effects of detrainment, entrainment, and relatively simple microphysics. It differs

from the original KF scheme in the following ways:

• A minimum entrainment rate is imposed to suppress widespread convection in marginally

unstable, relatively dry environments.

• Shallow (non precipitating) convection is allowed for any updraft that does not reach

minimum cloud depth for precipitating clouds; this minimum depth varies as a function

of cloud-base temperature.

• The entrainment rate is allowed to vary as a function of low-level convergence.

• Downdraft changes:

– Source layer is the entire 150 – 200 mb deep layer just above cloud base.

– Mass flux is specified as a fraction of updraft mass flux at cloud base. Fraction is

a function of source layer RH rather than wind shear or other parameters, i.e., old

precipitation efficiency relationship not used.

– Detrainment is specified to occur in updraft source layer and below.

55
8.2.2 Betts-Miller-Janjic

The Betts-Miller-Janjic (BMJ) scheme (Janjic, 1994, 2000) was derived from the Betts-Miller

(BM) convective adjustment scheme (Betts, 1986; Betts and Miller, 1986). However, the BMJ

scheme differs from the Betts-Miller scheme in several important aspects. The deep convection profiles and the relaxation time are variable and depend on the cloud efficiency, a nondimensional parameter that characterizes the convective regime (Janjic, 1994). The cloud efficiency depends on the entropy change, precipitation, and mean temperature of the cloud. The

shallow convection moisture profile is derived from the requirement that the entropy change

be small and nonnegative (Janjic, 1994). The BMJ scheme has been optimized over years of

operational application at NCEP, so that, in addition to the described conceptual differences,

many details and/or parameter values differ from those recommended in Betts (1986) and Betts

and Miller (1986). Recently, attempts have been made to refine the scheme for higher horizontal

resolutions, primarily through modifications of the triggering mechanism. In particular:

• A floor value for the entropy change in the cloud is set up below which the deep convection

is not triggered;

• In searching for the cloud top, the ascending particle mixes with the environment; and

• The work of the buoyancy force on the ascending particle is required to exceed a prescribed

positive threshold.

8.2.3 Grell-Devenyi ensemble

Grell and Devenyi (2002) introduced an ensemble cumulus scheme in which effectively multiple

cumulus schemes and variants are run within each grid box and then the results are averaged

to give the feedback to the model. In principle, the averaging can be weighted to optimize the

scheme, but the default is an equal weight. The schemes are all mass-flux type schemes, but with

differing updraft and downdraft entrainment and detrainment parameters, and precipitation

efficiencies. These differences in static control are combined with differences in dynamic control,

which is the method of determining cloud mass flux. The dynamic control closures are based on

convective available potential energy (CAPE or cloud work function), low-level vertical velocity,

or moisture convergence. Those based on CAPE either balance the rate of change of CAPE or

relax the CAPE to a climatological value, or remove the CAPE in a convective time scale. The

moisture convergence closure balances the cloud rainfall to the integrated vertical advection

of moisture. Another control is the trigger, where the maximum cap strength that permits

convection can be varied. These controls typically provide ensembles of 144 members.

8.3 Surface Layer

The surface layer schemes calculate friction velocities and exchange coefficients that enable

the calculation of surface heat and moisture fluxes by the land-surface models and surface

stress in the planetary boundary layer scheme. Over water surfaces, the surface fluxes and

surface diagnostic fields are computed in the surface layer scheme itself. The schemes provide no

tendencies, only the stability-dependent information about the surface layer for the land-surface

and PBL schemes. Currently, each surface layer option is tied to particular boundary-layer

options, but in the future more interchangeability and options may become available.

56
8.3.1 Similarity theory (MM5)

This scheme uses stability functions from Paulson (1970), Dyer and Hicks (1970), and Webb

(1970) to compute surface exchange coefficients for heat, moisture, and momentum. A convective

velocity following Beljaars (1994) is used to enhance surface fluxes of heat and moisture. No

thermal roughness length parameterization is included in the current version of this scheme.

A Charnock relation relates roughness length to friction velocity over water. There are four

stability regimes following Zhang and Anthes (1982). This surface layer scheme must be run in

conjunction with the MRF or YSU PBL schemes.

8.3.2 Similarity theory (Eta)

The Eta surface layer scheme (Janjic, 1996, 2002) is based on similarity theory (Monin and

Obukhov, 1954). The scheme includes parameterizations of a viscous sub-layer. Over water

surfaces, the viscous sub-layer is parameterized explicitly following Janjic (1994). Over land,

the effects of the viscous sub-layer are taken into account through variable roughness height for

temperature and humidity as proposed by Zilitinkevich (1995). The Beljaars (1994) correction

is applied in order to avoid singularities in the case of an unstable surface layer and vanishing

wind speed. The surface fluxes are computed by an iterative method. This surface layer scheme

must be run in conjunction with the Eta (Mellor-Yamada-Janjic) PBL scheme, and is therefore

sometimes referred to as the MYJ surface scheme.

8.4 Land-Surface Model

The land-surface models (LSMs) use atmospheric information from the surface layer scheme,

radiative forcing from the radiation scheme, and precipitation forcing from the microphysics

and convective schemes, together with internal information on the land’s state variables and

land-surface properties, to provide heat and moisture fluxes over land points and sea-ice points.

These fluxes provide a lower boundary condition for the vertical transport done in the PBL

schemes (or the vertical diffusion scheme in the case where a PBL scheme is not run, such as in

large-eddy mode). [Note that large-eddy mode with interactive surface fluxes is not yet available

in the ARW, but is planned for the near future.] The land-surface models have various degrees of

sophistication in dealing with thermal and moisture fluxes in multiple layers of the soil and also

may handle vegetation, root, and canopy effects and surface snow-cover prediction. The landsurface model provides no tendencies, but does update the land’s state variables which include

the ground (skin) temperature, soil temperature profile, soil moisture profile, snow cover, and

possibly canopy properties. There is no horizontal interaction between neighboring points in the

LSM, so it can be regarded as a one-dimensional column model for each WRF land grid-point,

and many LSMs can be run in a stand-alone mode. Table 8.3 summarizes the basic features of

the land-surface treatments in ARW.

8.4.1 5-layer thermal diffusion

This simple LSM is based on the MM5 5-layer soil temperature model. Layers are 1, 2, 4, 8, and

16 cm thick. Below these layers, the temperature is fixed at a deep-layer average. The energy

57
Table 8.3: Land Surface Options

Scheme Vegetation Soil Snow

Processes Variables (Layers) Scheme

5-layer N Temperature (5) none

Noah Y Temperature, Water+Ice, Water (4) 1-layer, fractional

RUC Y Temperature, Ice, Water + Ice (6) multi-layer

budget includes radiation, sensible, and latent heat flux. It also allows for a snow-cover flag, but

the snow cover is fixed in time. Soil moisture is also fixed with a landuse- and season-dependent

constant value, and there are no explicit vegetation effects.

8.4.2 Noah LSM

The Noah LSM is the successor to the OSU LSM described by Chen and Dudhia (2001). The

scheme was developed jointly by NCAR and NCEP, and is a unified code for research and

operational purposes, being almost identical to the code used in the NCEP North American

Mesoscale Model (NAM). This has the benefit of being consistent with the time-dependent soil

fields provided in the analysis datasets. This is a 4-layer soil temperature and moisture model

with canopy moisture and snow cover prediction. It includes root zone, evapotranspiration, soil

drainage, and runoff, taking into account vegetation categories, monthly vegetation fraction,

and soil texture. The scheme provides sensible and latent heat fluxes to the boundary-layer

scheme. The Noah LSM additionally predicts soil ice, and fractional snow cover effects, has an

improved urban treatment, and considers surface emissivity properties, which are all new since

the OSU scheme.

8.4.3 Rapid Update Cycle (RUC) Model LSM

This is a LSM with 6 sub-soil layers and up to two snow layers that is used operationally in the

RUC model (Smirnova et al., 1997, 2000). The model considers frozen soil processes, patchy

snow, with snow temperature and density variation, vegetation effects, and canopy water.

8.5 Planetary Boundary Layer

The planetary boundary layer (PBL) is responsible for vertical sub-grid-scale fluxes due to eddy

transports in the whole atmospheric column, not just the boundary layer. Thus, when a PBL

scheme is activated, explicit vertical diffusion is de-activated with the assumption that the PBL

scheme will handle this process. The most appropriate horizontal diffusion choices (Section

4.1.3) are those based on horizontal deformation or constant Kh values where horizontal and

vertical mixing are treated independently. The surface fluxes are provided by the surface layer

58
Table 8.4: Planetary Boundary Layer Options

Scheme Unstable PBL Entrainment PBL Top

Mixing treatment

MRF K profile + countergradient term part of PBL mixing from critical bulk Ri

YSU K profile + countergradient term explicit term from buoyancy profile

MYJ K from prognostic TKE part of PBL mixing from TKE

and land-surface schemes. The PBL schemes determine the flux profiles within the well-mixed

boundary layer and the stable layer, and thus provide atmospheric tendencies of temperature,

moisture (including clouds), and horizontal momentum in the entire atmospheric column. Most

PBL schemes consider dry mixing, but can also include saturation effects in the vertical stability

that determines the mixing. The schemes are one-dimensional, and assume that there is a clear

scale separation between sub-grid eddies and resolved eddies. This assumption will become less

clear at grid sizes below a few hundred meters, where boundary layer eddies may start to be

resolved, and in these situations the scheme should be replaced by a fully three-dimensional

local sub-grid turbulence scheme such as the TKE diffusion scheme (Section 4.1.4.) Table 8.4

summarizes the basic features of the PBL schemes in ARW.

8.5.1 Medium Range Forecast Model (MRF) PBL

The scheme is described by Hong and Pan (1996). This PBL scheme employs a so-called

counter-gradient flux for heat and moisture in unstable conditions. It uses enhanced vertical

flux coefficients in the PBL, and the PBL height is determined from a critical bulk Richardson

number. It handles vertical diffusion with an implicit local scheme, and it is based on local Ri

in the free atmosphere.

8.5.2 Yonsei University (YSU) PBL

The Yonsei University PBL is the next generation of the MRF PBL, also using the countergradient terms to represent fluxes due to non-local gradients. This adds to the MRF PBL

an explicit treatment of the entrainment layer at the PBL top. The entrainment is made

proportional to the surface buoyancy flux in line with results from studies with large-eddy

models. The PBL top is defined using a critical bulk Richardson number of zero (compared to

0.5 in the MRF PBL), so is effectively only dependent on the buoyancy profile which, in general,

lowers the calculated PBL top compared to MRF.

8.5.3 Mellor-Yamada-Janjic (MYJ) PBL

This parameterization of turbulence in the PBL and in the free atmosphere (Janjic, 1990, 1996,

2002) represents a nonsingular implementation of the Mellor-Yamada Level 2.5 turbulence clo59
Table 8.5: Radiation Options

Scheme Longwave/ Spectral CO2, O3, clouds

Shortwave Bands

RRTM LW 16 CO2, O3, clouds

GFDL LW LW 14 CO2, O3, clouds

GFDL SW SW 12 CO2, O3, clouds

MM5 SW SW 1 clouds

Goddard SW 11 CO2, O3, clouds

sure model (Mellor and Yamada, 1982) through the full range of atmospheric turbulent regimes.

In this implementation, an upper limit is imposed on the master length scale. This upper limit

depends on the TKE as well as the buoyancy and shear of the driving flow. In the unstable

range, the functional form of the upper limit is derived from the requirement that the TKE production be nonsingular in the case of growing turbulence. In the stable range, the upper limit

is derived from the requirement that the ratio of the variance of the vertical velocity deviation

and TKE cannot be smaller than that corresponding to the regime of vanishing turbulence. The

TKE production/dissipation differential equation is solved iteratively. The empirical constants

have been revised as well (Janjic, 1996, 2002).

8.6 Atmospheric Radiation

The radiation schemes provide atmospheric heating due to radiative flux divergence and surface

downward longwave and shortwave radiation for the ground heat budget. Longwave radiation

includes infrared or thermal radiation absorbed and emitted by gases and surfaces. Upward

longwave radiative flux from the ground is determined by the surface emissivity that in turn

depends upon land-use type, as well as the ground (skin) temperature. Shortwave radiation

includes visible and surrounding wavelengths that make up the solar spectrum. Hence, the only

source is the Sun, but processes include absorption, reflection, and scattering in the atmosphere

and at surfaces. For shortwave radiation, the upward flux is the reflection due to surface albedo.

Within the atmosphere the radiation responds to model-predicted cloud and water vapor distributions, as well as specified carbon dioxide, ozone, and (optionally) trace gas concentrations.

All the radiation schemes in WRF currently are column (one-dimensional) schemes, so each column is treated independently, and the fluxes correspond to those in infinite horizontally uniform

planes, which is a good approximation if the vertical thickness of the model layers is much less

than the horizontal grid length. This assumption would become less accurate at high horizontal

resolution. Table 8.5 summarizes the basic features of the radiation schemes in the ARW.

60
8.6.1 Rapid Radiative Transfer Model (RRTM) Longwave

This RRTM, which is taken from MM5, is based on Mlawer et al. (1997) and is a spectral-band

scheme using the correlated-k method. It uses pre-set tables to accurately represent longwave

processes due to water vapor, ozone, CO2, and trace gases (if present), as well as accounting for

cloud optical depth.

8.6.2 Eta Geophysical Fluid Dynamics Laboratory (GFDL) Longwave

This longwave radiation scheme is from GFDL. It follows the simplified exchange method of

Fels and Schwarzkopf (1975) and Schwarzkopf and Fels (1991), with calculation over spectral

bands associated with carbon dioxide, water vapor, and ozone. Included are Schwarzkopf and

Fels (1985) transmission coefficients for carbon dioxide, a Roberts et al. (1976) water vapor

continuum, and the effects of water vapor-carbon dioxide overlap and of a Voigt line-shape

correction. The Rodgers (1968) formulation is adopted for ozone absorption. Clouds are randomly overlapped. This scheme is implemented to conduct comparisons with the operational

Eta model.

8.6.3 Eta Geophysical Fluid Dynamics Laboratory (GFDL) Shortwave

This shortwave radiation is a GFDL version of the Lacis and Hansen (1974) parameterization.

Effects of atmospheric water vapor, ozone (both from Lacis and Hansen, 1974), and carbon

dioxide (Sasamori et al., 1972) are employed. Clouds are randomly overlapped. Shortwave

calculations are made using a daylight-mean cosine solar zenith angle over the time interval

(given by the radiation call frequency).

8.6.4 MM5 (Dudhia) Shortwave

This scheme is base on Dudhia (1989) and is taken from MM5. It has a simple downward

integration of solar flux, accounting for clear-air scattering, water vapor absorption (Lacis and

Hansen, 1974), and cloud albedo and absorption. It uses look-up tables for clouds from Stephens

(1978).

8.6.5 Goddard Shortwave

This scheme is based on Chou and Suarez (1994). It has a total of 11 spectral bands and

considers diffuse and direct solar radiation components in a two-stream approach that accounts

for scattered and reflected components. Ozone is considered with several climatological profiles

available.

61
Table 8.6: Physics Interactions. Columns correspond to model physical processes: radiation

(Rad), microphysics (MP), cumulus parameterization (CP), planetary boundary layer/vertical

diffusion (PBL), and surface physics (Sfc). Rows corresponds to model variables where i and o

indicate whether a variable is input or output (updated) by a physical process.

Rad MP CP PBL Sfc

Atmospheric Momentum i io

State or Pot. Temp. io io io io

Tendencies Water Vapor i io io io

Cloud i io o io

Precip i io o

Surface Longwave Up i o

Fluxes Longwave Down o i

Shortwave Up i o

Shortwave Down o i

Sfc Convective Rain o i

Sfc Resolved Rain o i

Heat Flux i o

Moisture Flux i o

Surface Stress i o

62
8.7 Physics Interactions

While the model physics parameterizations are categorized in a modular way, it should be

noted that there are many interactions between them via the model state variables (potential

temperature, moisture, wind, etc.) and their tendencies, and via the surface fluxes. Table 8.6

summarizes how the various physics processes interact in the model. In the table, i indicates

that the state variable or flux is required input for the physics scheme, and o indicates that the

tendency or flux is a probable output of the scheme. It can be seen that all the physical schemes

interact in some way with the surface physics (land-surface models, and, potentially, coupled

ocean models). The surface physics, while not explicitly producing tendencies of atmospheric

state variables, is responsible for updating the land-state variables.

Note also that, as mentioned, the microphysics does not output tendencies, but updates the

atmospheric state at the end of the model time-step. However, the rest of the o’s in the upper

half of the table are representative of the physical tendencies of these variables in the model.

The radiation, cumulus parameterization, and boundary-layer schemes all output tendencies,

but the tendencies are not added until later in the solver, so from this perspective the order

of call is not important. Moreover, these physics schemes do not have to be called at the

same frequency as each other or the model time step. When lower frequencies are used, their

tendencies are kept constant between calls. This is typically done for the radiation schemes,

which are too expensive to call every time, and for the cumulus schemes, for which it is also not

necessary. However, the surface/boundary-layer schemes are normally called every step in the

ARW because this is likely to give the best results.

The radiation is called first because of the required radiative fluxes that are input to the

land-surface scheme. The land-surface also requires rainfall from the microphysics and cumulus

schemes, but that is from the previous time-step. The boundary-layer scheme is necessarily after

the land-surface scheme because it requires the heat and moisture fluxes.

63
64
Chapter 9

Variational Data Assimilation

An introduction to the basic ideas of variational data assimilation and the WRF-Var system is

given in this chapter, followed by a brief overview of recent major improvements to WRF-Var.

This overview supplements the original description of the three-dimensional variational (3D-Var)

algorithm found in Barker et al. (2004). One of the most important additions to WRF-Var is a

new utility gen be, used to calculate background error covariances for a user’s own application;

it is discussed in the latter half of this chapter. The WRF-Var system is evolving rapidly, and a

future technical note will accompany the general release of the 4D-Var component of WRF-Var.

That technical note will include an extensive description of the entire WRF-Var system.

9.1 Introduction

The basic goal of any variational data assimilation system is to produce an optimal estimate

of the true atmospheric state at analysis time through iterative solution of a prescribed costfunction (Ide et al., 1997):

J(x) = Jb(x) + Jo(x) = 1

2

(x − x

b

)

T B−1

(x − x

b

) + 1

2

(y − y

o

)

T

(E + F)

−1

(y − y

o

). (9.1)

The variational problem can be summarized as the iterative minimization of (9.1) to find the

analysis state x that minimizes J(x). This solution represents the a posteriori maximum likelihood (minimum variance) estimate of the true state of the atmosphere given the two sources

of a priori data: the first guess (or background) x

b and observations yo

(Lorenc, 1986). The

fit to individual data points is weighted by estimates of their errors: B, E, and F are the background, observation (instrumental), and representivity error covariance matrices, respectively.

Representivity error is an estimate of inaccuracies introduced in the observation operator H

used to transform the gridded analysis x to observation space y = H(x) for comparison against

observations. This error will be resolution dependent and may also include a contribution from

approximations (e.g., linearizations) in H.

As described in Barker et al. (2004), the particular variational data assimilation algorithm

adopted in WRF-Var is a model-space, incremental formulation of the variational problem. In

this approach, observations, previous forecasts, their errors, and physical laws are combined to

produce analysis increments x

a

0

, which are added to the first guess x

b

to provide an updated

analysis.

65
Figure 9.1: Sketch showing the relationship between datasets (circles), and algorithms (rectangles) of the ARW system.

Figure 9.1 illustrates the relationship between WRF-Var, the various datasets, and the other

components of a typical NWP system (here ARW). The WRF-Var assimilation proceeds as

described in Barker et al. (2004). A number of recent upgrades to the WRF-Var algorithm will

be described in Section 9.2.

The three inputs to WRF-Var are:

a) First guess x

b— In cold-start mode, this is typically a forecast/analysis from another

model interpolated to the ARW grid (and variables) via the WRF SI and real programs. In

cycling mode, the first guess is a short-range (typically 1–6 hour) ARW forecast.

b) Observations y

o— In the current version of WRF-Var, observations may be supplied

either in a text (MM5 3D-Var) format or BUFR format (but not a combination of the two).

An observation preprocessor (3DVAR OBSPROC) is supplied with the code release to perform

basic quality control, assign observation errors, and reformat observations from the MM5 little r

text format into 3D-Var’s own text format. Details can be found in Barker et al. (2003, 2004).

c) Background error covariances B— used to define the spatial and multivariate response of

the analysis to an observation. In variational systems, these covariances are typically calculated

off-line, and significant tuning is required to optimize performance for a particular application

(e.g., Ingleby (2001); Wu et al. (2002)). The amount of work required to do this satisfactorily

is significant, and should not be underestimated. In order to assist the user, WRF developers

66
supply the following: i) a default set of statistics used for the initial set up of a domain; ii) a

utility gen be (described in Section 9.3) to process ensembles of forecasts into the appropriate

control variable space; and iii) diagnostic routines to assess the accuracy of observation and

background error statistics. These routines include both innovation vector-based approaches

(Hollingsworth and Lonnberg, 1986) and variational tuning approaches (Desroziers and Ivanov,

2001).

Following assimilation of all data, an analysis x

a

is produced that must be merged with the

existing lateral boundary conditions x

lbc (described in Barker et al. (2003)). Note: In cycling

mode, only the wrfbdy lateral boundary condition files (x

lbc) output of SI/real are used, and

not the wrfinput initial condition files (x

b

). In cold-start mode, both are required.

9.2 Improvements to the WRF-Var Algorithm

9.2.1 Improved vertical interpolation

The original WRF 3D-Var system described in Barker et al. (2004) used height interpolation

for all observation operators. If an observation is reported as a function of pressure, then height

is approximated using the hydrostatic relation. This step introduces an unnecessary source

of error. The new WRF-Var system performs vertical interpolation in terms of the original

observed coordinate, i.e., pressure or height.

9.2.2 Improved minimization and “outer loop”

The default WRF-Var cost function minimization uses a modified version of the limited memory

Quasi-Newton Method (QNM). Recently, an alternative Conjugate Gradient Method (CGM) has

been implemented. Unlike the QNM technique, the CGM method restricts 3D-Var’s inner loop

to be completely linear. This limitation is dealt with through the inclusion of an outer loop

in WRF-Var, the purpose of which is to iterate towards nonlinear solutions (e.g., observation

operators, balance constraints, and the forecast itself in 4D-Var) using the WRF-Var analysis

from the previous iteration as new background. The outer loop is also used as a form of

variational quality control as follows: observations are rejected if their O-B values are outside

a prescribed range (typically several times the observation error standard deviation). This

errormax test implicitly assumes the rejected large O-B values are due to a bad observation

(O) rather than poor background (B). However, if it is the background B that is incorrect then

the system will reject the most useful observations available to the assimilation system, i.e.,

those in areas where the first-guess is poor. The outer loop alleviates this effect by allowing

observations rejected in previous iterations to be accepted if their new O-B falls within the

required range in subsequent outer loops. The assimilation of nearby observations in previous

iterations essentially provides a “buddy check” to the observation in question.

9.2.3 Flexible choice of control variables

In practical variational data assimilation schemes, the background error covariance matrix B is

computed not in model space x

0

: u, v, T, q, ps, but in a control variable space v related to model

67
space via the control variable transform U, i.e.,

x

0 = Uv = UpUvUhv. (9.2)

The expansion U = UpUvUh represents the various stages of covariance modeling: horizontal

correlations Uh, vertical covariances Uv, and multivariate covariances Up.

The components of v are chosen so that their error cross-correlations are negligible, thus

permitting the matrix B to be block-diagonalized. The many varying applications (high/low

resolution, polar/tropical, etc.) of WRF-Var require a flexible choice of background error model.

This is achieved via a namelist option “cv options” as defined in Table 9.1.

2 3 4 5

cv options (original MM5) (NCEP) (Global) (Regional)

Analysis

x

0 u0

,v

0

,T

0

,q

0

,ps

0

(i, j, k)

Increment

Change of Up ψ

0

,χ

0

,p

0

u

,q

0 ψ

0

,χ

0

u

,T

0

u

,r

0

,p

0

Variable su

Vertical Uv B = EΛET RF B = EΛET

Covariances

Horizontal Uh RF Spectral RF Correlations

Control

v v(i, j, m) v(i, j, k) v(l, n, m) v(i, j, m)

Variables

Table 9.1: The definitions of the various stages of the control variable transform given by (9.2) for

the unified global/regional WRF-Var system. Indices (i, j, k) refer to grid-point space, index m

to vertical mode, and l, n to global spectral mode. The variables are: u, v: velocity components;

T: temperature; q: mixing ratio; ps: surface pressure; ψ: streamfunction; χ: velocity potential;

r: relative humidity. The subscript u indicates an unbalanced field. The acronym RF stands

for recursive filter.

Table 9.1 indicates that the only difference between global (cv options=4) and WRF regional

(cv options=5) versions of the WRF-Var control variable transform is in the horizontal error

correlations Uh. Note also, the only difference between the old MM5 background error model

(cv options=2) and WRF regional (cv options=5) is in the Up transform. The former imposes a

dynamical balance constraint via an unbalanced pressure control variable (Barker et al., 2004),

whereas in the new regional covariance model, balance is imposed via statistical regression (see

Section 9.3 for details). This choice of control variables is considered more appropriate for the

mass-based ARW solver.

68
9.2.4 First Guess at Appropriate Time (FGAT)

A First Guess at Appropriate Time (FGAT) procedure has been implemented in WRF-Var (Lee

et al., 2004). The FGAT capability results in a more accurate calculation of innovation vectors

(observation minus first guess difference), and hence a better use of observations when their valid

time differs from that of the analysis. FGAT is most effective for the analysis of observations

from asynoptic, moving platforms (e.g., aircraft and satellite data). Surface observations with

high temporal resolution also benefit from the use of FGAT.

9.2.5 Radar Data Assimilation

Numerous modifications have been made in order to assimilate Doppler radar radial velocity

and reflectivity observations. Firstly, vertical velocity increments are included in WRF-Var via

the “Richarson balance equation” that combines the continuity equation, adiabatic thermodynamic equation, and hydrostatic relation. Linear and adjoint codes of Richardson’s equation

have been incorporated into WRF-Var. In order to develop a capability for Doppler reflectivity

assimilation, we use the total water as a control variable, requiring a partitioning of the moisture

and water hydrometeor increments. A warm-rain parameterization is also included, which includes condensation of water vapor into cloud, accretion of cloud by rain, automatic conversion

of cloud to rain, and evaporation of rain to water vapor. Finally, the observation operators for

Doppler radial velocity and reflectivity are included in WRF-Var. Further details and results of

the radial velocity work can be found in Xiao et al. (2005). The radar reflectivity approach will

be described in a future paper.

9.2.6 Unified Regional/Global 3D-Var Assimilation

There are many benefits to having a single data assimilation system for regional and global

applications. The majority of the code is common to both (observation operators, minimization,

background error preconditioning, interpolation, etc.). Technically, the main changes required

to extend the regional application to global are related to the presence of a) the polar singularity,

and b) periodic East-West boundary conditions. Of course, there are also scientific questions

concerning the optimal mix of observations required for global/regional models, and the choice

of control variables and balance constraints. A unified global/regional 3D-Var system should

therefore be flexible to a variety of thinning/quality-control algorithms and also to alternative

formulations of the background error covariance matrix. This flexibility has been a key design

requirement throughout the WRF-Var project.

The major difference between regional and global options in WRF-Var is in the definition

of horizontal background error covariances. In regional mode, recursive filters (Purser et al.,

2003) are used to project observation information away from the observation location. In global

mode, a spectral decomposition is applied. The main benefits of the spectral technique are a)

the isotropic and homogeneous covariances that are implied neatly solve the problems associated

with the pole (the pole is not a special point in spectral space), and b) horizontal correlations

are defined in terms of a single function— the power spectrum (a function of total wavenumber).

However, the isotropy of correlation defined in spectral space is also a weakness— anisotropies

need to be defined in an alternative manner. One solution to this problem is to replace the

69
spectral correlations with grid-point correlations (e.g., in the Gridpoint Statistical Interpolation

scheme under development at NCEP). An alternative technique is to supplement the isotropic

spectral correlations with an anisotropic component derived via grid transformations, additional

control variables or 4D-Var. Research using the latter techniques is underway using the WRFVar system.

The WRF-Var code has been adapted to perform assimilation on a global, regular latitudelongitude grid. The major modifications are as follows.

a) Spectral to grid-point transformations (and their adjoints) have been included in 3D-Var

to represent the horizontal component (Uh) of the background error covariance model.

b) A new global WRF-Var registry was created to accommodate the output of global forecast

models (currently only the Korean Meteorological Administration’s (KMA) global model has

been tested). The final analysis increments are written in binary format and added back to the

global first guess to provide the global analysis in the native model format.

c) Changes have been made to allow for periodic boundary conditions in the East–West

direction.

d) A number of minor changes have been made to treat the polar rows as special points (e.g.,

in the grid-point ψ, χ to u, v wind conversion in the Up transform and the observation operators

for polar winds).

9.3 Background Error Covariances

Forecast (“first guess” or “background”) error covariances are a vital input to variational data

assimilation systems. They influence the analysis fit to observations and also completely define

the analysis response away from observations. The latter impact is particularly important

in data-sparse areas of the globe. Unlike ensemble filter data assimilation techniques (e.g.,

the Ensemble Adjustment Kalman Filter, the Ensemble Transform Kalman Filter), 3/4D-Var

systems do not implicitly evolve forecast error covariances in real-time. Instead, climatologic

statistics are usually estimated offline. The “NMC-method”, in which forecast error covariances

are approximated using forecast difference (e.g., T+48 minus T+24) statistics, is a commonly

used approach (Parrish and Derber, 1992). Experiments at ECMWF (Fisher, 2003) indicate

superior statistics may be obtained using a cycling analysis/forecast ensemble prediction system

based on perturbed observations/physics.

Recent advances permit the use of flow-dependent forecast error covariances in 3D/4D-Var

through, for example, grid transformations (Desroziers, 1997), anisotropic recursive filters (Wu

et al., 2002; Purser et al., 2003), or observation-space formulations of the variational problem

(Daley and Barker, 2001). Flow-dependence may be enhanced in 4D-Var through the use of the

forecast model to provide dynamical consistency to the evolving forecast state during 4D-Var’s

time window (Rabier et al., 1998). Still, the practical effort required to specify and implement

flow-dependent error covariances in 3/4D-Var is significant.

The NMC-method code developed for MM5 3D-Var (Barker et al., 2004) is nearing the end

70
of its useful life. The development of a unified global/regional WRF-Var system, and its application to a variety of models (e.g., ARW, MM5, KMA global model, Taiwan’s Nonhydrostatic

Forecast System [NFS]) has required a new, efficient, portable forecast background error covariance calculation code to be written. There is also a demand for such a capability to be available

and supported for the wider 3/4D-Var research community for application to their own geographic areas of interest (the default statistics supplied with the WRF-Var release are designed

only as a starting point). In this section, the new gen be code developed by NCAR/MMM to

generate forecast error statistics for use with the WRF-Var system is described.

The background error covariance matrix is defined as

B = T ' x0x0T

, (9.3)

where the overbar denotes an average over time and/or geographical area. The true background

error  is not known in reality, but is assumed to be statistically well-represented by a model state

perturbation x

0

. In the standard NMC-method (Parrish and Derber, 1992), the perturbation

x

0

is given by the difference between two forecasts (e.g., 24 hour minus 12 hour) verifying

at the same time. Climatological estimates of background error may then be obtained by

averaging such forecast differences over a period of time (e.g., one month). An alternative

strategy proposed by (Fisher, 2003) makes use of ensemble forecast output, defining the x

0

vectors as ensemble perturbations (ensemble minus ensemble mean). In either approach, the

end result is an ensemble of model perturbation vectors from which estimates of background

error may be derived. The new gen be utility has been designed to work with either forecast

difference, or ensemble-based, perturbations.

As described above, the WRF-Var background error covariances are specified not in model

space x

0

, but in a control variable space v, which is related to the model variables (e.g., wind

components, temperature, humidity, and surface pressure) via the control variable transform

defined in (9.2). Both (9.2) and its adjoint are required in WRF-Var. In contrast, the background error code performs the inverse control variable transform v = U−1

h U−1

v U−1

p x

0

in order

to accumulate statistics for each component of the control vector v.

Using the NMC-method, x

0 = xT2 − xT1 where T2 and T1 are the forecast difference times

(e.g., 48h minus 24h for global, 24h minus 12h for regional). Alternatively, for an ensemble-based

approach, xk

0 = xk − x ̄, where the overbar is an average over ensemble members k = 1, ne. The

total number of binary files produced by this stage is nf ×ne where nf is the number of forecast

times used (e.g., for 30 days with forecast every 12 hours, nf = 60). Using the NMC-method,

ne = 1 (1 forecast difference per time). For ensemble-based statistics, ne is the number of

ensemble members.

The background error covariance generation code gen be is designed to process data from a

variety of regional/global models (e.g., ARW, MM5, KMA global model, NFS, etc.), and process

it in order to provide error covariance statistics for use in variational data assimilation systems.

The initial, model-dependent “stage 0” is illustrated in Fig. 9.2.

Alternative models use different grids, variables, data formats, etc., and so initial converters

are required to transform model output into a set of standard perturbation fields (and metadata),

and to output them in a standard binary format for further, model-independent processing. The

standard grid-point fields are as follows.

• Perturbations: Streamfunction ψ

0

(i, j, k), velocity potential χ

0

(i, j, k), temperature T

0

(i, j, k),

71
Figure 9.2: Sketch of the role of Stage 0 converters in transforming model-specific data (e.g.,

ARW, KMA global model, etc.) to standard perturbation fields and relevant metadata (e.g.,

latitude, height, land/sea, etc.).

relative humidity r

0

(i, j, k), surface pressure p

0

s

(i, j).

• Full-fields: height z(i, j, k), latitude φ(i, j). (These are required for the production of

background error statistics stored in terms of physics variables, rather than tied to a

specific grid. This flexibility is included in gen be through a namelist option to define the

bins over which data is averaged in a variety of ways (e.g., latitude height, grid points).

Land-sea and orographic effects may also be represented in this way.

In general, the stage 0 converters are developed and maintained by those supporting individual models. Only the WRF-to-standard-fields converter gen be stage0 wrf is maintained and

supported by the ARW effort.

9.3.1 Removal of time-mean

In order to calculate covariances between fields, the average value must first be removed. This

is performed in the first stage utility gen be stage1.

9.3.2 Multivariate Covariances: Regression coefficients and unbalanced variables

The WRF-Var system permits a variety of background error covariance models to be employed,

as described in Section 9.2.3 above. The utility gen be is used to provide background error

statistics only for cv options 4 and 5.

72
The second stage of gen be (gen be stage2) provides statistics for the unbalanced fields χu,

Tu, and Psu used as control variables in WRF-Var. The unbalanced control variables are defined

as the difference between full and balanced (or correlated) components of the field. In this stage

of the calculation of background errors, the balanced component of particular fields is modeled

via a regression analysis of the field using specified predictor fields (e.g., streamfunction; see

Wu et al. (2002) for further details). The resulting regression coefficients are output for use

in WRF-Var’s Up transform. Currently, three regression analyses are performed resulting in

three sets of regression coefficients (Note: The perturbation notation has been dropped for the

remainder of this chapter for clarity.):

• Velocity potential/streamfunction regression: χb = cψ;

• Temperature/streamfunction regression: Tb,k1 =

P

k2 Gk1,k2ψk2; and

• Surface pressure/streamfunction regression: psb =

PkWkψk.

Data is read from all nf × ne files and sorted into bins defined via the namelist option

bin type. Regression coefficients G(k1, k2) and W(k) are computed individually for each bin

(bin type=1 is used here, representing latitudinal dependence) in order to allow representation

of differences between, for example, polar, mid-latitude, and tropical dynamical and physical

processes. In addition, the scalar coefficient c used to estimate velocity potential errors from

those of streamfunction is calculated as a function of height to represent, for example, the

impact of boundary-layer physics. Latitudinal/height smoothing of the resulting coefficients

may be optionally performed to avoid artificial discontinuities at the edges of latitude/height

boxes.

Having computed regression coefficients, the unbalanced components of the fields are calculated as χu = χ − cψ, Tu,k1 = Tk1 −

P

k2 Gk1,k2ψk2, and psu = ps −

P

k Wkψk. These fields are

output for the subsequent calculation of the spatial covariances as described below.

9.3.3 Vertical Covariances: Eigenvectors/eigenvalues and control variable projections

The third stage (gen be stage3) of gen be calculates the statistics required for the vertical component of the control variable transform. This calculation involves the projection of 3D fields

on model-levels onto empirical orthogonal functions (EOFs) of the vertical component of background error covariances (Barker et al., 2004). For each 3D control variable (ψ, χu, Tu, and r),

the vertical component of B, is calculated and an eigenvector decomposition performed. The

resulting eigenvectors E and eigenvalues Λ are saved for use in WRF-Var.

The gen be code calculates both domain-averaged and local values of the vertical component

of the background error covariance matrix. The definition of local again depends on the value

of the namelist variable bin type chosen. For example, for bin type=1, a kz × kz (where kz

is the number of vertical levels) vertical component of B is produced at every latitude (data

is averaged over time and longitude) for each control variable. Eigendecomposition of the resulting climatological vertical error covariances B = EΛET

results in both domain-averaged

and local eigenvectors E and eigenvalues Λ. Both sets of statistics are included in the dataset

supplied to WRF-Var, allowing the choice between homogeneous (domain-averaged) or local

(inhomogeneous) background error variances and vertical correlations to be chosen at run time

73
(Barker et al., 2004). Having calculated and stored eigenvectors and eigenvalues, the final part

of gen be stage3 is to project the entire sequence of 3D control variable fields into EOF space

vv = U

−1

v vp = Λ−1/2ET vp.

9.3.4 Horizontal Covariances: Recursive filter lengthscale (regional),

or power spectra (global)

The last aspect of the climatological component of background error covariance data required

for WRF-Var is the horizontal error correlations, the representation of which forms the largest

difference between running WRF-Var in regional and global mode. (It is however, still a fairly

local change.)

In a global application (gen be stage4 global), power spectra are computed for each of the kz

vertical modes of the 3D control variables ψ, χu, Tu, and r, and for the 2D control variable psu

data. In contrast, in regional mode, horizontal correlations are computed between grid-points

of each 2D field, binned as a function of distance. A Gaussian curve is then fitted to the data

as described in Barker et al. (2004) to provide correlation lengthscales for use in the recursive

filter algorithm.

74
Appendix A

Physical Constants

The following is a list of physical constants used in the model.

π = 3.1415926 Pi

k = 0.4 Von Karman constant

re = 6.370 × 106 m Radius of earth

g = 9.81 m s−2 Acceleration due to gravity

Ωe = 7.2921 × 10−5s

−1 Angular rotation rate of the earth

σB = 5.67051 × 10−8 W m−2 K−4 Stefan − Boltzmann constant

Rd = 287 J kg−1 K−1 Gas constant for dry air

Rv = 461.6 J kg−1 K−1 Gas constant for water vapor

cp = 7 × Rd/2 J kg−1 K−1 Specific heat of dry air at constant pressure

cv = cp − Rd J kg−1 K−1 Specific heat of dry air at constant volume

cpv = 4 × Rv J kg−1 K−1 Specific heat of water vapor at constant pressure

cvv = cpv − Rv J kg−1 K−1 Specific heat of water vapor at constant volume

cliq = 4190 J kg−1 K−1 Specific heat capacity of water

cice = 2106 J kg−1 K−1 Specific heat capacity of ice

Lv = 2.5 × 106 J kg−1 Latent heat of vaporization

Ls = 2.85 × 106 J kg−1 Latent heat of sublimation

Lf = 3.50 × 105 J kg−1 Latent heat of fusion

ρw = 1.0 × 103 kg m−3 Density of liquid water

75
76
Appendix B

List of Symbols

Symbols used in this document are listed in alphabatical order in this appendix.

Symbols Definition

a generic variable

A coefficient (Chapter 4), base-state lapse rate constant (Chapter 5)

B background error covariance matrix

c scalar coefficient

cs speed of sound

Ck a constant used in TKE closure

Cr Courant number

Crmax maximum Courant number

Crtheory Courant number from Table 3.1

Crβ activation Courant number in vertical velocity damping

Cs a constant used in eddy viscosity calculation

D deformation

Dnm deformation tensor, where n, m = 1, 2 and 3

e cosine component of the Coriolis term (Chapters 2, 3); turbulent kinetic energy

(Chapter 4)

E observation error covariance matrix

f sine component of the Coriolis term

F forcing terms for U, V , W, Θ and Qm

F representivity error covariance matrix

FXcor Coriolis forcing terms for X = U, V , and W

F1,2 coefficients for weighting functions in specified boundary condition

g acceleration due to gravity

Gk regression coefficient

H observation operator

J cost function

Kdh,dv horizontal and vertical eddy viscosity for gravity wave absorbing layer

Kh,v horizontal and vertical eddy viscosities

77
Symbols Definition

l0 minimum length scale for dissipation

lh,v horizontal and vertical length scales for turbulence

lcr critical length scale for dissipation

L latent heat of condensation

Lx,y periodicity length in x and y

m map scale factor

ns ratio of the RK3 time step to the acoustic time step

N Brunt-V ̈ais ̈al ̈a frequency

p pressure

p

0 perturbation pressure

p0 reference sea-level pressure

ph hydrostatic pressure

pht,hs hydrostatic pressure at the top and surface of the model

pdht,dhs dry hydrostatic pressure at the top and surface of the model

ps surface pressure

Pr Prandtl number

q generic scalar

qc,i,r,s mixing ratios for cloud water, ice, rain water and snow

qm generic mixing ratios for moisture

qv mixing ratio for water vapor

qvs saturation mixing ratio for water vapor

Qm generic coupled moisture variable

r relative humidity

re radius of earth

R remaining terms in equations

Rd gas constant for dry air

Rv gas constant for water vapor

t time

∆t a full model time step

T temperature

T0 reference sea-level temperature

u horizontal component of velocity in x-direction

U coupled horizontal component of velocity in x-direction (Chapters 2, 3, 6, 7);

control variable transform (Chapter 9)

Uh horizontal correlation

Up multivariate covariance

Uv vertical covariance

v horizontal component of velocity in y-direction

v three dimensional vector velocity

V coupled horizontal component of velocity in y-direction

V three dimensional coupled vector velocity

w vertical component of velocity

W coupled vertical component of velocity

78
Symbols Definition

Wk regression coefficient

z height

zd depth of damping layer

ztop height of model top

α inverse density of air

α

0 perturbation inverse density of air

α ̄ inverse density of air for the reference state

αd inverse density of dry air

αr local rotation angle between y-axis and the meridian

β off-centering coefficient for semi-implicit acoustic step

γ ratio of heat capacities for dry air at constant pressure and volume

γd divergence damping coefficient

γe external mode damping coefficient

γg damping coefficient for upper boundary gravity wave absorbing layer

γr Rayleigh damping coefficient

 molecular weight of water over the molecular weight of dry air (Chapter 4);

true background error (Chapter 9)

η terrain-following hydrostatic-pressure vertical coordinate

η ̇ contravariant ‘vertical’ velocity or coordinate velocity

θ potential temperature

θe equivalent potential temperature

θm moist potential temperature

Θ coupled potential temperature

μ hydrostatic pressure difference between surface and top of the model

μ ̄ reference state hydrostatic pressure difference between surface and top of the

model

μd dry hydrostatic pressure difference between surface and top of the model

τ acoustic time (Chapter 3), vertical structure function for Rayleigh damping

(Chapter 4)

τnm stress tensor (Chapter 4) where n.m = 1, 2 and 3

∆τ acoustic time step

φ geopotential (Chapters 2, 3, 5); latitude (Chapter 9)

φ ̄ geopotential for reference state

φ

0 perturbation geopotential

Φ generic prognostic variable (coupled)

ψ generic variable (Chapter 6)

ψ

0

streamfunction increment

χ

0 velocity potential increment

ω same as  ̇η

Ω coupled coordinate velocity

Ωe angular rotation rate of the earth

79
Subscripts/Superscripts Definition

()d dry

()h hydrostatic

()0 base state sea-level constant

() reference state

()0 perturbation from reference state

()t

∗

full value at a Runge-Kutta step

()00 perturbation from Runge-Kutta step value in acoustic steps

80
Appendix C

Acronyms

AFWA Air Force Weather Agency

API Application Program Interface

ARPS Advanced Regional Prediction System

ARW Advanced Research WRF

BUFR Binary Universal Form for the Representation of Meteorological Data

CAPE Convectively Available Potential Energy

CAPS Center for the Analysis and Prediction of Storms

CGM Conjugate Gradient Method

COAMPS Coupled Ocean / Atmosphere Mesoscale Prediction System

COMET Cooperative Program for Operational Meteorology, Education, and Training

DTC Developmental Testbed Center

ECMWF The European Centre for Medium-Range Weather Forecasts

EOF Empirical Orthogonal Function

ESMF Earth System Modeling Framework

FAA Federal Aviation Administration

FGAT First Guess at Appropriate Time

FSL Forecast System Laboratory

GFDL Geophysical Fluid Dynamics Laboratory

GFS Global Forecast System

GRIB Gridded Binary

KMA Korean Meteorological Administration

LSM Land Surface Model

MKS Meter Kilogram Second

MM5 Pennsylvania State / NCAR Mesoscale Model Version 5

MMM Mesoscale and Microscale Meteorology Division

MRF Medium Range Forecast Model

NAM North American Mesoscale Model

NCAR National Center for Atmospheric Research

NCEP National Centers for Environmental Prediction

NFS Non-hydrostatic Forecast System (Central Weather Bureau of Taiwan)

NMM Nonhydrostatic Mesoscale Model

NOAA National Oceanographic and Atmospheric Administration

NRL Navy Research Laboratory

81
NWP Numerical Weather Prediction

OSU Oregon State University

PBL Planetary Boundary Layer

PPM Piecewise Parabolic Method

QNM Quasi Newton Method

RHS Right Hand Side

RRTM Rapid Radiative Transfer Model

RUC Rapid Update Cycle

SI (WRF) Standard Initialization

TKE Turbulent Kinetic Energy

UCAR University Corporation for Atmospheric Research

YSU Yonsei University (Korea)

VAR Variational Assimilation

WRF Weather Research and Forecasting Model

WSF WRF Software Framework

82
Bibliography

Barker, D. M., W. Huang, Y.-R. Guo, and A. Bourgeois, 2003: A Three-Dimensional Variational

(3DVAR) Data Assimilation System For Use With MM5. NCAR Tech Note, NCAR/TN453+STR, 68 pp. [Available from UCAR Communications, P.O. Box 3000, Boulder, CO,

80307.].

Barker, D. M., W. Huang, Y.-R. Guo, A. Bourgeois, and X. N. Xiao, 2004: A Three-Dimensional

Variational Data Assimilation System for MM5: Implementation and Initial Results Mon.

Wea. Rev., 132, 897–914.

Beljaars, A.C.M., 1994: The parameterization of surface fluxes in large-scale models under free

convection, Quart. J. Roy. Meteor. Soc., 121, 255–270.

Betts, A. K., 1986: A new convective adjustment scheme. Part I: Observational and theoretical

basis. Quart. J. Roy. Meteor. Soc., 112, 677–691.

Betts, A. K., and M. J. Miller, 1986: A new convective adjustment scheme. Part II: Single column

tests using GATE wave, BOMEX, and arctic air-mass data sets. Quart. J. Roy. Meteor. Soc.,

112, 693–709.

Carpenter, R. L., K. K. Droegemeier, P. R. Woodward, and C. E. Hane, 1990: Application of

the Piecewise Parabolic Method (PPM) to meteorological modeling. Mon. Wea. Rev., 118,

586–612.

Chen, F., and J. Dudhia, 2001: Coupling an advanced land-surface/ hydrology model with the

Penn State/ NCAR MM5 modeling system. Part I: Model description and implementation.

Mon. Wea. Rev., 129, 569–585.

Chen, S.-H., and W.-Y. Sun, 2002: A one-dimensional time dependent cloud model. J. Meteor.

Soc. Japan, 80, 99–118.

Chou M.-D., and M. J. Suarez, 1994: An efficient thermal infrared radiation parameterization

for use in general circulation models. NASA Tech. Memo. 104606, 3, 85pp.

Cooper, W. A., 1986: Ice initiation in natural clouds. Precipitation Enhancement — A Scientific

Challenge, Meteor. Monogr., No. 43, Amer. Met. Soc., 29–32.

Daley, R., and E. Barker, 2001: NAVDAS: Formulation and Diagnostics. Mon. Wea. Rev., 129,

869–883.

83
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.

Deardorff, J. W., 1972 Numerical investigation of neutral and unstable planetary boundary

layers, J. Atmos. Sci., 29, 91–115.

Desroziers, G., 1997: A coordinate change for data assimilation in spherical geometry of frontal

structures. Mon. Wea. Rev., 125, 3030–3039.

Desroziers, G., and S. Ivanov, 2001: Diagnosis and adaptive tuning of observation-error parameters in a variational assimilation. Quart. J. Roy. Meteor. Soc., 127, 1433–1452.

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.

Dudhia, J., 1995: Reply, Mon. Wea. Rev., 123, 2571–2575.

Durran, R. D., and J. B. Klemp, 1983: A compressible model for the simulation of moist

mountain waves, Mon. Wea. Rev., 111, 2341–2361.

Dyer, A. J., and B. B. Hicks, 1970: Flux-gradient relationships in the constant flux layer, Quart.

J. Roy. Meteor. Soc., 96, 715–721.

Easter, R. C., 1993: Two modified versions of Bott’s positive definite numerical advection

scheme. Mon. Wea. Rev., 121, 297–304.

Fels, S. B. and M. D. Schwarzkopf, 1975: The Simplified Exchange Approximation: A New

Method for Radiative Transfer Calculations, J. Atmos. Sci., 32, 1475–1488.

Fisher, M., 2003: Background error covariance modeling. Seminar on Recent Development in

Data Assimilation for Atmosphere and Ocean, 45–63, ECMWF.

Fletcher, N. H., 1962: The Physics Of Rain Clouds. Cambridge University Press, 386 pp.

Grell, G. A., and D. Devenyi, 2002: A generalized approach to parameterizing convection combining ensemble and data assimilation techniques. Geophys. Res. Lett., 29(14), Article 1693.

Haltiner, G. J., and R. T. Williams, 1980: Numerical prediction and dynamic meteorology. John

Wiley & Sons, Inc., 477pp.

Hollingsworth, A., and P. Lonnberg, 1986: The Statistical Structure Of Short-Range Forecast

Errors As Determined From Radiosonde Data. Part I: The Wind Field. Tellus, 38A, 111–136.

Hong, S.-Y., and H.-L. Pan, 1996: Nonlocal boundary layer vertical diffusion in a medium-range

forecast model, Mon. Wea. Rev., 124, 2322–2339.

Hong, S.-Y., H.-M. H. Juang, and Q. Zhao, 1998: Implementation of prognostic cloud scheme

for a regional spectral model, Mon. Wea. Rev., 126, 2621–2639.

84
Hong, S.-Y., J. Dudhia, and S.-H. Chen, 2004: A Revised Approach to Ice Microphysical Processes for the Bulk Parameterization of Clouds and Precipitation, Mon. Wea. Rev., 132,

103–120.

Ide, K., P. Courtier, M. Ghil, and A. C. Lorenc, 1997: Unified notation for data assimilation:

Operational, sequential and variational. J. Met. Soc. Japan, 75, 181–189.

Ingleby, N. B., 2001: The statistical structure of forecast errors and its representation in the

Met. Office global 3-D variational data assimilation scheme. Quart. J. Roy. Meteor. Soc., 127,

209–232.

Janjic, Z. I., 1990: The step-mountain coordinate: physical package, Mon. Wea. Rev., 118,

1429–1443.

Janjic, Z. I., 1994: The step-mountain eta coordinate model: further developments of the convection, viscous sublayer and turbulence closure schemes, Mon. Wea. Rev., 122, 927–945.

Janjic, Z. I., 1996: The surface layer in the NCEP Eta Model, Eleventh Conference on Numerical

Weather Prediction, Norfolk, VA, 19–23 August; Amer. Meteor. Soc., Boston, MA, 354–355.

Janjic, Z. I., 2000: Comments on ”Development and Evaluation of a Convection Scheme for Use

in Climate Models”, J. Atmos. Sci., 57, p. 3686.

Janjic, Z. I., 2002: Nonsingular Implementation of the Mellor–Yamada Level 2.5 Scheme in the

NCEP Meso model, NCEP Office Note, No. 437, 61 pp.

Kain, J. S., and J. M. Fritsch, 1990: A one-dimensional entraining/ detraining plume model

and its application in convective parameterization, J. Atmos. Sci., 47, 2784–2802.

Kain, J. S., and J. M. Fritsch, 1993: Convective parameterization for mesoscale models: The

Kain-Fritcsh scheme, The representation of cumulus convection in numerical models, K. A.

Emanuel and D.J. Raymond, Eds., Amer. Meteor. Soc., 246 pp.

Kessler, E., 1969: On the distribution and continuity of water substance in atmospheric circulation, Meteor. Monogr., 32, Amer. Meteor. Soc., 84 pp.

Klemp, J. B., and D. K. Lilly, 1978: Numerical simulation of hydrostatic mountain waves, J.

Atmos. Sci., 35, 78–107.

Klemp, J. B., and R. Wilhelmson, 1978: The simulation of three-dimensional convective storm

dynamics, J. Atmos. Sci., 35, 1070–1096.

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.

Laprise R., 1992: The Euler Equations of motion with hydrostatic pressure as as independent

variable, Mon. Wea. Rev., 120, 197–207.

Lee, M.-S., D. Barker, W. Huang and Y.-H. Kuo, 2004: First Guess at Appropriate Time

(FGAT) with WRF 3DVAR. WRF/MM5 Users Workshop, 22–25 June 2004, Boulder, Colorado.

85
Lin, Y.-L., R. D. Farley, and H. D. Orville, 1983: Bulk parameterization of the snow field in a

cloud model. J. Climate Appl. Meteor., 22, 1065–1092.

Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux-form semi-Lagrangian transport

schemes. Mon. Wea. Rev., 124, 2046–2070.

Lorenc, A. C., 1986: Analysis methods for numerical weather prediction. Quart. J. Roy. Meteor.

Soc., 112, 1177–1194.

Mellor, G. L., and T. Yamada, 1982: Development of a turbulence closure model for geophysical

fluid problems. Rev. Geophys. Space Phys., 20, 851–875.

Michalakes, J., J. Dudhia, D. Gill, J. Klemp, and W. Skamarock, 1999: Design of a nextgeneration regional weather research and forecast model, Towards Teracomputing, World Scientific, River Edge, New Jersey, 117–124.

Michalakes, J., J. Dudhia, D. Gill, T. Henderson, J. Klemp, W. Skamarock, and W. Wang,

2004: The Weather Research and Forecast Model: Software Architecture and Performance.

To appear in proceeding of the Eleventh ECMWF Workshop on the Use of High Performance

Computing in Meteorology. 25–29 October 2004, Reading, U.K., Ed. George Mozdzynski.

Mlawer, E. J., S. J. Taubman, P. D. Brown, M. J. Iacono, and S. A. Clough, 1997: Radiative

transfer for inhomogeneous atmosphere: RRTM, a validated correlated-k model for the longwave. J. Geophys. Res., 102 (D14), 16663–16682.

Monin, A.S. and A.M. Obukhov, 1954: Basic laws of turbulent mixing in the surface layer of

the atmosphere. Contrib. Geophys. Inst. Acad. Sci., USSR, (151), 163–187 (in Russian).

Ooyama K. V., 1990: A thermodynamic foundation for modeling the moist atmosphere, J.

Atmos. Sci., 47, 2580–2593.

Parrish, D. F., and J. C. Derber, 1992: The National Meteorological Center’s Spectral Statistical

Interpolation analysis system. Mon. Wea. Rev., 120, 1747–1763.

Paulson, C. A., 1970: The mathematical representation of wind speed and temperature profiles

in the unstable atmospheric surface layer. J. Appl. Meteor., 9, 857–861.

Purser, R. J., W. -S. Wu, D. F. Parrish, and N. M. Roberts, 2003: Numerical aspects of the

application of recursive filters to variational statistical analysis. Part I: Spatially homogeneous

and isotropic Gaussian covariances. Mon. Wea. Rev., 131, 1524–1535.

Rabier, F., H., J. N. Thepaut, and P. Courtier, 1998: Extended assimilation and forecast

experiments with a four-dimensional variational assimilation system. Quart. J. Roy. Meteor.

Soc., 124, 1861–1887.

Reisner, J. R., R. M. Rasmussen, and R. T. Bruintjes, 1998: Explicit forecasting of supercooled

liquid water in in winter storms using the MM5 mesoscale model. Quart. J. Roy. Meteor.

Soc., 124, 1071–1107.

86
Roberts, R. E., J. E. A. Selby, and L. M. Biberman, 1976: Infrared continuum absorption by

atmospheric water-vapor in 8–12 um window. Applied Optics, 15 (9), 2085–2090.

Rodgers, C. D., 1968: Some extensions and applications of the new random model for molecular

band transmission. Quart. J. Roy. Meteor. Soc., 94, 99–102.

Rutledge, S. A., and P. V. Hobbs, 1984: The mesoscale and microscale structure and organization

of clouds and precipitation in midlatitude cyclones. XII: A diagnostic modeling study of

precipitation development in narrow cloud-frontal rainbands. J. Atmos. Sci., 20, 2949–2972.

Ryan, B. F., 1996: On the global variation of precipitating layer clouds. Bull. Amer. Meteor.

Soc., 77, 53–70.

Sasamori, T., J. London, and D. V. Hoyt, 1972: Radiation budget of the Southern Hemisphere.

Meteor. Monogr., 13, No. 35, 9–23.

Schwarzkopf, M. D., and S. B. Fels, 1985: Improvements to the algorithm for computing CO2

transmissivities and cooling rates. J. Geophys. Res., 90 (ND6), 541–550.

Schwarzkopf, M. D., and S. B. Fels, 1991: The simplified exchange method revisited — An

accurate, rapid method for computation of infrared cooling rates and fluxes. J. Geophys.

Res., 96 (D5), 9075–9096.

Skamarock W. C. and J. B. Klemp, 1992: The Stability of Time-Split Numerical Methods for

the Hydrostatic and the Nonhydrostatic Elastic Equations, Mon. Wea. Rev.: 120, 2109–2127.

Skamarock W. C., 2005: Timesplitting techniques for multidimensional transport, available at

http://www.mmm.ucar.edu/individual/skamarock/advect3d 2005.pdf, 26 pp.

Smirnova, T. G., J. M. Brown, and S. G. Benjamin, 1997: Performance of different soil model

configurations in simulating ground surface temperature and surface fluxes. Mon. Wea. Rev.,

125, 1870–1884.

Smirnova, T. G., J. M. Brown, S. G. Benjamin, and D. Kim, 2000: Parameterization of coldseason processes in the MAPS land-surface scheme. J. Geophys. Res., 105 (D3), 4077–4086.

Smolarkiewicz, P. K., and G. A. Grell, 1990: A class of monotone interpolation schemes, J.

Comp. Phys., 101, 431–440.

Stephens, G. L., 1978: Radiation profiles in extended water clouds. Part II: Parameterization

schemes, J. Atmos. Sci., 35, 2123–2132.

Tao, W.-K., J. Simpson, and M. McCumber 1989: An ice-water saturation adjustment, Mon.

Wea. Rev., 117, 231–235.

Thompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit forecasts of winter precipitation using an improved bulk microphysics scheme. Part I: Description and sensitivity analysis.

Mon. Wea. Rev., 132, 519–542.

87
Walko, R. L., W. R. Cotton, M. P. Meyers, and J. Y. Harrington, 1995: Newe RAMS cloud

microphysics parameterization. Part I: The single-moment scheme. Atmos. Res., 38, 29–62.

Wang, W., D. Barker, C. Bruy`ere, J. Dudhia, D. Gill, and J. Michalakes, 2004: WRF Version

2 modeling system user’s guide. http://www.mmm.ucar.edu/wrf/users/docs/user guide/.

Webb, E. K., 1970: Profile relationships: The log-linear range, and extension to strong stability,

Quart. J. Roy. Meteor. Soc., 96, 67–90.

Wicker, L. J. and W. C. Skamarock, 2002: Time splitting methods for elastic models using

forward time schemes, Mon. Wea. Rev., 130, 2088–2097.

Wicker, L. J., and R. B. Wilhelmson, 1995: Simulation and analysis of tornado development

and decay within a three-dimensional supercell thunderstorm. J. Atmos. Sci., 52, 2675–2703.

Wu, W. -S., R. J. Purser, and D. F. Parrish, 2002: Three-Dimensional Variational Analysis with

Spatially Inhomogeneous Covariances. Mon. Wea. Rev., 130, 2905–2916.

Xiao, Q. N., Y. H. Kuo, J. Sun, W. C. Lee, E. Lim, Y. R. Guo, and D. M. Barker, 2005:

Assimilation of Doppler Radar Observations with a Regional 3D–Var System: Impact of

Doppler Velocities on Forecasts of a Heavy Rainfall Case. J. Appl. Met., 44(6), 768–788.

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.

Zilitinkevich, S. S., 1995: Non-local turbulent transport: pollution dispersion aspects of coherent structure of convective flows, Air Pollution III — Volume I. Air Pollution Theory and

Simulation, Eds. H. Power, N. Moussiopoulos and C.A. Brebbia. Computational Mechanics

Publications, Southampton Boston, 53–60.

88