Documentation of the PMIP models (Bonfils et al. 1998)

PMIP Documentation for UIUC11

University of Illinois at Urbana-Champaign: Model UIUC MLAM-PMIP (4x5 L11) 1996

PMIP Representative(s)



Dr. Michael Schlesinger, Department of Atmospheric Sciences, University of Illinois at Urbana-Champaign, 105 South Gregory Avenue, Urbana, Illinois 61801; Phone: +1 217 333 2192; Fax: +1-217 244 4393; e-mail:

World Wide Web URL:

Model Designation

UIUC MLAM-PMIP (4x5 L11) 1996

Model Identification for PMIP


PMIP run(s)

0fix, 6fix

Number of days in each month: 31 28 31 30 31 30 31 31 30 31 30 31

Model Lineage

The 11-layer tropospheric/lower-stratospheric GCM is a descendent of the 2-layer atmospheric GCM developed in the late 1960's and early 1970's by Arakawa and Mintz at UCLA , and subsequently developed and used by Schlesinger during the past 24 years, first at the Rand Corporation from 1973 to 1976, second at Oregon State University (OSU) from 1976 to 1989 , and since 1989 at the University of Illinois at Urbana-Champaign. These modifications principally include an increase in vertical resolution from 2 to 11 layers and substantial changes in the treatment of atmospheric radiation, convection, cloud/precipitation formation, and land surface processes.


Model Documentation

main reference:

Schlesinger, M. E., N. G. Andronova, B. Entwistle , A. Ghanem , N. Ramankutty, W. Wang and F. Yang, 1997:

Modeling and Simulation of Climate and Climate Change, In Past and present variability of the

solar-terrestrial system: Measurements, data analysis and theoretical models, Proceedings of the

International School of Physics "Enrico Fermi", Course CXXXIII, 25 June - 5 July 1996, Varrena, Italy. G. Cini Castagnoli and A. Provenzale (Eds.), IOS Press, Amsterdam, p389-429.

secondary reference(s):

The dynamical structure and numerics of the UIUC model, as well as some of its surface schemes are as described in:

Ghan, S. J., J. W. Lingaas, M. E. Schlesinger, R. L. Mobley, and W. L. Gates, A documentation of the OSU two-level atmospheric general circulation model . 1982, Climatic Research Institute, Oregon State University: Corvallis, OR. p. 395.

The parameterizations of radiation, cloud formation, and related physics are discussed in:

Oh, J.-H., Physically-Based General Circulation Model Parameterization of Clouds and their Radiative Interaction in Department of Atmospheric Sciences, 1989, Oregon State University: Corvallis. p. 315.

Numerical/Computational Properties

Horizontal Representation

The horizontal distribution of dependent variables in the model is staggered according to the B-grid to simulate the process of geostrophic adjustment , and the model uses finite differences that conserve the total atmospheric mass, total energy under adiabatic and frictionless motion, and enstrophy (mean-square vorticity) for the nondivergent component of the wind field .

B-grid: Surface pressure, temperature and water-vapor are located in the center of each grid cell, and the horizontal velocity components are located at the corners of the cell.

Horizontal Resolution

4 x 5-degree latitude-longitude grid.

dim_longitude*dim_latitude: 72*46.

Vertical Domain

Surface to 50 mb (sigma=0.000) (model top). For a surface pressure of 1000 hPa, the lowest prognostic level is at 990 hPa and the highest is at 60 hPa.

Vertical Representation

Finite-difference sigma coordinates.

The model uses normalized pressure, sigma, as its vertical coordinate, such that the earth's surface is the coordinate surface sigma=1 and the top of the model is the coordinate surface sigma=0.

s = (p - pT)/(pS - pT), where p is pressure, pT is the constant pressure at the top of the model, and pS is the surface pressure, which varies with horizontal location and time.

Vertical Resolution

There are 11 unevenly spaced sigma layers between the surface and the model top at 50 hPa. The centers of the 11 layers are located at 0.0103, 0.0367, 0.0736, 0.1258, 0.2292, 0.3775, 0.5254, 0.6726, 0.8203, 0.9370, 0.9897.

For a surface pressure of 1000 hPa, 2 levels are below 900 hPa and 4 levels are above 200 hPa.

Computer/Operating System

For the PMIP simulation, the model was run on a Y_MP computer under the UNICOS operating system for 0fix and the model was run on a Cray C90 under the UNICOS operating system for 6fix.

Computational Performance

For the PMIP experiment, the time integration of the governing equations requires about 13 hours of Cray Y_MP computer time per simulated year.


Initial conditions of PMIP model for control simulation were taken from a control simulation of the UIUC 11-Layer AGCM.

The model was started from an equilibrium state of a 1*CO2 experiment.

