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

PMIP Documentation for CSIRO

Commonwealth Scientific and Industrial Research Organization: Model CSIRO V4-7 R21L9 199?

PMIP Representative(s)

Dr. Jozef Syktus, CSIRO Division of Atmospheric Research, Private Bag No 1, Aspendale, 3195 Victoria Australia, Phone: 61 3 9239 4548; Fax: 61 3 9239 4444; e-mail:


Model Designation

CSIRO V4-7 R21L9 199?

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 CSIRO model is derived from earlier two-level and four-level spectral models based on the primitive equations expressed in conservative flux form (cf. Gordon 1981 ,1993 ).

Model Documentation

main reference:

H. B. Gordon and S. P. O'Farrell, Transient Climate Change in the CSIRO Coupled Model with Dynamic Sea Ice , 1997, Mon. Wea. Rev., 125 : 875-907

secondary reference(s)

McGregor, J. L., H. B. Gordon, I. G. Watterson, M. R. Dix, and L. D. Rotstayn, 1993, The CSIRO 9-level AGCM. CSIRO Dar Tech. Paper No. 26, CSIRO, PMB1, Aspendale, Victoria 3195, Australia, 89 pp.

Watterson, I. G., Dix, M. R., Gordon, H. B., and McGregor, J. L. (1995). The CSIRO nine-level atmospheric general circulation model and its equilibrium present and doubled CO2 climates. Australian Meteorological Magazine, 44 (2): 111-125.

Documentation of the present version of the CSIRO model is provided by Gordon and O'Farell (1997).

Numerical/Computational Properties

Horizontal Representation

Spectral (spherical harmonic basis functions) with transformation to a Gaussian grid for calculation of nonlinear quantities and some physics. The atmospheric moisture field is represented only in gridded form.

Horizontal Resolution

Spectral rhomboidal 21 (R21), roughly equivalent to 3.2 x 5.6 degrees latitude-longitude.

dim_longitude*dim_latitude: 64*56

Vertical Domain

Surface to about 21 hPa. For a surface pressure of 1000 hPa, the lowest atmospheric level is at about 979 hPa.

Vertical Representation

Finite-difference sigma coordinates.

Vertical Resolution

There are 9 unevenly spaced sigma levels with the following levels: 0.0171 0.0802 0.1927 0.3381 0.5 0.6619 0.8073 0.9198 0.9829. For a surface pressure of 1000 hPa, 3 levels are below 800 hPa and 3 levels are above 200 hPa.

Computer/Operating System

The PMIP simulation was run on a Cray Y/MP computer using one processor in the UNICOS environment.

Computational Performance

For the PMIP experiment, about 35 seconds Cray Y/MP computation time per simulated day.


Initial conditions of PMIP model for control simulation were taken from previous integration of model with fixed SST run for few years. Initialization of the atmosphere, soil moisture, and snow cover/depth for 1 January 1979 is from an earlier model simulation with climatological sea surface temperatures.

Time Integration Scheme(s)

A semi-implicit leapfrog time scheme with an Asselin (1972) frequency filter is used for most calculations, with the momentum surface flux and vertical diffusion above the surface computed by split backward implicit integration. A time step of 30 minutes is used for dynamics and physics, except for full calculations of all radiative fluxes and heating rates, which are done every 2 hours.


Orography is truncated at the R21 resolution of the model (see Orography). To counter the negative values of atmospheric moisture that may otherwise develop, vertical transport of moisture is inhibited if the local water vapor mixing ratio drops below 2 x 10-6 kg (water) per kg (air). In addition, negative moisture values are removed by a proportional adjustment method while conserving the global mean (cf. Royer 1986 ). Cf. McGregor et al. 1993 for further details.

Sampling Frequency

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

Dynamical/Physical Properties

Atmospheric Dynamics

Primitive-equation dynamics are expressed in conservative flux form (i.e., weighting vorticity, divergence, temperature, and specific humidity by the prognostic surface pressure) as described by Gordon (1981) . Effects of frictional heating are included in the temperature tendency equation, and virtual temperature is used to compute geopotential height.


