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

PMIP Documentation for UGAMP

The UK Universities' Global Atmospheric Modelling Programme: Model UGAMP UGCM Version 2 (T42 L19) 1994


PMIP Representative(s)

Dr. Paul Valdes , Department of meteorology, University of Reading, Reading RG6 6BB, U.K.; Phone: 44 0118 9316517; Fax: 44 0118 9318905; e-mail:


Dr. Buwen Dong, Department of meteorology, University of Reading, Reading RG6 6BB, U.K.; Phone: 44 0118 9875123 ext:4384; Fax: 44 0118 9318316; e-mail: ;

World Wide Web URL:

Model Designation

UGAMP UGCM Version 2 (T42 L19) 1994

Model Identification for PMIP


PMIP run(s)

0fix, 6fix, 21fix, 0cal, 21cal

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

Model Lineage

The UGAMP model is based on the ECMWF (cycle 27) model (cf. Tiedtke et al. 1988 and Simmons et al. 1989 ), but with modifications principally in the treatment of radiation, convection, surface fluxes, vertical advection, and lateral and vertical dissipation (Slingo et al. 1994). The bottom boundary conditions over land was modified and coupled to a slab sea ice model (Valdes and Hall 1994) for paleoclimate modeling. The model atmosphere, soil moisture, and snow cover/depth were initiated from the ECMWF operational analysis for 12Z on 15 January 1987. Over the permanent ice caps, the snow depth was initiated at 10 meter water equivalence.

Model Documentation

main reference(s)

The 21fix and 21cal simulations were documented by :

Dong B. W. and P. J. Valdes, 1998: Simulations of the Last Glacial Maximum climates using a general circulation model: Prescribed versus computed Sea Surface temperatures. Climate dynamics, 14, 571-591.

and 6fix was documeted by :

Hall, N.M.J. and P. J. Valdes, 1997: A GCM simulation of the climate 6000 years ago. J. Clim., 10, 3-17.

secondary reference(s)

Valdes P. J., and N. M. J. Hall, 1994: Mid-latitude depressions during the ice age. NATO ASI Volume on Long Term Climatic Variations-Data and Modelling. J. C. Duplessy, Ed., 511-531.

Dong, B., and P. J. Valdes, 1995:Sensitivity studies of northern hemisphere glaciation using an atmospheric general circulation model. J. Climate, 8, 2471-2496.

Simmons, A. J., D. M. Burridge, M. Jarraud, C. Girard, and W. Wergen, 1989: The ECMWF medium-range prediction models: Development of the numerical formulations and the impact of increased resolution. Meteor. Atmos. Phys., 40, 28-60.

Subsequent modifications are described by Slingo et al. (1994) and references therein:

Slingo, J. M., M. Blackburn, A. Betts, R. Brugge, K. Hodges, B. Hoskins, M. Miller, L. Steenman-Clark, and J. Thuburn, 1994: Mean climate and transience in the tropics of the UGAMP GCM: Sensitivity to convection parameterization. Quart. J. Roy. Meteor. Soc., 120, 881-922

Documentation for the ECMWF(cycle 27) predecessor model is provided by Tiedtke et al. (1988) : Tiedtke, M., W. A. Heckley and J. M. Slingo, 1988: Tropical forecasting at ECMWF: The influence of physical parameterization on the mean structure of forecasts and analyses. Quart. J. Roy. Meteor. Soc., 114, 639-664

Numerical/Computational Properties

Horizontal Representation

Spectral (spherical harmonic basis functions) with transformation to a Gaussian grid for calculation of nonlinear quantities and some physics.

Horizontal Resolution

Spectral triangular 42 (T42), roughly equivalent to 2.8 x 2.8 degrees latitude-longitude. The transform grid is sufficient to prevent aliasing of quadratic quantities, with 128 equispaced longitudes and 64 Gaussian latitudes. The full radiative calculations are performed on a reduced longitudinal grid, retaining only the first 16 Fourier modes (see Radiation).