Time Integration Scheme(s)

The adiabatic, frictionless terms in the governing prognostic equations are marched forward in time in a one-hour cycle using a combination of Matsuno and leapfrog time integration steps, with a basic timestep of 6 minutes.

The adiabatic, frictional terms in the governing prognostic equations are calculated once per hour.


Orography is area-averaged on the model grid (see Orography). A longitudinal smoothing of the zonal pressure gradient and the zonal mass flux is performed polewards of 38 degrees latitude . It is unnecessary to fill spurious negative values of atmospheric moisture, since these are not generated by the numerical schemes.

Sampling Frequency

For the PMIP simulation, the model history is written every six hours.

Dynamical/Physical Properties

Atmospheric Dynamics

The dynamics of the Primitive equations are expressed in terms of u and v winds, temperature, and surface pressure. Specific humidity and cloud water are also a prognostic variables (see Cloud Formation).


Horizontal diffusion is not modeled.

Vertical diffusion of momentum, sensible heat, and moisture operates at all vertical levels. The diffusion depends on the vertical wind shear, but not on stability .

Gravity-wave Drag

Gravity-wave drag is not modeled.

A momentum drag is included in the top layer of the model that is proportional to air density and the square of the velocity (Hansen et. al 1983).

Solar Constant/Cycles

The solar constant is the PMIP-prescribed value of 1365 W/(m2). The orbital parameters and seasonal insolation distribution are calculated after PMIP recommendations. Both seasonal and diurnal cycles in solar forcing are included.


The carbon dioxide concentration is the PMIP-prescribed value of 345 ppm and 280 ppm for 0fix and 6fix, respectively. The daily horizontal distribution of column-integrated ozone is interpolated from prescribed monthly mean Total Ozone Mapping Spectrometer (TOMS) data . The radiative effects of water vapor , ozone , carbon dioxide , methane and nitrous oxide , and CFC-11 and CFC-12 . Direct radiative forcing by aerosols on solar radiation is included (see Radiation).




Radiative Transfer. - Upward and downward fluxes of solar (shortwave) radiation are calculated using a two-stream method with the delta-Eddington approximation. Scattering and absorption by both gases and cloud droplets are included. The spectrum of solar radiation is divided into three intervals: (1) 0 to 0.44 µm (2) 0.44 to 0.69 µm and (3) 0.69 to 3.85 µm.The first two intervals are for the treatment of Rayleigh scattering and absorption by ozone and carbon dioxide . The third interval is for water-vapor absorption and is subdivided into six k-distribution intervals . The optical depth and single-scattering albedo for cloud droplets are determined separately for non-ice clouds and for cirrus clouds .

Upward and downward fluxes of longwave (terrestrial) radiation are calculated using a two-stream formulation. Absorption is treated in four spectral bands: (1) 0 - 340 cm-1 and 1380 - 1900 cm-1, for the water-vapor line centers; (2) 340 - 540, 800 - 980, 1100 - 1380, and 1900 - 3000 cm-1 for the water-vapor line wings and part of the water-vapor continuum; (3) 540 - 800 cm-1, for carbon dioxide and part of the water-vapor continuum; and (4) 980 - 1100 cm-1, for ozone and part of the water-vapor continuum. Absorption by the trace gases methane, nitrous oxide, and chlorofluorocarbon compounds CFC-11 and CFC-12 is included in spectral bands 2-4. Pressure-broadening effects are included in the absorption calculations for water vapor , ozone , carbon dioxide , methane and nitrous oxide , and CFC-11 and CFC-12 . Absorption by cloud droplets is treated by an emissivity formulation for non-ice clouds , extratropical clouds and tropical cirrus clouds .

The radiation parameterization includes cloud-cover feedback by calculating separately the radiative fluxes for the cloudy and clear portions of each grid cell, and includes cloud optical-depth feedback by linking the radiative properties to the prognostic cloud-water content. Clouds located in contiguous vertical layers comprise a cloud group. The contiguous cloud layers within each group are overlapped fully in the vertical. The noncontiguous cloud groups, separated from each other by at least one layer of clear air, are overlapped randomly .

The directive radiative forcing of sulfate aerosol due to its backscattering of incoming solar radiation is included by the model. It is assumed that the sulfate aerosol particles consist of 75% H2SO4 and 25% H2O, and that their size distribution is log-normal, with a geometric mean radius of 0.05 µm and a standard deviation of 2.0 µm . The optical properties of the sulfate aerosol - extinction efficiency, single-scattering albedo, and asymmetry factor - are computed once and for all outside the GCM using our Mie scattering model.


The model has three parameterizations for convection: (1) dry-convective adjustment, (2) middle-level convection, and (3) penetrating convection.

