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

PMIP Documentation for GFDL

Geophysical Fluid Dynamics Laboratory : Model GFDL CDG (R30 L20) 1997

PMIP Representative(s)

Anthony J. Broccoli, Geophysical Fluid Dynamics Laboratory/NOAA, Princeton University, P.O. Box 308, Princeton, New Jersey 08540; Phone: +1-609-452-6671; Fax: +1-609-987-5063; e-mail:

World Wide Web URL:

Model Designation

GFDL CDG (R30 L20) 1997

Model Identification for PMIP


PMIP run(s)

0fix, 6fix, 0cal, 21cal

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

Model Lineage

One of several versions of global spectral models in use at the Geophysical Fluid Dynamics Laboratory, ithis version of the GFDL model is applied primarily to climate studies. In its treatment of atmospheric dynamics, the model is similar to the GFDL Dynamical Extended-Range Forecasting (DERF) model (cf. Gordon and Stern 1982 and Miyakoda and Sirutis 1986 ), but it displays a number of differences in numerical features (e.g., use of rhomboidal rather than triangular spectral truncation and a somewhat higher vertical resolution), as well as in some model parameterizations (e.g., radiation, vertical diffusion, land surface submodel, sea ice submodel).

Model Documentation

Key documentation of dynamical features is given by Gordon and Stern (1982) and Manabe and Hahn (1981) . Physical parameterizations are described by Manabe (1969) , Holloway and Manabe (1971) , Manabe and Holloway (1975) , Manabe et al. (1965) , Broccoli and Manabe (1992) , Wetherald and Manabe (1988) , and Wetherald et al. (1991) .

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 rhomboidal 30 (R30), roughly equivalent to 2.25 x 3.75 degrees latitude-longitude.

dim_longitude*dim_latitude : 96*80

Vertical Domain

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

Vertical Representation

Finite-difference sigma coordinates.

Vertical Resolution

!The model uses asigma-coordinate in the vertical with 20 levels:0.0091, 0.0308, 0.0600, 0.0976, 0.1438, 0.1982, 0.2600, 0.3281, 0.4008, 0.4765, 0.5533, 0.6291, 0.7020, 0.7701, 0.8315, 0.8847, 0.9284, 0.9618, 0.9845, 0.9966. For a surface pressure of 1000 hPa, six levels are below 800 hPa and six levels are above 200 hPa.

Computer/Operating System

The PMIP simulations were run on a Cray C90 computer using a single processor in a UNICOS operating environment.

Computational Performance

For the PMIP simulatons, about 1.4 minutes of Cray C90 computation time per simulated day.


For the PMIP simulations, the model atmosphere is initialized from an idealized, isothermal, dry, motionless state. Soil moisture is initialized at saturation, and snow cover is absent.

Time Integration Scheme(s)

A leapfrog semi-implicit scheme similar to that of Bourke (1974) with an Asselin (1972) frequency filter is used for time integration. The time step is 13 minutes and 20 seconds for dynamics and physics, except for full calculations of all radiative fluxes, which are done once every 24 hours.


Orography is smoothed (see Orography). Negative moisture values (that arise because of spectral truncation) are filled by borrowing from nearest neighbors in the vertical, and horizontally in the east-west direction only.

Sampling Frequency

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

Dynamical/Physical Properties

Atmospheric Dynamics

Primitive-equation dynamics are expressed in terms of vorticity, divergence, surface pressure, specific humidity, and temperature (with a linearized correction for virtual temperature in diagnostic quantities, where applicable).


Linear fourth-order (Ñ4) horizontal diffusion is applied on constant pressure surfaces to vorticity, divergence, temperature, and specific humidity.

Second-order vertical diffusion with coefficients derived from mixing-length considerations is applied to momentum, heat, and moisture is applied at all levels up to a height of about 5 km. The vertical diffusion is not stability-dependent.

Gravity-wave Drag

The parameterization of orographic gravity-wave drag follows linear theory, as formulated by Y. Hayashi (cf. Broccoli and Manabe 1992) . The drag is given by the vertical divergence of the wave stress, which near the surface is equal to the product of the subgrid-scale orographic variance and a representative mountain wavenumber with a mass-weighted average of the atmospheric density, the Brunt-Vaisalla frequency, and the wind in the first three layers above the surface. The wave stress is assumed to be zero above a critical sigma level where the projection of the wind on the surface wind vector vanishes. Below thisi level, the stress increases linearly to its surface value as the sigma-distance from the critical level increases.

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. A seasonal, but not a diurnal cycle in solar forcing, is simulated.


