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

PMIP Documentation of ECHAM3

Model Max-Planck-Institut fuer Meteorologie: Model MPI ECHAM3 (T42 L19) 1994

PMIP Representative(s)

Dr. habil. Bjoern Grieger, formerly at the University of Bremen, now at the Max-Planck-Institute for Aeronomy, Max-Planck-Str. 2, D-37191 Katlenburg-Lindau, Germany, Phone: +49 5556 979 466, Fax: +49 2561 9595 6137, e-mail:


Stephan Lorenz, University of Bremen, Fachbereich 5, Postfach 330440, 28334 Bremen, Germany; Phone: +49 40 41173 175; Fax: +49 40 41173; e-mail:


Dr. Ulrike Wyputta, formerly at the University of Bremen, e-mail:



Model Designation

MPI ECHAM3 (T42 L19) 1994

Model Identification for PMIP


PMIP run(s)

0fix, 6fix, 21fix

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

Model Lineage

This is the third in a series of models developed at MPI that derive from an earlier version (cycle 17) of the operational forecast model of the European Centre for Medium-Range Weather Forecasts (ECMWF). The present MPI model retains some features (more in numerics and dynamics than in physics) of the ECMWF model. The model is identical to latest PMIP model exept for different initial conditions and the Earth's orbital parameters.

Model Documentation

main reference:

Deutsches Klimarechenzentrum (DKRZ) Modellbetreuungsgruppe, 1992: The ECHAM3 atmospheric general circulation model. DKRZ Tech. Report No. 6, ISSN 0940-9237, Deutsches Klimarechenzentrum, Hamburg, Germany, 184 pp.

Roeckner, E., K. Arpe, L. Bengtsson, S. Brinkop, L. Duemenil, M. Esch, E. Kirk, F. Lunkeit, M. Ponater, B. Rockel, R. Sausen, U. Schlese, S. Schubert, and M. Windelband, 1992: Simulation of the present-day climate with the ECHAM model: Impact of model physics and resolution. MPI Report No. 93, ISSN 0937-1060, Max-Planck-Institut fuer Meteorologie, Hamburg, Germany, 171 pp.

secondary reference:

Lorenz, S., B. Grieger, P. Helbig and K. Herterich 1996: Investigating the sensitivity of the Atmospheric General Circulation Model ECHAM 3 to paleoclimatic boundary conditions, Geologische Rundschau 85, 513-524

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. 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 a pressure of about 996 hPa.

Vertical Representation

Hybrid sigma-pressure coordinates after Simmons and Burridge (1981) and Simmons and Struefing (1981) .

Vertical Resolution

There are 19 irregularly spaced hybrid levels. The values of each level are:

0.0, 0.0, 0.0003389933, 0.0033571866, 0.0130700434, 0.0340771467, 0.0706498323, 0.1259166826, 0.2011954093, 0.2955196487, 0.4054091989, 0.5249322235, 0.6461079479, 0.7596983769, 0.8564375573, 0.9287469142, 0.9729851852, 0.9922814815, 1.0.

For a surface pressure of 1000 hPa, five levels are below 800 hPa and seven levels are above 200 hPa. The two lowest layers are pure sigma layers and the top two layers are at constant pressures.

Computer/Operating System

The PMIP simulation were run on a Cray 90 computer using a multi processors in the UNICOS environment.

Computational Performance

For the PMIP experiments, about 120 seconds of Cray 90 computation time per simulated day.


For the PMIP simulations, the model atmosphere is intialized from the ECMWF analysis for 1 January 1979, and the soil moisture and snow cover/depth from the ECMWF January climatology. Initial conditions of PMIP model for control simulation were taken from AMIP simulation.

Time Integration Scheme(s)

The semi-implicit time itegration scheme of Robert et al. (1972) and Robert (1981 , 1982 ) is applied with an Asselin (1972) frequency filter. The time step is 24 minutes for dynamics and physics, except for radiation which is calculated at 2-hour intervals.