Dry-convective adjustment occurs if the temperature lapse rate between any two adjacent vertical layers is absolutely unstable, that is, exceeds the dry-adiabatic lapse rate. If this occurs, the instability is instantaneously removed by adjusting the temperatures of the two layers such that their lapse rate is dry adiabatic. This is done by transferring heat vertically between the layers under the constraint that their total enthalpy is conserved. Dry-convective adjustment is performed from the lowest to the highest model layer, iteratively.

Middle-level convection occurs if the temperature lapse rate between any two adjacent vertical layers is conditionally unstable and the lower-layer air is sufficiently near saturation that it would be positively buoyant if displaced to the higher layer . This condition occurs when the moist static energy of the lower layer exceeds the saturated moist static energy of the upper layer. When the instability exists, an upward convective mass flux occurs between the layers within a convective tower, and a compensating downward mass flux occurs between the layers in the environment outside the convective tower. Because the air within the convective tower is saturated, the convective mass flux therein generates liquid water, part of which is converted into convective precipitation that falls out of the cloud. The subsiding mass flux in the environment modifies the environmental temperature, water vapor and horizontal momentum. This modification of the environment reduces the instability at a rate that depends on the convective mass flux. The latter is calculated such that the instability is removed with an e-folding time of one hour. The fractional cloudiness for middle-level convective cloud is a function of the convective mass flux and the relative humidity of the higher layer.

Penetrating convection occurs if the temperature lapse rate between the PBL and any layer above is conditionally unstable and the PBL air is sufficiently near saturation that it would be positively buoyant if displaced to the higher layer . This condition occurs when the moist static energy of the PBL exceeds the saturated moist static energy of the layer above. Also, the relative humidity of the PBL must equal or exceed 95 percent for penetrative convection to be initiated . The treatment of penetrating convection is essentially the same as the treatment for middle-level convection, except that: (1) as many convective towers may coexist as there are layers above the PBL, one tower extending from the PBL to each layer for which the instability exists; (2) environmental air is entrained into each convective tower from all layers through which it passes, and this mass transport modifies the temperature and water vapor within the cloud, as well as the temperature, water vapor and horizontal momentum within the environment; and (3) when the initiating instability for any cloud tower ceases to exist, the cloudiness at its top level evaporates with a prescribed e-folding time.

Large-scale condensation occurs in a layer not only when the cell is everywhere saturated, but also when only part of the grid cell is saturated . The rate of condensation depends on the large-scale convergence rates of moisture, heat and mass, and the time rate of change of fractional relative humidity of the layer, U. The latter is determined from U = bUs + (1 - b)Uo, where b is the fractional cloudiness, Us (=1) is the saturated relative humidity within the cloud, (1 - b) is the cloud-free fraction, and Uo the fractional relative humidity of the clear air. Closure is achieved by assuming: (1) the moisture convergence is partitioned between the cloud and clear air in proportion to b and 1-b, respectively; and (2) Uo = Uoo + b(Uo - Uoo), where Uoo is the relative humidity at which condensation can begin. The result is that b = 1 - {(Us - U)/(Us - Uoo)}1/2, which increases from zero for U = Uoo to unity for U = Us = 1. Uoo is taken to be 0.90.

Cloud Formation

Stratiform and cumuloform clouds can coexist within the same vertical atmospheric column, albeit not in the same layer. A cloud in any vertical layer is identified as either a stratiform or cumuloform cloud depending on the preceding cloud type, the large-scale condensation and the convective mass flux (see Convection) in the layer. If there is convective mass flux, the cloud type is taken to be cumuloform regardless of whether the preceding cloud was stratiform or cumuloform. If there is large-scale condensation and no convection, the cloud is taken to be stratiform. If there is neither convection nor large-scale condensation, the cloud maintains its cloud type until it dissipates by evaporation .

Cloud in the PBL (see Planetary Boundary Layer) is diagnostically computed on the basis of a cloud-topped mixed-layer model .


For clouds with temperatures below 0oC, a fraction of the cloud water is taken to be ice, with the fraction increasing linearly from zero at 0oC to unity at -15oC. Precipitation occurs in the ice phase if the cloud temperature is less than 0oC, and in the liquid phase otherwise. The rate of cloud-water conversion to precipitation is ten times larger for the ice phase than for the liquid phase. Large-scale precipitation beneath cloud base evaporates (sublimates) at a rate that is proportional to the product of the precipitation rate, the relative humidity deficit from saturation, and the cloud-free fraction of the grid cell. Convective precipitation beneath cloud base evaporates at a rate that is proportional to the product of the relative humidity deficit and the cloud-water content .

Planetary Boundary Layer

The top of the planetary boundary layer (PBL) is taken to be the height of the lowest two atmospheric layers. Cloud in the PBL is diagnostically computed on the basis of a cloud-topped mixed-layer model of Lilly and Guinn and Schubert .