The carbon dioxide concentration used in the PMIP control simulations is 300 ppm. Zonally averaged seasonal mean ozone distributions are specified as a function of height after data of Hering and Borden (1965) . These data are linearly interpolated for intermediate times. Radiative effects of water vapor, but not those of aerosols, also are included (see Radiation).


Shortwave Rayleigh scattering, and absorption in ultraviolet (wavelengths less than 0.35 micron) and visible (wavelengths between 0.5 to 0.7 micron) spectral bands by ozone, and in the near-infrared (wavelengths 0.7 to 4.0 microns) by water vapor are treated by a method similar to that of Lacis and Hansen (1974) . Pressure corrections and multiple reflections between clouds and the surface are treated, but radiative effects of aerosols are not included.

Longwave radiation follows the method of Rodgers and Walshaw (1966) , as modified by Stone and Manabe (1968) . Absorption by water vapor, carbon dioxide, and ozone is included. The 6.3-micron band, the rotation band, and the continuum of water vapor are subdivided into 19 subintervals, two of which contain the 15-micron carbon dioxide and 9.6-micron ozone absorption bands, with transmissions of the latter two gases multiplied together in overlapping bands. A random model is used to represent the absorptivity for each subinterval, and the Curtis-Godson approximation is used to estimate the effective pressure for absorption. The temperature dependence of line intensity is also incorporated. Carbon dioxide absorptivity is obtained from data of Burch et al. (1961) , with its temperature dependence following the model of Sasamori (1959) . The ozone absorptivity is obtained from data of Walshaw (1957) .

Cloud-radiative interactions are treated as described by Wetherald and Manabe (1988) and Wetherald et al. (1991) . Cloud optical properties (absorptivity/ reflectivity in the ultraviolet-visible and near-infrared intervals, and emissivity in the longwave) depend on cloud height (high, middle, and low) and thickness. The values of shortwave properties are assigned following Rodgers (1967a) ; in the longwave, the emissivity of thin high clouds (above 10.5 km) is prescribed as 0.6 after Kondratiev (1972) , while all other clouds are treated as blackbodies (emissivity = 1.0). No partial cloudiness is accounted for in each grid box; clouds, therefore, are treated as fully overlapped in the vertical. See also Cloud Formation.


Convection is simulated via the convective adjustment processes. If the lapse rate exceeds dry adiabatic, a convective adjustment restores the lapse rate to dry adiabatic, with conservation of dry static energy in the vertical. A moist convective adjustment scheme after Manabe et al. (1965) also operates when the lapse rate exceeds moist adiabatic and the air is supersaturated. (For supersaturated stable layers, nonconvective large-scale condensation takes place--see Precipitation). In moist convective layers it is assumed that the intensity of convection is strong enough to eliminate the vertical gradient of potential temperature instantaneously, while conserving total moist static energy. It is further assumed that the relative humidity in the layer is maintained at 100 percent, owing to the vertical mixing of moisture, condensation, and evaporation from water droplets. Shallow convection is not explicitly simulated.

Cloud Formation

Cloud forms when the relative humidity of a vertical layer exceeds a height- dependent threshold, whether this results from large-scale condensation or moist convective adjustment (see Convection). (The threshold relative humidity varies linearly between 100 percent at the bottom of the model atmosphere to 99 percent at its top.) It is assumed that no condensed water is retained as cloud liquid water, but immediately precipitates (see Precipitation). Cloud type is defined according to height: high cloud forms above 10.5 km, middle cloud between 4.0 to 10.5 km, and low cloud between 0.0 to 4.0 km. Cloud may form in a single layer or in multiple contiguous layers; cloud occupying a single layer is treated as radiatively thin, and otherwise as radiatively thick. Clouds are assumed to fill the whole grid box (cloud fraction = 1), and therefore to be fully overlapped in the vertical. Cf. Wetherald and Manabe (1988) and Wetherald et al. (1991) for further details. See also Radiation.


Precipitation from large-scale condensation and from the moist convective adjustment process (see Convection) forms under supersaturated conditions. Subsequent evaporation of falling precipitation is not simulated.

Planetary Boundary Layer

The PBL is simulated only in a rudimentary way (e.g., above the surface layer vertical diffusion of momentum, heat, and moisture are not stability-dependent --see Diffusion). The PBL top is not computed. See also Surface Characteristics and Surface Fluxes.