dim_longitude*dim_latitude: 128*64

Vertical Domain

Surface to 10 hPa; for a surface pressure of 1000 hPa, the lowest atmospheric level is at about 996 hPa.

Vertical Representation

Hybrid sigma-pressure coordinates after Simmons and Burridge (1981) and Simmons and Strüfing (1981) . To avoid oscillations in the profile of an advected quantity with rapidly changing gradient, vertical advection is treated by the Total Variation Diminishing (TVD) scheme of Thuburn (1993) .

Vertical Resolution

There are 19 irregularly spaced hybrid pressure-sigma levels. 5 levels are below 800hPa and 7 levels are above 200hPa. The transform grid is with 128 equispaced longitudes and 64 Gaussian latitudes.

The UGAMP GCM uses a hybrid coordinate in the vertical. The half-level values are given by


where P(0) is a constant pressure. This coordinate is identical to the usual

sigma when A(K+1/2)=0, and in general equals sigma when p(s)=p(0)*eta=p/p(0)

at levels where coordinate surfaces are surfaces of constant pressure. Values

of eta in between half-levels are given by


for p in the range between p(k-1/2) and p(k+1/2).

The value for p(0) is 101325 Pa. A(k+1/2) and B(k+1/2) are as following

k A(k+1/2) B(k+1/2)

0 0.000000 0.0000000000

1 2000.000000 0.0000000000

2 4000.000000 0.0000000000

3 6046.110595 0.0003389933

4 8267.927560 0.0033571866

5 10609.513232 0.0130700434

6 12851.100169 0.0340771467

7 14698.498086 0.0706498323

8 15861.125180 0.1259166826

9 16116.236610 0.2011954093

10 15356.924115 0.2955196487

11 13621.460403 0.4054091989

12 11101.561987 0.5249322253

13 8127.144155 0.6461079479

14 5125.141747 0.7596983769

15 2549.969411 0.8564375573

16 783.195032 0.9287469142

17 0.000000 0.9729851852

18 0.000000 0.9922814815

19 0.000000 1.000000000

Computer/Operating System

The PMIP simulation was run on a Cray ymp8 computer using a single processor in a UNICOS environment.

Computational Performance

For the PMIP experiment, about 8 minutes Cray ymp8 computer time per simulated day (including data-archiving and storage time).


For the control experiment, the model atmosphere, soil moisture, and snow cover/depth were initialized from the ECMWF operational analysis for 12Z on 15 January 1987.

Time Integration Scheme(s)

The time integration is by a semi-implicit Hoskins and Simmons (1975) scheme with an Asselin (1972) time filter. Advection of vorticity and moisture by a zonally symmetric flow is also treated implicitly. The time step is 30 minutes for dynamics and physics, except for full radiation/cloud calculations once every 3 hours (on a reduced grid, at every fourth point in longitude only--see Radiation). To ensure mass conservation, the global mean value of the logarithm of surface pressure is rescaled at each time step (but with mass sources/sinks associated with evaporation/precipitation neglected).


Orography is smoothed (see Orography). Negative values of atmospheric specific humidity (due to numerical truncation errors in the discretized moisture equation) are filled by borrowing moisture from successive vertical levels below until all specific humidity values in the column are nonnegative. Any moisture which must be borrowed from the surface does not affect the hydrological budget there.

Sampling Frequency

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

Dynamical/Physical Properties

Atmospheric Dynamics

Primitive-equation dynamics are expressed in terms of vorticity, divergence, temperature, the logarithm of surface pressure, and specific humidity. Variations of the gas constant and specific heat capacity with water vapor content are also included.


Sixth-order (Ñ6) hyperdiffusion is applied in spectral space to vorticity, divergence, temperature, and moisture on the hybrid coordinate surfaces (see Vertical Representation). A correction is also applied to the temperature term to approximate dissipation on constant pressure surfaces. The diffusion time scale is 4 hours at the horizontal truncation limit (see Horizontal Resolution), but this is successively halved on the top four model levels, beginning at approximately 73 hPa.

