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

PMIP Documentation for CCSR1

Center for Climate System Research: Model CCSR/NIES 5.4 02 T21/L20 1995


PMIP Representative(s)

Dr. Ayako Abe-Ouchi, Center for Climate System Research (CCSR), University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo, 153 Japan, Phone: +81-3-5453-3955; Fax: +81-3-5453-3964; e-mail:

WWW URL: (in Japanese);

Model Designation

CCSR/NIES AGCM 5.4 02 (T21 L20) 1995

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

Model CCSR/NIES AGCM (T21 L20) 1995 is based on a simple global atmospheric model first developed at the University of Tokyo (cf. Numaguti 1993), and further refined as a collaboration between CCSR and the National Institute of Environmental Studies (NIES). It is intended for use as a community climate model.

The model is identical to the latest AMIP model exept for different initial conditions and the Earth's orbital parameters.

Model Documentation

Numaguti and others (1996) Development of an Atmospheric General Circulation Model.(prepared for submittion. available upon request)

A summary of model features including fundamental equations is provided by Numaguti et al. (1995). The spectral formulation of atmospheric dynamics follows closely Bourke (1988). The radiation scheme is described by Nakajima and Tanaka (1986) and Nakajima et al. (1996). The convective parameterization is based on the work of Arakawa and Schubert (1974) and Moorthi and Suarez (1992). Cloud formation is treated prognostically after the method of Le Treut and Li (1991). Gravity-wave drag is parameterized as in McFarlane (1987). The planetary boundary layer (PBL) is simulated by the turbulence closure scheme of Mellor and Yamada (1974, 1982)[10,11]. The representation of surface fluxes follows the approach of Louis (1979), with inclusio of adjustments recommended by Miller et al. (1992) for low winds over the oceans.

Numerical/Computational Properties

Horizontal Representation

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

Horizontal Resolution

Spectral triangular 21 (T21), roughly equivalent to a 5.6 x 5.6 degree latitude/longitude grid. dim_longitude*dim_latitude: 64*32

Vertical Domain

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

Vertical Representation

Sigma coordinates with discretization following the vertical differencing scheme of Arakawa and Suarez (1983) that conserves global mass integrals of potential temperature and total energy for frictionless adiabatic flow .

Vertical Resolution

There are 20 unevenly spaced sigma levels. For a surface pressure of 1000 hPa, 5 levels are below 800 hPa and 8 levels are above 200 hPa. The model uses a sigma coordinate in the vertical as: 1: 0.99500 2: 0.97999 3: 0.94995 4: 0.89988 5: 0.82977 6: 0.74468 7: 0.64954 8: 0.54946 9: 0.45447 10: 0.36948 11: 0.29450 12: 0.22953 13: 0.17457 14: 0.12440 15: 0.08468 16: 0.05980 17: 0.04493 18: 0.03491 19: 0.02488 20: 0.00830

Computer/Operating System

The PMIP simulation was run on a HITAC S-3800 computer using a single processor in the VOS3 operational environment.

Computational Performance

For the PMIP experiment, about 0.3 minutes of HITAC S-3800 computation time per simulated day.


Initial conditions of model atmosphere, atmospheric state, soil moisture, and snow cover/depth were obtained from a AMIP run.

Time Integration Scheme(s)

Semi-implicit leapfrog time integration with an Asselin (1972) time filter. The time step length is 40 minutes. Shortwave and longwave radiative fluxes are recalculated every 3 hours, but with the longwave fluxes assumed constant over the 3-hour interval, while the shortwave fluxes are assumed to vary as the cosine of the solar zenith angle.


Orography is smoothed (see Orography). Spurious negative atmospheric moisture values are filled by borrowing from the vertical level immediately below, subject to the constraint of conservation of global moisture.

Sampling Frequency

For the PMIP simulation, the model history is written once per 24-hour period.

Dynamical/Physical Properties

Atmospheric Dynamics

Primitive equation dynamics are expressed in terms of vorticity and divergence,temperature, specific humidity, cloud liquid water, and surface pressure.


Eighth-order linear (Ñ8) horizontal diffusion is applied to vorticity, divergence, temperature, specific humidity, and cloud liquid water on constant sigma surfaces. Stability-dependent vertical diffusion of momentum, heat, and moisture in the planetary boundary layer (PBL) as well as in the free atmosphere follows the Mellor and Yamada (1974, 1982)[10,11 ] level-2 turbulence closure scheme. The eddy diffusion coefficient is diagnostically determined as a function of a Richardson number modified to include the effects of condensation. The diffusion coefficient also depends on the vertical wind shear and on the square of an eddy mixing length with an asymptotic value of 300 m. Cf. Numagati et al. (1995) for further details. See also Planetary Boundary Layer and Surface Fluxes.