In the control integrations, raw orography data obtained from the U.S. Navy dataset (cf. Joseph 1980 ) are smoothed by application of the method of Lindberg and Broccoli (1996). The same procedure is applied to the PMIP LGM topography for the 21K integration. Subgrid-scale orographic variances required for the gravity-wave drag parameterization are obtained from the same dataset (see Gravity-wave Drag) and used in all PMIP integrations.


Prescribed SST simulations of the modern climate use the AMIP monthly sea surface temperature fields, with daily values determined by linear interpolation. For the computed SST simulations, a static, isothermal "mixed layer" with a depth of 50 m is employed. An interactive surface heat budget constitues the forcing for this layer, with a heat convergence, varying both spatially and seasonally, added to mimic the process of ocean heat transport.

Sea Ice

Sea ice concentrations and thicknesses used in the prescribed SST simulation are based on analyses from the NOAA/Navy Joint Ice Center and submarine ice thickness data. In the computed SST simulations, the ice model of Flato and Hibler (1992) is used. Snow may accumulate on sea ice, but does not alter its thermodynamic properties. The surface temperature of sea ice is prognostically determined after Holloway and Manabe (1971) from a surface energy balance (see Surface Fluxes) that includes a conduction heat flux from the ocean below. The conduction flux is proportional to the difference between the surface temperature of the ice and the subsurface ocean temperature (assumed to be at the melting temperature of sea ice, or 271.2 K), and the flux is inversely proportional to the constant ice thickness (2 m). The heat conductivity is assumed to be a constant equal to the value for pure ice, and there is no heat storage within the ice.

Snow Cover

Precipitation falls as snow if the air temperature at 350 m above the surface is < 0 degrees C. Snow accumulates on both land and sea ice (see Sea Ice), and snow mass is determined prognostically from a budget equation that accounts for accumulation, melting, and sublimation. Sublimation is calculated as part of the surface evaporative flux. Snowmelt, determined from the excess surface energy available for a surface at temperature 273.2 K, contributes to soil moisture (cf. Holloway and Manabe 1971) . Snow cover affects the albedo of the surface, but not its thermodynamic properties (see also Surface Characteristics, Surface Fluxes, and Land Surface Processes).

Surface Characteristics

Roughness lengths are prescribed constants for ocean (3.5 x 10-4 m) and land (0.045 m).

Over oceans, the surface albedo depends on solar zenith angle (cf. Payne 1972 ). The albedo of sea ice is a function of ice thickness and surface temperature. Albedos for snow-free land are obtained from the modern albedo data base of CLIMAP Project Members (1981) and do not depend on solar zenith angle or spectral interval.

Snow cover modifies the local background albedo of the surface according to its depth, following Holloway and Manabe (1971) . The albedo also depends on surface temperature, with higher albedos at lower temperatures.

Longwave emissivity is prescribed as unity (blackbody emission) for all surfaces.

Surface Fluxes

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

Surface turbulent eddy fluxes follow Manabe (1969). The momentum flux is proportional to the product of a drag coefficient, the wind speed, and the wind velocity vector at the lowest atmospheric level. Surface sensible heat flux is proportional to the product of a transfer coefficient, the wind speed at the lowest atmospheric level, and the vertical difference between the temperature at the surface and that of the lowest level. The drag and transfer coefficients are functions of surface roughness length (see Surface Characteristics) but are not stability-dependent.

The surface moisture flux is the product of potential evaporation and evapotranspiration efficiency beta. Potential evaporation is proportional to the product of the same transfer coefficient as for the sensible heat flux, the wind speed at the lowest atmospheric level, and the difference between the specific humidity at the lowest level and the saturated specific humidity for the local surface temperature and pressure. The evapotranspiration efficiency beta is prescribed to be unity over oceans, snow, and ice surfaces; over land, beta is a function of the ratio of soil moisture to the constant field capacity (see Land Surface Processes).

Land Surface Processes

Ground temperature is determined from a surface energy balance (see Surface Fluxes) without provision for soil heat storage.

Soil moisture is represented by the single-layer "bucket" model of Manabe (1969) , with field capacity everywhere 0.15 m. Soil moisture is increased by precipitation and snowmelt; it is depleted by surface evaporation, which is determined from a product of the evapotranspiration efficiency beta and the potential evaporation from a surface saturated at the local surface temperature and pressure (see Surface Fluxes). Over land, beta is given by the ratio of local soil moisture to a critical value that is 75 percent of field capacity, and is set to unity if soil moisture exceeds this value. Runoff occurs implicitly if soil moisture exceeds the field capacity.

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