Second-order vertical diffusion is applied below a hybrid model coordinate level of 0.650 to parameterize the PBL (see Planetary Boundary Layer). In addition, the TVD vertical advection scheme (see Vertical Representation) includes some dissipation of kinetic energy where sharp changes in gradient are encountered.

Gravity-wave Drag

Momentum transports associated with gravity waves are simulated by the method of Palmer et al. (1986) , using directionally dependent subgrid-scale orographic variances obtained from the U.S. Navy dataset (cf. Joseph 1980 and see Orography). Surface stress due to gravity waves excited by stably stratified flow over irregular terrain is calculated from linear theory and dimensional considerations. Gravity-wave stress is a function of atmospheric density, low-level wind, and the Brunt-Vaisalla frequency. The vertical structure of the momentum flux induced by gravity waves is calculated from a local wave Richardson number, which describes the onset of turbulence due to convective instability and the turbulent breakdown approaching a critical level. See also Orography.

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 correct annual calendar is used, including Leap Years 1980, 1984, and 1988.)


Carbon dioxide concentration is prescribed value of 345, 280 and 200 ppm for control, 6ka and 21ka, respectively. The specified ozone profile depends on pressure, total ozone in a column, the height of maximum concentration, latitude, longitude, and season. Total ozone is obtained from London et al. (1976) data, and the altitude of maximum concentration from Wilcox and Belmont (1977) . Mie radiative parameters of five types of aerosol are provided for (concentration depending only on height) from WMO-ICSU (1984) data. Radiative effects of water vapor, carbon monoxide, methane, nitrous oxide, and oxygen are also included (see Radiation).


Atmospheric radiation is simulated after Morcrette (1989 , 1990 , 1991 ). Absorption by water vapor, ozone, carbon monoxide, carbon dioxide, methane, nitrous oxide, and oxygen is accounted for, with shortwave/longwave absorption coefficients calculated from line parameters of Rothman et al. (1983) .

For clear-sky conditions, shortwave radiation is modeled by a two-stream formulation in two spectral wavelength intervals (0.25 to 0.68 micron and 0.68 to 4.0 microns), using a photon path distribution method to separate the contributions of scattering and absorption processes to radiative transfer. Rayleigh scattering and Mie scattering/absorption by five aerosol types (see Chemistry) are treated by a delta-Eddington approximation.

The clear-sky longwave scheme employs a broad-band flux emissivity method in six spectral intervals between wavenumbers 0 and 2.6 x 105 m-1, with continuum absorption by water vapor included between wavenumbers 3.5 x 104 to 1.25 x 105 m-1. The temperature/pressure dependence of longwave gaseous absorption follows Morcrette et al. (1986) . Aerosol absorption is also modeled by an emissivity formulation.

Shortwave scattering and absorption by cloud droplets are treated by a delta-Eddington approximation; radiative parameters include optical thickness, single-scattering albedo linked to cloud liquid water path, and prescribed asymmetry factor. Cloud types are distinguished by also defining shortwave optical thickness as a function of effective droplet radius. Clouds are treated as graybodies in the longwave, with emissivity depending on cloud liquid water path after Stephens (1978) . Longwave scattering by cloud droplets is neglected, and droplet absorption is modeled by an emissivity formulation in terms of the cloud liquid water path. For purposes of the radiation calculations, clouds of different types are treated as randomly overlapped in the vertical; convective cloud and the same type of nonconvective cloud in adjacent layers are treated as fully overlapped.

The full radiation calculations are performed every 3 hours on a reduced horizontal grid (every fourth point in longitude only), but with effective transmissivities and emissivities returned on the T42 Gaussian grid (see Horizontal Resolution). For intermediate time steps, the effective transmissivities are scaled by the instantaneous incoming solar radiation to represent correctly the diurnal cycle; the effective emissivities are scaled by the instantaneous Planck function to treat temperature variations. However, the influence of clouds remains fixed between full-radiation steps. See also Cloud Formation.