Gravity-wave Drag

Orographic gravity-wave drag is parameterized after McFarlane (1987). Deceleration of the resolved flow by dissipation of orographically excited gravity waves is a function of the rate at which the parameterized vertical component of the gravity-wave momentum flux decreases in magnitude with height. This momentum-flux term is the product of local air density, the component of the local wind in the direction of that at the near-surface reference level, and a displacement amplitude. At the surface, this amplitude is specified in terms of the mesoscale orographic variance, and in the free atmosphere by linear theory, but it is bounded everywhere by wave saturation values. 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 carbon dioxide concentration is the AMIP-prescribed value of 345, 280 and 200 ppm for 0fix, 6 fix and 21fix run, respectively. Radiative effects of water vapor, oxygen, ozone, nitrous oxide (0.3 ppm, globally uniform), and methane (1.7 ppm, globally uniform) are included. Monthly zonal ozone profiles are specified from data of Keating and Young (1985) and Dütsch (1978), and they are linearly interpolated for intermediate time points. Although the model is able to treat radiative effects of aerosols, they are not included for the PMIP simulation. See also Radiation.


The radiative transfer scheme is based on the two-stream discrete ordinate method (DOM) and the k-distribution method described in detail by Nakajima et al. (1996)).The radiative fluxes at the interfaces of each vertical layer are calculated considering solar incidence, absorption, emission, and scattering by gases and clouds, with the flux calculations being done in 18 wavelength regions. Band absorption by water vapor, carbon dioxide, ozone, nitrous oxide, and methane is considered in from 1 to 6 subchannels for each wavelength region. Continuum absorption by water vapor, oxygen, and ozone also is included. Rayleigh scattering by gases and absorption by clouds are considered as well. See also Chemistry.

The radiative flux in each wavelength region is calculated as a sum of the products of the fluxes over the subchannels and their respective k-distribution weights, where each subchannel's fluxes are calculated by the two-stream DOM. The optical depth of each subchannel is estimated as the sum of the optical thicknesses of band absorption and of continuum absorption by gases. The transmissivity, reflectivity, and source function in each layer then are calculated as functions of optical depth, single-scattering albedo, asymmetry factor, cutoff factor, Planck function, solar incidence, and solar zenith angle. At each layer interface, fluxes are computed by the adding technique.

In the presence of clouds, radiative fluxes are weighted according to the convective and large-scale cloud fractions of each grid box. The fluxes are computed by treating clouds as a mixture of scattering and absorbing water and ice particles in the shortwave; shortwave optical properties and longwave in the emissivity are functions of optical depth. Fluxes therefore depend on the prognostic liquid water content (LWC) as well as on the fraction of ice cloud (see Cloud Formation). Radiative transfer in large-scale and convective clouds is treated separately, assuming random and full overlap, respectively, in the vertical. Cf. Numaguti et al. (1995) for further details.


Penetrative and shallow cumulus convection are simulated by the Relaxed Arakawa-Schubert (RAS) scheme of Moorthi and Suarez (1992) , a modification of the Arakawa and Schubert (1974) parameterization. The RAS scheme predicts mass fluxes from a spectrum of clouds that have different entrainment/detrainment rates and levels of neutral buoyancy (i.e., different cloud-top heights). The thermodynamic properties of the convective clouds are determined from an entraining plume model and the vertical profile of cloud liquid water (see Cloud Formation) is calculated from the difference between the adiabatic total water mixing ratio (a function of the grid-scale specific humidity) and the saturated specific humidity at the same level, given a prescribed vertical profile of precipitation.

The predicted convective mass fluxes are used to solve budget equations that determine the impact of convection on the grid-scale fields of temperature (through latent heating and compensating subsidence) and moisture (through precipitation and detrainment). The vertical mass flux at the base of each cloud type is predicted from the cloud work function A, defined as the integral over the cloud depth of the product of the mass flux (with a linear vertical profile assumed) and the buoyancy (proportional to the difference between the cloud virtual temperature and that of the grid-scale environment at the same height). Because a nonzero cloud-base mass flux implies a positive-definite work function, the former is determined assuming that the work function vanishes in a specified time scale T > that is longer than the convective time step t of 80 minutes.