Orography is smoothed (see Orography). Negative moisture values arising from truncation of the spherical harmonic basis functions are filled for purposes of the radiation calculations, but negative moisture values are tolerated in transport algorithms (advection, convection, and diffusion). Negative cloud water values are avoided by invoking a suitable condensation term (see Cloud Formation).

Sampling Frequency

For the PMIP simulations, the model history is written every 12 hours.

Dynamical/Physical Properties

Atmospheric Dynamics

Primitive-equation dynamics are expressed in terms of vorticity, divergence, temperature, and the logarithm of surface pressure, specific humidity, and cloud water (a model prognostic variable). Virtual temperature is also where applicable, for diagnostic variables.


Scale-selective fourth-order (Ñ4) horizontal diffusion is applied to all atmospheric prognostic variables for total spectral wave numbers >15 after the method of Laursen and Eliasen (1989) . To avoid unrealistic results near mountains, horizontal diffusion operates only on the deviation of the temperature field from that of a standard atmosphere. The horizontal diffusion coefficient varies for different variables and vertical levels, and is selected to make the slope of the atmospheric spectral kinetic energy approximate that observed.

Vertical diffusion operates above the planetary boundary layer (PBL) only in conditions of static instability. In the PBL, vertical diffusion of momentum, heat, cloud water, and moisture is proportional to the vertical gradient of the appropriate field. The vertically variable diffusion coefficient depends on stability (expressed as a bulk Richardson number) and the vertical shear of the u-wind, following standard mixing-length theory. See also Planetary Boundary Layer and Surface Fluxes.

Gravity-wave Drag

Drag associated with orographic gravity waves is simulated after the method of Palmer et al. (1986) , as modified by Miller et al. (1989) , using directionally dependent subgrid-scale orographic variances obtained from the U.S. Navy dataset (cf. Joseph 1980) . 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.

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 simulated. corrected anaytische loesung nicht reiehnentwicklung


The carbon dioxide concentrations are the PMIP-prescribed values of 345, 280 and 200ppm for 0fix, 6fix and 21fix, respectively. The ozone profile is determined from total ozone in a column (after data by London et al. 1976 ) and the height of maximum concentration (after data by Wilcox and Belmont 1977 ), and depends on pressure, latitude, longitude, and season. Radiative effects of water vapor and of three types of aerosol (oceanic, desert, urban) are also included (see Radiation).


The radiation schemes follow Hense et al. (1982) , Eickerling (1989) , and Rockel et al. (1991) . Radiative transfer equations are based on the two-stream approximation of Kerschgens et al. (1978) and Zdunkowski et al. (1980) , and are solved in 4 shortwave spectral intervals (with boundaries at 0.215, 0.685, 0.891, 1.273, and 3.580 microns) and in 6 longwave spectral intervals (with boundaries at 3.96, 7.98, 8.89, 10.15, 11.76, 20.10, and 100 microns).

Shortwave equations are expressed in terms of optical depth, single-scattering albedo, diffuse and direct backscattering parameters, and diffusivity factor; longwave equations (with scattering effects neglected) are expressed in terms of the Planck function, optical depth, and diffusivity factor. Absorbers determining shortwave and longwave optical depths include ozone, carbon dioxide, and water vapor (with continuum absorption included). Optical depths for diffuse shortwave and longwave absorption are calculated using coefficients derived from "exact" reference models. Effects of pressure broadening on longwave absorption are also treated.

In the shortwave, optical depths for Rayleigh scattering are determined from the molecular cross section for each gas, and those for absorption and Mie scattering by ocean, desert, and urban aerosols are from data of Shettle and Fenn (1975) . Single scattering albedo is derived from the optical depths for diffuse shortwave scattering/absorption. At a given vertical level, the total optical depth for direct shortwave radiation (which is dependent on solar zenith angle) is obtained from linear superposition of optical depths for scattering/absorption by gases and aerosols summed over all the levels above.