Convection follows the scheme of Betts and Miller (1993) , and consists of a relaxed convective adjustment towards calculated temperature and humidity reference profiles based on observations. The relaxation times are 2 hours for precipitating deep convection and 4 hours for nonprecipitating shallow convection, regarded as mutually exclusive processes. The convection is treated as shallow if the cloud top, defined by the level of nonbuoyancy, is below about 725 hPa for land or 810 hPa for ocean (for a surface pressure of 1000 hPa), or if there is insufficient moisture for precipitation to form with deep convection. Otherwise, the deep-convection scheme (including the possibility of midlevel convection with a cloud base above 725 hPa) is operative.

The temperature and humidity reference profiles for deep convection are based on relevant observational data (cf. Betts 1986 ). The temperature reference profile is a lapse rate that is slightly unstable with respect to the wet virtual adiabat below the freezing level, and that returns at cloud top to the moist adiabat of the cloud base. For energy conservation, this reference profile is corrected (with a second iteration) in order to remove the vertically integrated difference between the total moist enthalpy of the environment and that of the reference profile. The humidity reference profile is derived from the temperature reference by linearly interpolating between the humidities for specified values of subsaturation pressure deficit at cloud base, freezing level, and cloud top. Below the cloud base, cooling/drying by convective downdrafts is parameterized by specifying reference profiles for air parcels originating near 850 hPa that descend at constant subsaturation and equivalent potential temperature.

Nonprecipitating shallow convection is parameterized as a mixing of enthalpy and moisture of air below cloud base with air at and just above the capping inversion top. The reference profile is a mixing line structure joining the conserved saturation pressure and potential temperature points of all mixtures of the two sources of air (cf. Betts 1983 , 1986 ). Reference temperature and humidity profiles are computed after specifying a partial degree of mixing within the cloud, and mixing that is a function of the inversion strength at cloud top. Cf. Betts and Miller (1993) and Slingo et al. (1994) for further details. See also Cloud Formation and Precipitation.

Cloud Formation

Cloud formation is simulated following the diagnostic method of Slingo (1987) . Clouds are of three types: shallow and deep convective cloud (see Convection); stratiform cloud associated with fronts/tropical disturbances that forms in low, middle, or high vertical layers; and low cloud associated with temperature inversions.

The fraction of shallow convective cloud (typically about 0.30) is related to the moisture tendencies within the cloud layer (cf. Betts and Miller 1993 ). The fraction of deep convective cloud (ranging between 0.20 to 0.80) is determined from the scaled convective precipitation rate (see Precipitation). If deep convective cloud forms above 400 hPa and the fractional area is > 0.4, anvil cirrus and shallow convective cloud also form.

Stratiform cloud is present only when the local relative humidity is > 80 percent, the amount being a quadratic function of this humidity excess. Low stratiform cloud is absent in regions of grid-scale subsidence, and the amount of low and middle stratiform cloud is reduced in dry downdrafts around subgrid-scale convective clouds. Low cloud forms below a temperature inversion if the relative humidity is > 60 percent, the cloud amount depending on this humidity excess and the inversion strength. See also Radiation for treatment of cloud-radiative interactions.


Precipitation is obtained from deep convection as part of the relaxed adjustment to the reference temperature and humidity profiles (see Convection). Subsequent evaporation of this precipitation is implicitly treated through inclusion of effects of convective downdrafts in the lowest three atmospheric layers. There is additional evaporation below elevated convective cloud bases that are situated above these downdraft layers.

In the absence of convective adjustment, precipitation also results from gridscale condensation when the local specific humidity exceeds the saturated value at the ambient temperature and pressure; the amount of precipitate depends on the new equilibrium specific humidity resulting from the accompanying latent heat release. Before falling to the surface, grid-scale precipitation must saturate all layers below the condensation level by evaporation. Melting of falling snow (see Snow Cover) occurs for air temperatures > +2 degrees C.

Planetary Boundary Layer