See also Cloud Formation, Diffusion, Surface Characteristics and Surface Fluxes


Orography, obtained from the 1 x 1-degree data of Gates and Nelson (1975) is area-averaged over each 4 x 5-degree model grid square.


For control: PMIP datasets have been used for SSTs and sea-ice have been prepared at PCMDI and were calculated by averaging the 10-year AMIP datasets (1979-1988).

For 6 fix: SSTs and sea-ice prescribed at their present day value, as in the control run.

Sea Ice

Sea ice is prescribed for both 0fix and 6fix, as described above.

Snow Cover

The snow mass is determined from a prognostic budget equation that includes the rates of snow accumulation, melting and sublimation. Precipitation falls as snow if the temperature of the lowest model layer is less than 0°C. The snowmelt rate is computed over land from the difference between the downward heat fluxes at the surface and the upward heat fluxes that would occur for a ground temperature equal to the melting temperature of snow (0°C). Snowmelt contributes to the soil moisture (see Land Surface Processes). Accumulation and melting of snow may also occur on sea ice (see Sea Ice). The surface sublimation rate is equated to the evaporative flux from snow, unless sublimation removes all the snow mass (see Surface Fluxes) in less than one hour, in which case the sublimation rate is set equal to the snow-mass removal rate.

Snow cover also alters the surface albedo (see Surface Characteristics).

Surface Characteristics

Surface roughness is specified as in Hansen et al. . Over land, the roughness length is a fit to the data of Fiedler and Panofsky as a function of the standard deviation of the orography. The maximum of this value and that of the roughness of the local vegetation, including a "zero plane displacement" value for tall vegetation types determines the roughness length over land. Over sea ice, the roughness is a constant 4.3 x 10-4 m after Doronin . Over ocean, the roughness length is a function of the surface wind speed, following Garratt .

Snow-free surface albedo is updated monthly by interpolation using values for January, April, July, and October specified from data of Matthews . The albedo of snow-covered surfaces is determined as a linear weighted (by snow depth) interpolation of snow-free and snow-covered values. The albedo of snow is a function of its temperature ; it also depends on solar zenith angle , but not on spectral interval.

Longwave emissivity is specified to be unity (blackbody emission) for all surfaces.

Surface Fluxes

The absorbed surface solar flux is determined from the surface albedo, and surface longwave emission from the Planck function with constant surface emissivity of 1.0 (see Surface Characteristics).

The turbulent surface fluxes of momentum, sensible heat and moisture are parameterized by bulk formulas that depend on the differences of the momentum, temperature and moisture between the ground and surface air, the surface-air wind speed, and aerodynamic drag and transfer coefficients. The surface-air wind is taken as a fraction of the winds extrapolated from the lowest two model layers. Under the assumption that the PBL is well mixed in potential temperature and moisture, the surface-air temperature is extrapolated from the temperature at the lowest level (about 80 meters) with the dry-adiabatic lapse rate (9.8°C/km) and the surface-air moisture is taken to be the same as that at the lowest atmospheric level. The aerodynamic drag and transfer coefficients depend on the vertical stability and surface roughness length, with the same transfer coefficient used for the fluxes of sensible heat and moisture. The surface moisture flux (see Land Surface Processes) depends on an evapotranspiration efficiency, taken as unity over snow, ice and water, and as a function of the soil wetness over land.

Land Surface Processes

The ground temperature is taken to be the average temperature over the diurnal skin depth, calculated from a prognostic budget equation whose source and sink terms include the surface fluxes of radiation, sensible heat and latent heat, and the heat transfer into the ground (see Surface Fluxes). The latter depends on the thermal conductivity and bulk heat capacity of the ground.

The soil wetness is determined from a prognostic budget equation that includes the rates of precipitation, snowmelt, surface evaporation and runoff. Soil wetness is the ratio of the soil moisture content to the field capacity, the latter prescribed for each of the 35 combinations of the AGCM?s 5 soil textures (Sandy; sandy loam; light loam; loamy; and heavy loam, clay) and 7 surface types ((1) Evergreen wood and forest; (2) mixed and deciduous wood, and forest; (3) grassland; (4) cropland; (5) shrub and semi-desert; (6) desert; and (7) tundra, mountain, arctic flora) after data of Vinnikov and Yeserkepova . The evapotranspiration efficiency over land is taken as the minimum of 4/3 the soil wetness and unity. The runoff rate is a nonlinear function of the soil wetness and the combined rates of precipitation and snowmelt. If the predicted soil wetness exceeds unity, the excess moisture is taken as additional runoff.

Last update November 9, 1998. For further information, contact: Céline Bonfils ( )