the PMIP models (Bonfils et al. 1998)
PMIP Documentation for MRI2
Meteorological Research Institute:
Model MRI GCM-IIb (4x5 L15) 1995
Dr. Akio Kitoh, Climate Research Department, Meteorological Research Institute,
1-1, Nagamine, Tsukuba, Ibaraki, 305-0052 Japan; Phone:+81-298-53-8594;
Fax: +81-298-55-2552; e-mail: email@example.com
Dr Hiroshi Koide, Climate Research Department, Meteorological Research
Institute, 1-1, Nagamine, Tsukuba, Ibaraki, 305-0052 Japan; Phone:+81-298-53-8597
(8681); Fax: +81-298-55-2552; e-mail: firstname.lastname@example.org;
World Wide Web URL: http://www.mri-jma.go.jp.
MRI GCM-IIb (4x5 L15) 1995
Model Identification for PMIP
0fix, 6fix, 21fix, 0cal, 21cal
Number of days in each month: 31 28 31 30 31 30 31 31 30 31 30 31
The MRI model is derived from an earlier version of the University of California
at Los Angeles (UCLA) model (cf. Arakawa and Mintz 1974 and Arakawa and
Lamb 1977 ). Subsequent modifications include increases in vertical resolution,
and changes in parameterizations of gravity-wave drag, atmospheric radiation,
convection, surface characteristics, and land surface processes. The model
is identical to latest AMIP model MRI GCM-IIb (4x5 L15) 1995, except for
different initial conditions and the Earth's orbital parameters.
Documentation for an earlier five-level version of the MRI model is provided
by Tokioka et al. (1984) . Descriptions of subsequent modifications are
given by Yagai and Tokioka (1987) , Yagai and Yamazaki (1988) , Kitoh et
al. (1988) , Tokioka et al. (1988) , Noda and Tokioka (1989), and Shibata
and Aoki (1989) .
Results from the AMIP simulation by the MRI model MRI GCM-IIb (4x5 L15)
1995 are discussed by Kitoh et al. (1995).
Results from the couled atmosphere/mixed-layer ocean model are given
by Kitoh (1997) and Kitoh et al. (1998).
Finite differences on a C-grid with conservation of mass, momentum, energy,
and potential enstrophy (cf. Arakawa and Lamb 1977) .
4 x 5-degree latitude-longitude grid.
Surface to 1 hPa (with the highest prognostic level at 1.39 hPa). For a
surface pressure of 1000 hPa, the lowest atmospheric level is at a pressure
of about 912 hPa. In addition, the depth of the boundary layer is a prognostic
variable--see Planetary Boundary Layer.
Finite differences in hybrid coordinates. Above 100 hPa log-pressure coordinates
are used, and below 100 hPa modified sigma coordinates: sigma = (P - PI)/(PS
- PI), where P and PS are atmospheric and surface pressure, respectively,
and PI is a constant 100 hPa. The vertical differencing scheme is after
Tokioka (1978) .
The 15 hybrid levels in pressure (modified sigma) coordinate are as
the following :
912 (0.902), 741 (0.712), 589 (0.543), 457 (0.397), 335 (0.261), 237
(0.152), 168 (0.076), 119 (0.021),
72 (-0.283), 37.3 (-0.633), 19.3 (-0.815), 10 (-0.909), 5.18 (-0.958),
2.68 (-0.983), 1.39 (-0.996).
15 hybrid levels. For a surface pressure of 1000 hPa, 1 level is below
800 hPa and 9 levels are above 200 hPa. See also Vertical Representation.
The PMIP simulation was run on a HITACHI S3800 using a single processor
in the VOS3 environment.
For the PMIP experiment, about 35 minutes of HITAC S3800 computation time
per simulated month.
Initial conditions of model atmosphere, soil temperature and moisture,
and snow depth/cover for the PMIP 0fix experiment were taken from a time
step of the AMIP simulation, then a time step of the last year of 0fix
run was used for 6fix and 21fix run. The model was span up for several
months for 0fix and 6fix run, and a year for 21fix run.
Time Integration Scheme(s)
The model is integrated by the leapfrog scheme with a time step of 5 minutes,
and with a Matsuno step (an Euler forward integration followed by a backward
integration) inserted hourly. At the forward stage of the Matsuno step,
diabatic and dissipative terms, sources and sinks in atmospheric water
vapor, and the change in the depth of the prognostic planetary boundary
layer (PBL) are calculated. The PBL entrainment rate and the turbulent
fluxes at the PBL top and at the surface are solved iteratively by backward
implicit differencing (see Planetary Boundary Layer and Surface Fluxes).
Shortwave and longwave radiation are computed hourly, but with transmission
functions for longwave fluxes calculated every 3 hours.
Orography is area-averaged (see Orography). A longitudinal smoothing of
the zonal pressure gradient and of the zonal and meridional mass flux is
also performed in high latitudes (cf. Tokioka et al. 1984) . A precise
determination of the horizontal flux of atmosphere moisture and use of
vertical interpolation at half-levels prevents the occurrence of negative
humidity values (cf. Tokioka et al. 1984) , thereby avoiding the need for
For the PMIP simulation, the model history is written every 6 hours.
Primitive-equation dynamics are expressed in flux form in terms of u and
v winds, temperature, specific humidity, and surface pressure.
Nonlinear horizontal diffusion is applied only to momentum (cf. Holloway
and Manabe 1971) .
Vertical diffusion is not applied above the well-mixed PBL (see Planetary
Boundary Layer). However, momentum is redistributed vertically by cumulus
convection (see Convection).
Orographic gravity-wave drag is simulated following Palmer et al. (1986)
with quantitative adjustments described by Yagai and Yamazaki (1988) .
The dependence of the surface momentum flux on the surface wind direction
is considered (see Orography on the computation of the required subgrid-scale
orographic covariances). Surface stress due to gravity waves excited by
stably stratified flow over irregular terrain is calculated from linear
theory and dimensional considerations. The 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. The gravity-wave generation factor kappa
is set to 1.2 x 10-5 m-1, which is same as the latest
AMIP simulation with MRI GCM-IIb (4x5 L15) 1995.
The solar constant is the PMIP-prescribed value of 1365 W/(m2).
The orbital parameters and seasonal insolation distribution are calculated
after PMIP recommendation. Both seasonal and diurnal cycles in solar forcing
The carbon dioxide concentration is the PMIP-prescribed value of 345, 280,
200 ppm for control, 6k and 21k run, respectively. Zonal profiles of ozone
concentration are prescribed monthly from data of McPeters et al. (1984)
. The radiative effects of water vapor are also treated, but not those
of aerosols (see Radiation).
Shortwave Rayleigh scattering and absorption in ultraviolet (wavelengths
less than 0.35 micron) and visible (wavelengths 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 and carbon dioxide, follows the method of Lacis and Hansen
(1974) . Pressure corrections and multiple reflections between randomly
overlapped partial clouds and the surface are treated.
Longwave calculations are based on the multiparameter random model of
Shibata and Aoki (1989) applied in four spectral regions (with boundaries
at 2.0 x 103, 5.5 x 104, 8.0 x 104, 1.2
x 105, and 2.2 x 105 m-1). Absorption
bands of carbon dioxide, ozone, and water vapor are included. Continuum
absorption of water vapor is treated after the method of Clough et al.
(1980) , but with the temperature dependency used by Roberts et al. (1976)
. Mean diffuse transmittances in each region are determined from line-by-line
calculations approximated by an exponential function including 10 or 12
parameters. Transmittances through inhomogeneous atmospheres are computed
by a modified Godson (1953) method.
In the shortwave, clouds are treated by a delta Eddington approximation
with prescribed single-scattering albedo and asymmetry factor, and with
cloud optical depth a function of height. In the longwave, clouds behave
as blackbodies, except for high (above 400 hPa) clouds, whose emissivity
is related to shortwave optical depth. For purposes of the radiation calculations,
all clouds are assumed to be randomly overlapped in the vertical. See also
Simulation of cumulus convection is by a modified Arakawa and Schubert
(1974) scheme, as implemented by Lord (1978) , Lord and Arakawa (1980)
, and Lord et al. (1982) . Convective mass fluxes are assumed to originate
in the PBL (see Planetary Boundary Layer); these are predicted from mutually
interacting cumulus subensembles, which have different entrainment rates
and levels of neutral buoyancy that define the tops of the clouds and their
associated convective updrafts. In turn, the convective mass fluxes feed
back on the large-scale fields of temperature (through latent heating and
compensating subsidence), moisture (through precipitation and detrainment),
and momentum (through cumulus friction). The effects on cloud buoyancy
of phase changes from water to ice and the drying and cooling effects of
convective-scale downdrafts on the environment are not included explicitly.
The Arakawa-Schubert scheme is modified to impose an additional constraint
between the minimum entrainment rate and the prognostic depth of the PBL
(cf. Tokioka et al. 1988) . The mass flux for each cumulus subensemble
is predicted from an integral equation that includes a positive-definite
work function (defined by the tendency of cumulus kinetic energy for the
subensemble) and a negative-definite kernel which expresses the effects
of other subensembles on this work function. The predicted mass fluxes
are optimal solutions of this integral equation under the constraint that
the rate of generation of conditional convective instability by the large-scale
environment is balanced by the rate at which the cumulus subensembles suppress
this instability via large-scale feedbacks (cf. Lord et al. 1982) . The
mass fluxes are computed by the "exact direct method," which guarantees
an exact solution within roundoff errors.
If the lapse rate becomes dry convectively unstable at any level, moisture
and enthalpy are redistributed vertically. In addition, a moist convective
adjustment simulates midlevel convection originating above the PBL. When
the lapse rate exceeds moist adiabatic under supersaturated conditions,
mass is mixed such that either the lapse rate is restored to moist adiabatic,
or the supersaturation is eliminated by formation of convective precipitation.
Cf. Tokioka et al. (1984) for further details.
Five types of cloud are simulated: penetrative and midlevel convective
cloud, cirrus anvil cloud, large-scale condensation cloud, and PBL stratus
cloud. The PBL stratus is the only cloud type that does not interact with
radiation. The penetrative and midlevel convective cloud are associated,
respectively, with the cumulus convection and moist convective adjustment
schemes (see Convection). Convective clouds in a vertical column are assumed
to overlap randomly, with a total cloud fraction of 0.3. Cirrus anvil cloud
forms if cumulus convection penetrates to levels above 400 hPa, and large-scale
condensation cloud is present if the local relative humidity of a layer
is at least 100 percent. The latter two cloud types cover the grid box
(cloud fraction = 1). Stratus cloud capping the PBL forms if the specific
humidity at the PBL top is greater than saturation (see Planetary Boundary
Layer). The base of the cloud is determined from the vertical distribution
of temperature and specific humidity in the PBL. See also Radiation for
treatment of cloud-radiative interactions.
Large-scale precipitation results from condensation under supersaturated
conditions. Precipitation from cumulus convection originating in the PBL
is simulated by the Arakawa-Schubert (1974) scheme, and precipitation from
midlevel convection by a moist adjustment process (see Convection). Both
large-scale and convective precipitation may evaporate in falling through
lower layers; the amount of evaporation is equal to that required to saturate
each of these layers in turn.
Planetary Boundary Layer
The PBL is parameterized as a well-mixed layer after Randall (1976) . The
depth and mean structure of the PBL are prognostic functions of the surface
momentum, heat, and moisture fluxes (see Surface Fluxes), as well as horizontal
mass convergence, entrainment, and cumulus mass fluxes determined by the
cumulus convection scheme (see Convection). PBL u-v winds, temperature,
and specific humidity are predicted. Discontinuities in atmospheric variables
that may exist at the PBL top because of the absence of vertical diffusion
in the free atmosphere (see Diffusion) are also predicted. In addition,
PBL stratus cloud is produced under saturated conditions (see Cloud Formation),
and its evaporation by drier entrained air is treated.
The turbulence kinetic energy (TKE) also is determined. (TKE is generated
by wind shear and convective buoyancy fluxes, and is depleted by surface
dissipation, by work done against the free atmosphere, and by newly turbulent
air that is entrained into the PBL.) Because of the mutual dependence of
the entrainment rate and the turbulent fluxes at the PBL top and at the
surface, these quantities are solved iteratively under the assumption that
the generation of TKE balances its dissipation. See also Surface Characteristics
and Time Integration Scheme(s).
Orography, originally from US Navy 10-minute resolution dataset (Joseph
1980), is area averaged over each 4 x 5-degree model grid box is used 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 covariances required for the parameterization of gravity-wave
drag (see Gravity-wave Drag) are computed from the U.S. Navy 10-minute
resolution topography dataset (cf. Joseph 1980) . The same method was used
for 21ka after addition of the orographic difference.
The prescribed monthly climatorogical SSTs for 0fix and 6fix integrations
were made by averaging AMIP monthly sea surface temperature fields for
1979 to 1988, with daily values determined by linear interpolation.
For 21fix run, the CLIMAP SSTs ( for February and August ) are used
as a mid-month values, and interpolated with sine function. The SST of
the region where temporally covered with sea ice is appropriately interpolated
with assuming SST of the ice edge is -1.8 degrees C (see Sea Ice).
For 0cal and 21 cal, SST is calculated in every 5 day from the heat
balance of the 50-m depth mixed-layer (slab) ocean using the heat fluxes
at the sea surface and prescribed heat flux value (Q-flux). Monthly mean
Q-flux value is obtainedby runing the model for 3 years restoring the model
SST and sea-ice amount to observations. The same Q-flux value is used for
0cal and 21cal with some adjustment over the sea-ice area in 21k.
For 0fix and 6fix run, sea ice concentration fraction (ranging between
0 and 1) for each grid box is prescribed monthly by averaging over corresponding
weekly U.S. Navy and National Oceanic and Atmospheric Administration (NOAA)
Joint Ice Center data for 1979 to 1988.
For 21fix run, see Ocean. The spatially variable thickness of sea ice
is given by the local concentration fraction multiplied by 3.0 meters in
the Northern Hemisphere and by 0.5 meter in the Southern Hemisphere. Daily
values of the above quantities are determined by linear interpolation.
The ice surface temperature is predicted from the net flux of energy
(see Surface Fluxes), including a subsurface conduction heat flux that
is proportional to the difference between the ice surface temperature and
that prescribed (271.3 K) for the ocean below. Snow may accumulate on sea
ice or melt if the snow surface temperature exceeds 0 degrees C (see Snow
Cover). Snow also alters the heat capacity /conductivity of the ice, but
the heat capacity of snow is independent of depth.
For 0cal and 21cal, sea-ice concentration (compactness) and thinckness
are also predicted in every 5 day following Mellor and Kantha (1989). The
sea-ice model is basically the same one used in MRI coupled GCM (Tokioka
et al., 1996), but sea surface salinity is fixed at the present climatology
and advection of sea-ice is not included.
If the surface air temperature is <0 degrees C, precipitation falls
as snow. Snow may accumulate on both land and sea ice, with complete coverage
of a grid box assumed (i.e., there is no fractional snow coverage). Prognostic
snow mass is depleted both by snowmelt (which contributes to soil moisture)
and by sublimation (which contributes to surface evaporation). Snow melts
if its surface temperature exceeds 0 degrees C. Snow cover affects both
surface albedo and the thermal properties of the surface. See also Sea
Ice, Surface Characteristics, Surface Fluxes, and Land Surface Processes.
The surface type is defined as open ocean, sea ice, glacial ice, lake,
or land. (Lakes are treated as mixed layers of depth 10 m, with prognostic
The surface roughness length is a fixed value over oceans (2 x 10-4
m), sea ice (1 x 10-4 m), glacial ice (5 x 10-3 m),
and land (0.45 m); the surface drag coefficient over land is a function
of orographic variances, however (see Orography and Surface Fluxes).
Surface albedos depend on solar zenith angle (cf. Paltridge and Platt
1976) , but not on spectral interval. The ocean albedo is a maximum of
0.07. The albedo of bare sea ice ranges between a maximum of 0.50 to 0.64
as a function of surface temperature; the albedo of snow-covered sea ice
is a maximum of 0.70, depending on snow depth. Snow-covered glacial ice
albedos range between 0.70 and 0.85. Snow-free land albedos are obtained
from the data of Matthews (1983) ; with snow cover, the land albedo ranges
from its snow-free value to a maximum between 0.60 to 0.70, depending on
topographic height (see Orography). On all surfaces, the albedo of melting
snow is 0.60, and that of frost is 0.30. Longwave emissivity is prescribed
as unity (blackbody emission) for all surfaces.
Surface solar absorption is determined from albedos, and longwave emission
from the Planck equation with emissivity of 1.0 (see Surface Characteristics).
Turbulent eddy fluxes of momentum, heat, and moisture are parameterized
as bulk formulae with drag and transfer coefficients that depend on vertical
stability (bulk Richardson number) and the (locally variable) depth of
the PBL normalized by the surface roughness length, following Deardorff
(1972) . The drag coefficient over land is also increased as a function
of orographic variances (see Orography and cf. Yagai and Tokioka 1987 ).
The surface atmospheric values of wind, dry static energy, and humidity
required for the bulk formulae are taken to be those predicted in the PBL
(see Planetary Boundary Layer). Because of the mutual dependence of the
turbulent fluxes at the surface and at the PBL top, these (as well as the
PBL entrainment rate) are solved by mutual iteration. See also Surface
Characteristics and Time Integration Scheme(s). The surface moisture flux
depends also on an evapotranspiration efficiency factor beta, which is
set to unity over ocean and ice surfaces and in areas of continental dew
formation; otherwise, over land, beta is a function of soil moisture (see
Land Surface Processes).
Land Surface Processes
Land-surface parameterizations follow Katayama (1978) , with subsequent
modifications in the treatment of soil moisture described by Kitoh et al.
(1988) and inclusion of multi-layer model (Noda et al., 1995). Soil heat
and moisture move along the gradients of temperature and wetness, respectively.
The thermodynamic and hydrologic effects of a vegetation canopy are not
explicitly modeled, however. The temperature of bare and snow-covered land
is predicted in four layers, and that of glacial ice in a single-layer.
The upper boundary condition is the net balance of surface energy fluxes
(see Surface Fluxes), while there is zero heat transfer at the lower boundary.
Heat exchange associated with the freezing or melting of soil moisture
and interstitial ice is taken into account. Snow cover and soil moisture/ice
also affect the heat capacity /conductivity of the surface, but the heat
capacity of snow is assumed to be independent of its depth. Soil moisture
(in both liquid and frozen form) is predicted in four layers (with bottom
boundaries at 0.10, 0.50, 1.50, and 10.0 m below the surface); it is augmented
by precipitation and snow/ice melt, and is depleted by surface evaporation
and runoff. The evapotranspiration efficiency beta (see Surface Fluxes)
is set to unity if the fractional soil moisture (the ratio of soil moisture
to the field capacity, assumed to be 20 percent of the volume of a soil
column) is at least 0.5; beta is set to twice the fractional soil moisture
otherwise. When all of the moisture in a soil layer is completely frozen
(freezing begins when the layer temperature falls below 0 degrees C), no
deeper penetration of moisture is allowed. Runoff occurs in any layer that
is either saturated or completely frozen.
Last update November 9, 1998. For further information, contact: Céline