The shortwave optical depth for clouds is parameterized after the method of Stephens (1978) , and the single-scattering albedo and backscattering parameter follow Kerschgens et al. (1978) . Longwave emissivity is an exponential function of geometrical cloud thickness, prognostic cloud water content, and a mass absorption coefficient, following Stephens (1978). For purposes of the radiation calculations, clouds in contiguous vertical layers are treated as fully overlapped, and as randomly overlapped otherwise. Cf. DKRZ (1992) for further details. See also Cloud Formation.


The mass-flux convective scheme of Tiedtke (1989) accounts for shallow, midlevel and penetrative convection, as well as the effects of cumulus-scale downdrafts. Stratocumulus convection is parameterized as an extension of the model's vertical diffusion scheme (cf. Tiedtke et al. 1988 ). The closure assumption for midlevel/penetrative convection is that large-scale moisture convergence determines the bulk cloud mass flux; for shallow convection, the mass flux is maintained instead by moisture from surface evaporation. Entrainment and detrainment of mass in convective plumes occurs both through turbulent exchange and organized inflow and outflow. Momentum transport by convective circulations is also included, following Schneider and Lindzen (1976) .

Cloud Formation

The (stratiform) cloud-formation scheme is based on prognostic cloud water, following Sundqvist (1978) , Roeckner and Schlese (1985) , and Roeckner et al. (1991) . Subgrid-scale condensation and cloud formation in a fraction c of each grid box are governed by transport equations for water vapor and cloud water. The threshold relative humidity of the cloud-free fraction (1 - c) is prescribed as a function of height and stability: it is a minimum of 50 percent at tropopause level when penetrative convection occurs (cf. Xu and Krueger 1991 ); otherwise the threshold humidity decreases linearly from 99 percent at the surface to 85 percent at the top of the PBL (see Planetary Boundary Layer), above which it remains constant.

Cloud droplets grow by condensation if the grid-mean relative humidity exceeds this threshold relative humidity or if the relative humidity in the c fraction of the grid box exceeds 100 percent. The condensation rate depends on the moisture convergence into the grid box, with a fraction c of the convergence producing more cloud, and the fraction (1 - c) increasing the relative humidity of the cloud-free part of the grid box.

The cloud water in the cloud-free fraction (1 - c) that increases (decreases) as a result of advective/diffusive transports into the grid cell or through numerical effects (e.g., spectral truncation) is assumed to evaporate (condense) instantaneously; explicit filling of negative cloud water values is therefore unnecessary (see Smoothing/Filling). Convective cloud water is not advected (rather, there is instantaneous precipitation and/or evaporation to the environment--see Precipitation). Cf. DKRZ (1992) for further details. See also Radiation for treatment of cloud-radiative interactions.


For convective precipitation, freezing and melting processes are not considered. No cloud water is stored in convective cloud (see Cloud Formation); once detrained, it evaporates instantaneously, with any portion not moistening the environment falling out as precipitation. Conversion from convective cloud droplets to rain/snow is proportional to the product of cloud water content and upward convective mass flux at cloud base, weighted by an empirical function of height (cf. Yanai et al. 1973) . Convective snow forms if the temperature of the cloud layer is <0 degrees C. Following Kessler (1969) , evaporation of convective precipitation is assumed to be proportional to the saturation deficit and the rainfall intensity.

For stratiform mixed-phase precipitation formation (i.e., in a temperature range from 0 to -40 degrees C), the ice and liquid phases are treated independently. Growth of cloud droplets (see Cloud Formation) to precipitating raindrops occurs by autoconversion, following the exponential relationship of Sundqvist (1978) , and by collisions with larger drops, following the parameterization of Smith (1990) . Partitioning of cloud liquid vs ice is according to the temperature-dependent relation of Matveev (1984) , and the loss of ice crystals by sedimentation follows Heymsfield (1977) . Evaporation of stratiform precipitation in a layer below cloud is proportional to the saturation deficit, but cannot exceed the precipitation flux at the layer top. Stratiform snow forms if the cloud layer temperature snow forms if the cloud layer temperature is <0 degrees C. Falling convective and stratiform snow melts if the temperature of a layer is >2 degrees C. Cf. DKRZ (1992) for further details. See also Snow Cover.