Linear second-order (Ñ2) horizontal diffusion of the temperature, vorticity, divergence, and moisture fields is computed via a split implicit time integration. Diffusion (with coefficient 106 m2/s) is applied to the upper-half of the rhomboid for the spectral temperature, vorticity, and divergence. Diffusion of the gridded moisture is calculated from a temporary spectral representation of the moisture field without surface pressure weighting, and is applied to the entire rhomboid with the diffusion coefficient halved. Diffusion for temperature and moisture contains a first-order correction to constant pressure surfaces.

Vertical diffusion of momentum, heat, and moisture is parameterized in terms of stability-dependent K-theory following Blackadar (1962) , with the choice of asymptotic mixing length after Louis (1979) . The calculation of vertical momentum diffusion is via backward implicit time differencing (see Time Integration Scheme(s)).

Gravity-wave Drag

Under conditions of vertical stability, orographic gravity-wave drag is simulated after the method of Chouinard et al. (1986) . The drag at the surface is dependent on sub-gridscale orographic variance (see Orography), and it is parameterized by means of a "launching" height which is defined to be twice the local standard deviation of the surface heights. Following Palmer et al. (1986) , the maximum launching height is limited to 800 m in order to prevent two-grid noise near steep mountains. At a particular sigma level the frictional drag on the atmosphere from breaking gravity waves depends on the projection of the wind on the surface wind and on the Froude number, which in turn is a function of the launching height, the atmospheric density, the Brunt-Vaisalla frequency, and the wind shear. Gravity-wave drag is assumed to be zero above a critical level, which is taken to be the top sigma level of the model (see Vertical Domain).

Solar Constant/Cycles

The solar constant is the AMIP-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 simulated.


The carbon dioxide concentration is the AMIP-prescribed value of 345 and 280ppm for 0fix and 6fix respectively. Ozone concentrations, specified as a function of latitude and pressure, are interpolated from the Dopplick (1974) seasonal climatology. Radiative effects of water vapor, but not of aerosols, also are included (see Radiation).


The radiation code is after Fels (1985) , Fels and Schwarzkopf (1975 , 1981 ), and Schwarzkopf and Fels (1991). The shortwave calculations are based on a modified Lacis and Hansen (1974) approach. The shortwave spectrum is divided into 9 bands, the first band covering the ultraviolet (wavelengths 0.1 to 0.4 micron) and visible (0.4 to 0.7 micron) spectral intervals, while the other 8 bands are in the near-infrared (0.7 to 4.0 microns). Rayleigh scattering by air molecules and absorption by ozone and water vapor are treated in the first band. In the 8 near-infrared bands, variable absorption by water vapor is included, and carbon dioxide absorption is calculated after a modified Sasamori et al. (1972) method. Pressure corrections and multiple reflections between clouds and the surface are treated, but not the radiative effects of aerosols.

Longwave calculations follow the simplified exchange method of Fels and Schwarzkopf (1975) and Schwarzkopf and Fels (1991) applied over seven spectral bands (with wavenumber boundaries at 0, 4.0 x 104, 5.6 x 104, 8.0 x 104, 9.9 x 104, 1.07 x 105, 1.20 x 105, and 2.20 x 105 m-1). Absorption by the vibrational and rotational lines of water vapor, carbon dioxide, and ozone, as well as continuum absorption by water vapor are treated, but some weak absorption bands of ozone and carbon dioxide are neglected. Carbon dioxide transmission coefficients are calculated for the actual temperature and pressure profile of each vertical column after the interpolation method of Fels and Schwarzkopf (1981). Longwave ozone and water vapor absorption (including temperature effects) are computed by a random-band model.

Cloud optical properties are prescribed. In the visible, cloud absorptivity is assumed to be zero; infrared absorptivity depends on cloud height, as do the visible and infrared cloud reflectivities. The absorptivity and reflectivity of clouds are also proportional to cloud amount. Longwave emissivity is prescribed as unity (blackbody emission) for all clouds. For purposes of calculating radiation and top-of-atmosphere cloud cover, clouds in different vertical layers are assumed to be randomly overlapped. Cf. McGregor et al. (1993) for further details. See also Cloud Formation.