Vertical diffusion of momentum, heat, and moisture (proportional, respectively, to the vertical gradients of the wind, the dry static energy, and the specific humidity) is operative only below a hybrid-coordinate vertical level of 0.650 (about 650 hPa for a surface pressure of 1000 hPa). The vertically varying diffusion coefficient depends on stability (bulk Richardson number) and the vertical shear of the wind, following standard mixing-length theory (cf. Louis 1979 and Louis et al. 1981 ). See also Diffusion, Surface Characteristics, and Surface Fluxes.


For 0fix, 0cal and 6fix, orography is obtained from a U.S. Navy dataset (cf. Joseph 1980 ) with resolution of 10 minutes arc on a latitude/longitude grid. The mean terrain heights are then calculated for a T106 Gaussian grid, and the square root of the corresponding subgridscale orographic variance is added. The resulting "envelope orography" (cf. Wallace et al. 1983 ) is smoothed by applying a Gaussian filter with a 50 km radius of influence (cf. Brankovic and Van Maanen 1985 ). This filtered orography is then spectrally fitted and truncated at the T42 resolution of the model. See also Gravity-wave Drag.

In 21k simulations, the anomalies of 21k and 0k orography based on Peltier's reconstruction was added to UGAMP 0k orography. The land-sea mask is also based on Peltier's reconstruction. The ice sheet reconstruction for 21k is from Peltier et al. (1994).


For control: SSTs and sea ice extent are specified from the Alexander and Mobley (1976) dataset.

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

For 21 fix: SSTs and sea-ice: The change in SSTs (LGM minus present-day) given by CLIMAP (1981), available at NGDC rather than the LGM absolute values in order to avoid differences due to differences in present day climatologies, was used. To obtain seasonally varying SSTs and sea ice edge from data for February and August, a simple sinusoidal variations, with extrems in February and August, is used.

For computed SSTs experiments, the model was coupled to a mixed layer ocean model (Dong and Valdes, 1995), with prescribed, seasonally varying ocean heat flux transport. The ocean heat transport is diagnosed seperately over ocean and sea ice from the contral simulation.

For 21 cal: The 21 kyr BP experiments with computed SSTs and sea ice should be performed in the same way as the CO2 experiments. Typically this means that a coupled atmosphere-mixed layer ocean model will be used with prescribed present day ocean heat flux.

Sea Ice

For 0fix and 6fix, sea ice extents are prescribed monthly by AMIP, but ice surface temperatures are calculated using a simple sea ice model. This is a purely thermodynamic model in which heat is stored and diffused vertically across a slab of ice which is 2 meters thick. (Points with surface temperatures < -2 degrees C that are not on land are identified as sea ice; the masking procedure is described by Brugge 1993). Snow does not accumulate on sea ice.

For 21fix, sea ice edge of CLIMAP data is used for February and August.

Snow Cover

Grid-scale precipitation falls as snow if the temperature of the cloud layer is below 0 degrees C and that of intervening layers is below +2 degrees C (thereby inhibiting the melting of falling snow--see Precipitation). Snow depth (in meters of equivalent liquid water) is determined prognostically from a budget equation, but with accumulation only on land. The fractional area of a grid box covered by snow is given by the ratio of the snow depth to a critical depth (0.015 m), or is set to unity if the depth exceeds the critical value. Over the permanent ice caps, the snow depth was initiated at 10 meter water equivalence. Sublimation of snow is calculated as part of the surface evaporative flux (see Surface Fluxes), and snowmelt (occurring for ground temperatures > 0 degrees C) contributes to soil moisture (see Land Surface Processes). Snow cover also alters the surface albedo (see Surface Characteristics).

Surface Characteristics

The land surface is modeled as bare or with snow cover. Vegetation is not explicitly specified, but is accounted for in the prescribed surface properties described below.