In the RAS scheme, the new cloud-base mass flux at each time step is estimated by the method of virtual displacement: the amount of grid-scale warming and drying expected from a unit mass flux is calculated, and a new cloud work function A' is determined; the new cloud-base mass flux M' then is derived from a simple proportionality relation. The grid-scale mass flux is obtained by summing over the contributions from the spectrum of cloud types. The profile of the mass flux associated with convective downdrafts is also simulated from a fixed fraction of the evaporation of convective precipitation (see Precipitation). Updated values of convective cloud fraction and convective liquid water content (LWC) for the grid box also are determined from the grid-scale mass flux (see Cloud Formation and Radiation). Cf. Numaguti et al. (1995) for further details.

Cloud Formation

The convective cloud fraction in a grid box is estimated as proportional to the grid-scale convective mass flux. The grid-scale liquid water content (LWC) at a given height due to convective cloud is determined by a sum over the cloud-type spectrum of the products of LWC and mass flux for each cloud type (see Convection) Large-scale (stratiform) cloud formation is determined from prognostic cloud liquid water content (LWC) following Le Treut and Li (1991) . The stratiform LWC follows a conservation equation involving rates of large-scale water vapor condensation, evaporation of cloud droplets, and the transformation of small droplets to large precipitating drops (see Precipitation). The stratiform LWC (including ice content) also determines the large-scale cloud fraction (see below) and cloud optical properties (see Radiation).

The fraction of stratiform cloud C in any layer is determined from the probability that the total cloud water (liquid plus vapor) is above the saturated value, where a uniform probability distribution with prescribed standard deviation is assumed. (For purposes of the radiation calculations, the square root of C is taken as the cloud fraction). At each time step, new values of LWC and vapor are determined by iteration, subject to conservation of moist internal energy. The portion of C that is ice cloud is assumed to vary as a linear function of the temperature departure below the freezing point 273.15 K, with all of C being ice cloud if the temperature is < 258.18 K. Cf. Numaguti et al. (1995) and Le Treut and Li (1991) for further details.


The autoconversion of cloud liquid water into precipitation is estimated from the prognostic liquid water content (LWC) divided by a characteristic precipitation time scale which is an exponential function of temperature (see Cloud Formation). Precipitation conversion is distinguished for liquid vs ice particles. Snow is assumed to fall when the local wet-bulb temperature is less than the freezing point of 273.15 K, with melting of falling snow occurring if the wet-bulb temperature exceeds this value. See also Snow Cover. Falling liquid precipitation evaporates proportional to the difference between the saturated and ambient specific humidities and inversely proportional to the terminal fall velocity (cf. Kessler (1969)). Falling ice and snow melts if the ambient wet-bulb temperature exceeds the freezing point (273.15 K); evaporation may follow, as for falling liquid precipitation.

Planetary Boundary Layer

The Mellor and Yamada (1974, 1982) [10,11] level-2 turbulence closure scheme represents the effects of the PBL. The scheme is used to determine vertical diffusion coefficients for momentum, heat, and moisture from the product of the squared mixing length (whose asymptotic value is 300 m), the vertical wind shear, and a Richardson number that is modifed to include the effect of condensation on turbulent fluxes. (A diffusion coefficient is never allowed to fall below 0.15 m2/s.) Cf. Numaguti et al. (1995) for further details. See also Surface Fluxes.


Raw orography is obtained from the ETOPO5 dataset (cf. NOAA/NGDC, 1989) at a resolution of 5 x 5 minutes for 0fix and 6fix run. It is modified after PMIP recommendation for 21fix run using the difference of the orographic difference between present and 21ka. Orographic variances required for the gravity-wave drag scheme are obtained from the same dataset. Orography is smoothed by first expanding the grid point data in spectral space, then filtering according to the formula [1-(n/N)4], where n is the spectral wavenumber and N = 21 corresponds to the horizontal resolution of the model. Finally, the smoothed spectral data is returned to the T21 Gaussian grid. See also Gravity-wave Drag.


The prescribed monthly climatological SSTs for 0fix and 6fix integrations were made by averaging AMIP monthly sea surface temperature fields, with daily values determined by linear interpolation.

For 21fix ? ? ? ?

Sea Ice

For 0fix and 6fix run, monthly AMIP sea ice extents are prescribed. The thickness of the ice can vary: the local thickness is determined from the observed fractional coverage multiplied by a constant 1 m. The surface temperature of the ice is predicted from a surface energy balance that takes account of conduction heating from the ocean below. The temperature of the underlying ocean is assumed to be 273.15 K, the freezing point of the sea ice. Snow may accumulate on sea ice, but modifies only the thermal conductivity of the ice. See also Surface Fluxes and Snow Cover.