Planetary Boundary Layer

The PBL is typically represented by the first 5 vertical levels above the surface (see Vertical Resolution). The PBL top (usually at elevations <2000 m) is determined as the greater of the height predicted from Ekman theory vs a convective height that depends on the dry static energy in the vertical. Heat, moisture, cloud water, and momentum fluxes follow Monin-Obukhov similarity theory at the surface and standard mixing-length theory above the surface. See also Diffusion and Surface Fluxes.


Orography is obtained from the U.S. Navy dataset with resolution of 10 minutes arc on a latitude/longitude grid (cf. Joseph 1980) for 0fix and 6fix runs. It is modified after PMIP recommendation for 21 fix run, using the difference of the orographic difference between present and 21ka. These data are smoothed using a Gaussian filter with radius of influence 50 km. The resulting heights are transformed into spectral space and truncated at T42 resolution.


For control: The SSTs (sea surfacce temperatures) and the sea ice distribution are taken from the COLA/CAC AMIP data set, which was prepared by the Center for Ocean-Land-Atmosphere Interactions (COLA) and the Climate Analysis center within the Atmospheric Modeling Intercomparison Project (AMIP), (Reynolds, 1988).

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: a global reconstruction of the surface conditions was supplied by the CLIMAP project (1981). It includes the sea ice distribution and SSTs for February and August. The annual cycle of SSTs was approximated by a sine function, taking the values for February and August like extrema.

Sea Ice

For 0fix and 6fix runs, monthly AMIP sea ice extents are prescribed. The temperature of the upper 0.10 meter of sea ice is computed from energy fluxes at the surface (see Surface Fluxes) and from the ocean below. The ocean heat flux depends on the ice thickness and the difference between the temperature of the underlying ocean and that of the ice. Snow does not accumulate on sea ice. Cf. DKRZ (1992) for further details. See also Snow Cover.

For 21fix run, sea ice edge of CLIMAP data is used for February and August. The area covered by the glacial ice sheets is also taken from CLIMAP, but the the thickness of the ice sheets is altered according to Tushingham and Peltier (1991), which implies a thickness reduction of Laurentide ice sheet by about one third compared with the CLIMAP thickness and a corresponding increase of the thickness of the other ice sheets.

Snow Cover

Snow falls to the surface if air temperatures at levels below where it forms (see Precipitation) are <2 degrees C. Snow accumulates only on land to a depth that is determined prognostically from a budget equation and melts if the temperatures of the snow and of the uppermost soil layer exceed 0 degrees C (see Land Surface Processes). Snow cover affects surface albedo (see Surface Characteristics) as well as heat transfer and (after melting) soil moisture. The fractional area of a grid box that is snow covered is given by the ratio of the water-equivalent snow depth to a critical value of 0.015 m, with complete coverage if this depth is exceeded. Sublimation of snow is calculated as part of the surface evaporative flux (see Surface Fluxes).

Surface Characteristics

The fraction of each grid box covered by vegetation is determined from Matthews (1983) 1 x 1-degree data, but is reduced when soil moisture available in the root zone becomes less than 40 percent of field capacity (0.20 m)--see Land Surface Processes.

The surface roughness is calculated in differnet ways over land and over sea, where the land/sea distribution is taken from a US navy data set. The surface roughness length is prescribed as a uniform 1 x 10-3 meter over sea ice; it is computed prognostically over open ocean from the surface wind stress by the method of Charnock (1955) , but is constrained to be a minimum of 1.5 x 10-5 m. Over sea ice, the roughness length is assumed to be constant. Over land, the roughness length is a function of the orography, where the vegetation is taken from Wilson and Henderson-Sellers (1985). (see Orography) There exist 3 types of vegetation: forest, no vegetation and land ice.