Roughness length is prescribed as 1.0 x 10-3 m over sea ice. Over open ocean the roughness is computed from the surface wind stress following Charnock (1955) , but it is constrained to be at least 1.5 x 10-5 m. The roughness length over land is prescribed as a blended function of local orographic variance (Tibaldi and Geleyn 1981 ), vegetation (Baumgartner et al. 1977 ), and urbanization (from the U.S. Navy data set described by Joseph 1980 ) that is interpolated to the model grid; the logarithm of local roughness length is also smoothed by the same Gaussian filter used for orography (see Orography).

Annual means of satellite-observed surface albedo (range 0.07 to 0.80) from data of Preuss and Geleyn (1980) and Geleyn and Preuss (1983) are interpolated to the model grid and smoothed by the same Gaussian filter as used for orography (see Orography). Snow cover alters this background albedo, with a limiting value of 0.80 for snow depths > 0.01 m equivalent water. Sea ice albedo is prescribed as 0.55, and ocean albedo as 0.07. All albedos are also functions of solar zenith angle.

Longwave emissivity is prescribed as 0.996 for all surfaces. See also Sea Ice, Snow Cover, Surface Fluxes, and Land Surface Processes.

Surface Fluxes

Surface solar absorption is determined from surface albedos, and longwave emission from the Planck equation with prescribed constant surface emissivity (see Surface Characteristics).

Surface turbulent eddy fluxes are simulated as stability-dependent diffusive processes, following Monin-Obukhov similarity theory. Fluxes of momentum/heat/moisture are calculated from bulk formulae that include the product of a drag/transfer coefficient, the low-level wind speed, and the vertical difference between winds/dry static energy/specific humidity at the surface and their values at the lowest atmospheric level (996 hPa for a surface pressure of 1000 hPa). The low-level wind speed includes an imposed minimum of 3 m/s and an additional 3 m/s (added quadratically) in the presence of convection. (The former quantity increases surface fluxes in the limit of low wind speed, while the latter accounts for subgrid-scale convective circulations--cf. Slingo et al. 1994 .) The surface drag/exchange coefficients are functions of stability (bulk Richardson number) and roughness length (see Surface Characteristics) following the formulation of Louis (1979 ) and Louis et al. (1981) . The same transfer coefficient is used for the surface heat and moisture fluxes.

The surface moisture flux is also equivalent to the potential evaporation from a saturated surface multiplied by an evapotranspiration efficiency factor beta (cf. Budyko 1974 ). The factor beta is specified as unity over oceans and regions of dew formation (where the lowest atmospheric level is supersaturated); otherwise, beta varies with the snow cover and soil moisture content (see Snow Cover and Land Surface Processes).

Land Surface Processes

Vegetation is not explicitly specified, but is accounted for in the prescribed surface properties, such as Roughness length, bare land albedo. Soil temperature and moisture are calculated using a three layer diffusive model. It consists of a surface layer 0.15 m thick, middle layer 0.75 m thick, and deep layer 2.5 m thick. At the bottom of the soil model a no-flux boundary condition is used. Soil temperature is determined by simulating heat diffusion with an upper boundary condition specified by the net balance of surface energy fluxes (see Surface Fluxes). Soil heat capacity and diffusivity are prescribed constants: the density weighted heat capacity is 2.4 x 106 J/(m3 K) and heat diffusivity is 7.5 x 10-7 m2/s).

Soil moisture also obeys a diffusion equation (with diffusivity one-seventh that of the heat diffusivity). The upper boundary condition is specified from the combined rainfall and snowmelt, and from surface evaporation that is reduced by the presence of (fractional) snow cover. Runoff occurs if the soil moisture exceeds the layer capacity (scaled according to thickness: 0.02 m for the surface layer and 0.10 m for the middle and 0.333 m for the bottom layer). The evapotranspiration efficiency factor beta (see Surface Fluxes) is a composite of values determined for the snow-covered and bare-land fractions of a grid box. For snow-covered surfaces (see Snow Cover), beta is unity. Over bare land, beta is the ratio of the surface layer moisture to a prescribed fraction (0.75) of field capacity, but is constrained to be at most unity. There is also a temperature-dependent correction to account for limitation of evaporation due to lack of shortwave radiation.

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