For 21fix run, sea ice edge of CLIMAP data is used for Feb. and Aug., the shape of sea ice edge for other months is determined by considering summer SST over both hemisphere. Grid with lower summer SST is freeze earlier and melt later.

Snow Cover

Precipitation falling on a surface with skin temperature < 273.15 K accumulates as snow, and a snowpack melts if the skin temperature exceeds this value. Fractional coverage of a grid box is determined by the ratio of the local snow mass to a critical threshold of 200 kg/(m2). Sublimation of snow contributes to the surface evaporative flux (see Surface Fluxes), and snowmelt augments soil moisture and runoff (see Land Surface Processes). Snow cover alters the evaporation efficiency and permeability of moisture, as well as the albedo, roughness, and thermal properties of the surface (see Surface Characteristics).

Surface Characteristics

The surface is classified according to the 32 vegetation types of Matthews (1983), but with only the locally dominant type specified for each grid box. The stomatal resistance of the vegetation is a prescribed spatially uniform value, but is set to zero in desert areas. Over ice surfaces, the roughness length is a constant 1 x 10-3 m. Over ocean, the roughness length is a function of the surface momentum flux, following the formulation of Miller et al. (1992). Over land, roughness lengths are assigned according to vegetation type following Takeuchi and Kondo (1981). IN areas with snow cover, the roughness length is decreased proportional to the square root of the fractional snow cover. The roughness length for calculation of surface momentum fluxes is 10 times the corresponding value for heat and moisture fluxes.

Over ice surfaces, the albedo is a constant 0.7 (unaffected by snow accumulation). Over ocean, the albedo depends on sun angle and the optical thickness of the atmosphere. The albedos of the land surface are specified according to vegetation type from the data of Matthews (1983). For snow-covered land, the albedo increases over that of the background proportional to the square root of the fractional snow cover. Longwave emissivity is everywhere specified to be 1.0 (i.e., blackbody emission). See also Surface Fluxes and Land Surface Processes.

Surface Fluxes

Solar absorption at the surface is determined from the albedo, and longwave emission from the Planck equation with prescribed emissivities (see Surface Characteristics). The representation of turbulent surface fluxes of momentum, heat, and moisture follows Monin-Obukhov similarity theory as expressed by the bulk formulae of Louis (1979). The requisite wind, temperature, and humidity values are taken to be those at the lowest atmospheric level (see Vertical Domain).

The associated drag/transfer coefficients are functions of the surface roughness (see Surface Characteristics) and vertical stability expressed as a function of a modified Richardson number (see Planetary Boundary Layer). The effect of free convective motion is incorporated into the surface wind speed following Miller et al. (1992), and the surface wind speed also is not allowed to fall below 4 m/s.

For calculation of the moisture flux over ocean, ice, and snow-covered surfaces, the evaporation efficiency beta is unity. Over partially snow-covered grid boxes, beta increases as the square root of the snow fraction (see Snow Cover). The evaporation efficiency over vegetation is limited by the specified stomatal resistance. Cf. Numaguti et al. (1995) for further details. See also Land Surface Processes).

Land Surface Processes

The skin temperature of soil and land ice is predicted by a heat diffusion equation that is discretized in 3 layers with a zero-flux lower boundary condition; heat capacity and conductivity are spatially uniform values. Surface snow is treated as part of the uppermost soil layer, and thus modifies its heat content, as well as the heat conduction to lower layers.

Soil liquid moisture is predicted in a single layer according to the "bucket" formulation of Manabe et al. (1965)) The moisture field capacity is a spatially uniform 0.15 m, with surface runoff occurring if the predicted soil moisture exceeds this value. Snowmelt contributes to soil moisture, but if snow covers a grid box completely, the permeability of the soil to falling liquid precipitation becomes zero. For partial snow cover, the permeability decreases proportional to increasing snow fraction (see Snow Cover).

Soil moisture is depleted by surface evaporation; the evaporation efficiency beta (see Surface Fluxes) is not determined solely by the ratio of soil moisture to its saturation value, but is limited by the specified stomatal resistance of the vegetation. Other effects of vegetation, such as the interception of precipitation by the canopy and its subsequent reevaporation, are not included. See also Surface Characteristics and Surface Fluxes.

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