After checking for supersaturation and attendant release of precipitation (see Precipitation) and performing a dry convective adjustment if needed, a modified Arakawa (1972) "soft" moist adjustment scheme predicts any subsequent precipitation release and the redistribution of moisture and momentum that may occur within sub-gridscale cumulus towers. These form when a layer (other than the lowest layer) is moist unstable with respect to at least one layer above, and when the relative humidity in the lowest unstable layer is > 75 percent. It is assumed that a constant convective mass flux effects a vertical redistribution of heat within the cumulus tower, such that the moist instability at each level (the difference between the moist static energy at cloud base and the saturation value at each level) decays with an e-folding time of one hour. (The heating at cloud base is assumed to be zero to ensure closure for the convective scheme.)

The attendant moisture redistribution and removal (see Precipitation) results in drying of the environment if the ambient relative humidity is at least 60 percent. The convective mass flux also transfers momentum upward through the cumulus tower, and downward via the surrounding large-scale descent. Shallow cumulus convection is parameterized by a modified Geleyn (1987) scheme that operates as an extension of the vertical diffusion of heat and momentum. The stability dependence of this diffusion is defined by a modified moist bulk Richardson number. Cf. McGregor et al. (1993) for further details.

Cloud Formation

Model description is to a large degree similar to AMIP database, however there are few important changes in cloud formulation (based on NCAR CCM2 for stratiform clouds. If convective activity occurs in the previous timestep (see Convection), convective cloud fraction is determined according to a formula by Slingo (1987).

Large-scale cloud amounts are determined from a modified form of the Rikus (1991) diagnostic, which is a quadratic function of the difference between the relative humidity of a layer and critical humidities that depend on cloud height (low, middle, and high cloud). Following Saito and Baba (1988) , maximum cloud fractions are also specified for each cloud type (0.70 for low cloud, 0.53 for middle cloud, and 0.50 for high cloud). Middle and high clouds are restricted to single layers, while low clouds can occupy two layers.

Stability-dependent low cloud associated with temperature inversions also may form if the relative humidity at cloud base is at least 60 percent. The cloud fraction is a function of the intensity of the temperature inversion, following Slingo (1987) . The overall fraction of low cloud is then set to the largest value predicted by either this stability-dependent diagnostic or by another operative mechanism (0.55 in the case of convective activity, or 0.70 for large-scale condensation).


Precipitation forms as a result of supersaturation and/or the moist convective adjustment process (see Convection). There is no subsequent evaporation of precipitation.

Planetary Boundary Layer

There are typically two atmospheric levels in the model PBL (whose top is not explicitly determined, however). In order to validate the simulation of surface atmospheric temperature against observations, a 2-meter (screen height) temperature is calculated from the bulk Richardson number determined for the lowest model layer (by applying the Monin-Obukhov assumption of constant momentum and heat fluxes in the surface layer). For unstable conditions, this requires an iterative solution. For purposes of computing surface fluxes, however, atmospheric winds, temperatures, and humidities at the first full atmospheric level above the surface (at sigma = 0.979) are used (see Surface Fluxes). Cf. McGregor et al. (1993) .


Orography from the 1 x 1-degree data of Gates and Nelson (1975) is transformed to spectral coefficients and truncated at the R21 resolution of the model. Orographic variance data (supplied by the United Kingdom Meteorological Office) are also used for the parameterization of gravity-wave drag (see Gravity-wave Drag).


For control: PMIP datasets will be 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

For 0fix and 6fix run, a sea ice model is used. It is a fully dynamic model described in Gordon and O'Farell 1997. Ocean currents are prescribed for fixed SST in both cases. Sea Ice Mass derived from sea ice thicknes(m) x concentration (fraction) x 900kg/m**3 (assumed density of sea ice).

SEAICE % is calculated from mean monthly concentration files.

Snow Cover

Precipitation falls to the surface as snow if the temperature of the second atmospheric vertical level above the surface (at sigma = 0.914) is below 0 degrees C. The latent heat is incorporated into the surface temperature prognostic for non-ocean surfaces (see Sea Ice and Land Surface Processes). Prognostic snow mass, with accumulation and melting over both land and sea ice, is included, but the allowable snow depth is limited to 4 m. Snow cover affects the heat capacity and conductivity of the land surface (see Land Surface Processes), and sublimation of snow contributes to the surface evaporative flux (see Surface Fluxes). Melting of snow, which contributes to soil moisture, occurs when the surface (top soil layer or sea ice) temperature is > 0 degrees C. Snow cover affects the albedo of the surface, but with less impact if the snow is melting (see Surface Characteristics). No fractional snow cover calculated in the model. Value of 3 mm of H2O was used as a cutoff for snow/nosnow in calculating this variable from snow thickness data. Snow mass. Maximum depth of snow is imposed in the model and is set to 400 mm H2O (snow calculated in the CSIRO GCM as water equivalent). If snow exceedes 400 mm it is turned to sea ice (over ice) and on land it is assumed to vanish (glaciers flow etc), however GCM take into account this into energy budget.

Surface Characteristics

There are 12 vegetations types taken from Dorman and Sellers (1989). J. Appl. Meteor., 28, 833-855. The roughness length over land is prescribed as a monthly values from Dorman and Sellers (1989) dataset at each point. However, the roughness length for other surface types is the same for heat and momentum fluxes: over ice, the roughness is a uniform 0.001 m, while over the ocean it is a function of surface wind stress, following Charnock (1955) .

Surface albedo and surface charactersitics are prescribed as a monthly values (and interpolated to a daily values) from Dorman and Sellers dataset.

Longwave emissivity is set to unity (blackbody emission) for all surface types. Cf. McGregor et al. (1993) for further details.

Surface Fluxes

Surface solar absorption is determined from surface albedos, and surface longwave emission from the Planck equation with prescribed emissivity of 1.0 (see Surface Characteristics). Following Monin-Obukhov similarity theory applicable to a constant-flux layer, surface turbulent eddy fluxes are expressed as bulk formulae. The requisite atmospheric surface winds, potential temperatures, and specific humidities are taken to be those at the first vertical level above the surface (at sigma = 0.979). The effective ground value of specific humidity needed for determination of the surface moisture flux from the bulk formula is obtained as a fraction alpha of the saturated humidity at the ground temperature, where alpha is a function of soil moisture (see Land Surface Processes) or is prescribed to be 1 over oceans and ice. The bulk drag and transfer coefficients are functions of roughness length (see Surface Characteristics) and vertical stability (bulk Richardson number) following Louis (1979) , with the same transfer coefficient used for the heat and moisture fluxes. Over the oceans, the neutral transfer coefficient for surface heat and moisture fluxes is a constant 8.5 x 10-4.

Land Surface Processes

There are 9 soil types is used from Zobler dataset. The soil uses following parameters; saturated moisture content, willing moisture content, saturation value of matrix potential, specific heat, density, albedo, emissivity, and field capacity. Vegetation uses following parameters; min observed conopy stomal resistance, fraction of vegetation, LAI, roughness lenght, albedo and emissivity. SVAT model is after Kowalczyk et al. The model used in PMIP experiment is exactly the same as in Gordon and O'Farell (1997), except dynamic ocean is not used.

E. A. Kowalczyk and J. R. Garratt and P. B. Krummel, 1991, A Soil-Canopy Scheme for Use in a Numerical Model of the Atmosphere - 1D Stand-Alone Model, CSIRO Division of Atmospheric Research, Technical Paper, 23

E. A. Kowalczyk and J. R. Garratt and P. B. Krummel, 1994, Implementation of a Soil-Canopy Scheme into the CSIRO GCMs, CSIRO Division of Atmospheric Research, Technical Paper, 32.

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