The annual-mean surface albedo is obtained from satellite data of Geleyn and Preuss (1983) . Over land, this background albedo is altered by snow cover as a linear function of the ratio of the water-equivalent snow depth to a critical value (0.01 m). Albedos of snow (range 0.30 to 0.80), sea ice (range 0.50 to 0.75), and continental ice (range 0.6 to 0.8) vary as a function of surface temperature and forested area, as given by Robock (1980) and Kukla and Robinson (1980) . The albedo of ocean is a constant 0.065 for diffuse radiation, while that for the direct beam depends on solar zenith angle, but never exceeds 0.15.

Longwave emissivity is prescribed as 0.996 for all surfaces. Cf. DKRZ 1992 for further details.

Surface Fluxes

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

Surface turbulent eddy fluxes of momentum, dry static energy (sensible heat), cloud water, and moisture are simulated as stability-dependent bulk formulae, following Monin-Obukhov similarity theory. The required near-surface values of wind, temperature, cloud water, and humidity are taken to be those at the lowest atmospheric level (sigma = 0.996). (At the surface, cloud water is assumed to be zero.) Surface drag and transfer coefficients in the bulk formulae are functions of stability and roughness length (see Surface Characteristics), following Louis (1979) and Louis et al. (1981) , but with modifications by Miller et al. (1992) for calm conditions over the oceans. The stability criterion is the moist bulk Richardson number, which includes the impact of cloud processes on buoyancy (cf. Brinkop 1992) .

The surface moisture flux depends on the surface specific humidity; over ocean, snow, ice, and wet vegetation fractions of each grid box, this is taken as the saturated humidity at the surface temperature and pressure (i.e., potential evaporation is assumed). Over the bare soil fraction, the surface specific humidity is the product of relative humidity (that is a function of soil moisture--see Land Surface Processes) and the saturated specific humidity. For a dry vegetation canopy, the potential evaporation is reduced by an evapotranspiration efficiency factor beta that is the inverse sum of aerodynamic resistance and stomatal resistance; the latter depends on radiation stress, canopy moisture, and soil moisture stress in the stress in the vegetation root zone (cf. Sellers et al. 1986 , Blondin 1989 , and Blondin and Boettger 1987 ).

Land Surface Processes

Soil temperature is determined after Warrilow et al. (1986) from the heat conduction in 5 layers (proceeding downward, layer thicknesses are 0.065, 0.254, 0.913, 2.902, and 5.70 m), with net surface heat fluxes (see Surface Fluxes) as the upper boundary condition and zero heat flux as the lower boundary condition at 10 m depth.

Snow pack temperature is also computed from the soil heat equation using heat diffusivity/capacity for ice in regions of permanent continental ice, and for bare soil where water-equivalent snow depth is <0.025 m. For snow of greater depth, the temperature of the middle of the snow pack is solved from an auxiliary heat conduction equation (cf. Bauer et al. 1985 ). The temperature at the upper surface is determined by extrapolation, but it is constrained not to exceed the snowmelt temperature of 0 degrees C.

There are separate prognostic moisture budgets for snow, vegetation canopy, and soil reservoirs. Snow cover is augmented by snowfall and is depleted by sublimation and melting (see Snow Cover). Snow melts (augmenting soil moisture) if the temperatures of the snow pack and of the uppermost soil layer exceed 0 degrees C. The canopy intercepts precipitation and snow (proportional to the vegetated fraction of a grid box), which is then subject to immediate evaporation or melting.

Soil moisture is represented as a single-layer "bucket" model (cf. Manabe 1969 ) with field capacity 0.20 m that is modified to account for vegetative and orographic effects. Direct evaporation of soil moisture from bare soil and from the wet vegetation canopy, as well as evapotranspiration via root uptake, are modeled (see Surface Fluxes). Surface runoff includes effects of subgrid-scale variations of field capacity related to the orographic variance (see Orography); in addition, wherever the soil is frozen, moisture contributes to surface runoff instead of soil moisture. Deep runoff due to drainage processes also occurs independently of infiltration if the soil moisture is between 5 and 9 percent of field capacity (slow drainage), or is larger than 90 percent of field capacity (fast drainage). Cf. Duemenil and Todini (1992) for further details.

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