the PMIP models (Bonfils et al. 1998)
PMIP Documentation for LMD4
Laboratoire de Météorologie
Dynamique: Model LMCE LMD4.3 (sin(lat)x7.5 L11) 1991
for 0fix and 6fix:
Dr. Pascale Braconnot, Laboratoire des Sciences du Climat et de l´Environnement,
Bat 709 CEA-DSM, Orme des Merisiers, F 91191 Gif-sur-Yvette cedex, France,
Phone: 33 (1) 69 08 77 11; Fax 33 (1) 69 08 77 16;
Dr. Gilles Ramstein Laboratoire des Sciences du Climat et de l´Environnement,
Bat 709 CEA-DSM, Orme des Merisiers, F 91191 Gif-sur-Yvette cedex, France,
Phone: 33 (1) 69 08 77 11; Fax 33 (1) 69 08 77 16;
LMD4.3 (sin(lat)x7.5 L11) 1991
Model Identification for PMIP
0fix, 6fix, 21fix, 0cal, 21cal
Number of days in each month: 30 30 30 30 30 30 30 30 30 30 30 30
The LMD model derives from an earlier version developed for climate studies
(cf. Sadourny and Laval 1984) . Subsequent modifications principally include
changes in the representation of radiation and horizontal diffusion, Inclusion
of parameterizations of gravity-wave drag and prognostic cloud formation
and cloud water pronostic equation (Le Treut 1994 and Ramstein 1998).
Sadourny R., Laval K., 1984, January and July performance of the LMD general
circulation model. In New Perspectives in Climate Modelling, A. Berger
and C. nicolis (eds), Developments in Atmospheric Science, 16, Elsevier,
To see the improvements between LMD4ter and LMD4 (previous version)
Le Treut H., Z.X. Li and M. Forichon, 1994 : Sensitivity of the LMD
general circulation model to greenhouse forcing associated with two different
cloud water parametrizations, J. Climate, Vol. 7, 1827-1841. .
Z.X. Li and Le Treut, H, 1992, Cloud radiation feedbacks in a general
circulation model and their dependence on cloud modelling assumptions,
Climate Dynamics, Vol. 7, 133-139.
Ramstein, G., Serafini, Y. V., Le Treut, H., 1998, Cloud processes associated
with past and future climate changes, Climate Dynamics, in press.
Finite differences on a uniform-area (constant area), staggered C-grid
(cf. Arakawa and Lamb 1977) .
The model has a regular grid in longitude and in sine of the latitude
(ie. 7.5° for longitude and latitude between 3° (equator) and 13°
There are 36 grid points equally spaced in the sine of latitude and 48
points equally spaced in longitude. (The mesh size is 335 km north-south
and 837 km east-west at the equator, and is about 536 x 592 km at 45 degrees
Surface to about 4 hPa. For a surface pressure of 1000 hPa, the lowest
atmospheric level is at 979 hPa.
Finite-difference sigma coordinates.
There are 11 unevenly spaced sigma levels: the model uses a sigma-coordinate
in the vertical with the following levels: 0.991, 0.967, 0.914, 0.830,
0.708, 0.566, 0.411, 0.271, 0.150 0.07, 0.01.
For a surface pressure of 1000 hPa, 3 levels are below 800 hPa and 2
levels are above 200 hPa.
The PMIP simulations were performed using a CRAY C94 under the UNICOS operating
For the PMIP experiment, about 1.5 hours Cray C94 computation time per
For the PMIP experiment, atmosphere, soil moisture, and snow cover/depth
are initialized for 1 January 1979 from a previous model simulation.
For 21k, the initialisation is done with a second year first day run
of the control run. Then the physics part is inhibit, only the dynamics
part is run for 10 days which is the time needed by the model to dump the
waves generated by the ice-sheet LGM. Then the physics part is switch on
at the eleventh day. The run lasts 16 years, we keep only the 15 last years.
For 0k computed SSTs, the flux are computed from the fixed SSts control
run and prescribed to the computed SSts control run for 30 years, when
the equilibrium is reached. Then 15 years are performed to be analysed
as control SST computed run.
For 21k comuted SSTs, the flux from the fixed SSTs control run are corrected
for the sea level drop. For a first step, 30 years are performed to reach
the equilibrium and 15 new years are performed to be analysed as LGM SST
Time Integration Scheme(s)
The time integration scheme for dynamics combines 4 leapfrog steps with
a Matsuno step, each of length 6 minutes. In the model, physics is updated
every 30 minutes, except for shortwave/longwave radiative fluxes, which
are calculated every 6 hours. For computation of vertical turbulent surface
fluxes and diffusion, an implicit backward integration scheme with 30-minute
time step is used, but with all coefficients calculated explicitly. See
also Surface Fluxes and Diffusion.
Orography is averaged on the model grid (see Orography). At the four latitude
points closest to the poles, a Fourier filtering operator from Arakawa
and Mintz (1974) is applied to the momentum, thermodynamics, continuity,
and water vapor tendency equations to slow the longitudinally propagating
gravity waves for numerical stability. Negative moisture values (arising
from vertical advection by the centered non diffusive scheme) are filled
by borrowing moisture from the level below.
For the PMIP simulation, the model history is written once every 24 hours.
Primitive-equation dynamics are expressed in terms of u and v winds, potential
enthalpy, specific humidity, and surface pressure. The advection scheme
is designed to conserve potential enstrophy for divergent barotropic flow
(cf. Sadourny 1975a , b ). Total energy is also conserved for irrotational
flow (cf. Sadourny 1980) . The continuity and thermodynamics equations
are expressed in flux form, conserving mass and the space integrals of
potential temperature and its square. The water vapor tendency is also
expressed in flux form, thereby reducing the probability of spurious negative
moisture values (see Smoothing/Filling).
Linear horizontal diffusion is applied on constant-pressure surfaces to
potential enthalpy, divergence, and rotational wind via a biharmonic operator
, where Ñ denotes a first-order difference
on the model grid, while Ñ * is a formal
differential operator on a regular grid without geometrical corrections.
Because of the highly diffusive character of the flux-form water vapor
tendency equation (see Atmospheric Dynamics), no further horizontal diffusion
of specific humidity is included. Cf. Michaud (1987) for further details.
Second-order vertical diffusion of momentum, heat, and moisture is applied
only within the planetary boundary layer (PBL). The diffusion coefficient
depends on a diagnostic estimate of the turbulence kinetic energy (TKE)
and on the mixing length (which decreases up to the prescribed PBL top)
that is estimated from Smagorinsky et al. (1965) . Estimation of TKE involves
calculation of a countergradient term from Deardorff (1966) and comparison
of the bulk Richardson number with a critical value. Cf. Sadourny and Laval
(1984) for further details. See also Planetary Boundary Layer and Surface
The formulation of gravity-wave drag closely follows the linear model described
by Boer et al. (1984) . The drag at any level is proportional to the vertical
divergence of the wave momentum stress, which is formulated as the product
of a constant aspect ratio, the local Brunt-Vaisalla frequency, a launching
height determined from the orographic variance over the grid box (see Orography),
the local wind velocity, and its projection on the wind vector at the lowest
model level. The layer where gravity-wave breakdown occurs (due to convective
instability) is determined from the local Froude number; in this critical
layer the wave stress decreases quadratically to zero as a function of
The solar constant is the AMIP-prescribed value of 1365 W/(m2).
The orbital parameters and seasonal insolation distribution are calculated
from PMIP recommendations. A seasonal, but not a diurnal cycle in solar
forcing, is simulated.
Carbon dioxide concentration is prescribed. The values of 345, 280 and
200 ppm for control, 6ka and 21ka, respectively have been chosen for PMIP.
A value of 280 ppm (pre-indusrtrial value) has been chosen for the computed
SSTs control. Three-dimensional ozone concentration is diagnosed as a function
of the 500 hPa geopotential heights following the method of Royer et al.
(1988) . Radiative effects of water vapor and cloud water, but not those
of aerosols, are also included (see Radiation).
Shortwave radiation is derived from an updated scheme of Fouquart and Bonnel
(1980) . Upward/downward shortwave irradiance profiles are evaluated in
two stages. First, a mean photon optical path is calculated for a scattering
atmosphere including clouds and gases. The reflectance and transmittance
of these elements are calculated by, respectively, the delta-Eddington
method (cf. Joseph et al. 1976) and by a simplified two-stream approximation.
The scheme evaluates upward/downward shortwave fluxes for two reference
cases: a conservative atmosphere and a first-guess absorbing atmosphere;
the mean optical path is then computed for each absorbing gas from the
logarithm of the ratio of these reference fluxes. In the second stage,
final upward/downward fluxes are computed for two spectral intervals (0.30-0.68
micron and 0.68-4.0 microns) using more accurate gas transmittances (Rothman
1981) and with adjustments made for the presence of clouds (see Cloud Formation).
For clouds, the asymmetry factor is prescribed, and the optical depth and
single-scattering albedo are functions of cloud liquid water content after
Stephens (1978) .
Longwave radiation is modeled in six spectral intervals between wavenumbers
0 and 2.82 x 105 m-1 from the method of Morcrette
(1990, 1991 ). Absorption by water vapor (in two intervals), by the water
vapor continuum (in two intervals in the atmospheric window, following
Clough et al. 1980) , by the carbon dioxide and the rotational part of
the water vapor spectrum (in one interval), and by ozone (in one interval)
is treated. The temperature and pressure dependence of longwave absorption
by gases is included. 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 from the cloud liquid water path.
For purposes of the radiation calculations, all clouds are assumed to overlap
randomly in the vertical. See also Cloud Formation.
When the temperature lapse rate is conditionally unstable, subgrid-scale
convective condensation takes place. If the air is supersaturated, a moist
convective adjustment after Manabe and Strickler (1964) is carried out:
the temperature profile is adjusted to the previous estimate of the moist
adiabatic lapse rate, with total moist static energy in the column being
held constant. The specific humidity is then set to a saturated profile
for the adjusted temperature lapse, and the excess moisture is rained out
If the temperature lapse rate is conditionally unstable but the air
is unsaturated, condensation also occurs following the Kuo (1965) cumulus
convection scheme, provided there is large-scale moisture convergence.
In this case, the lifting condensation level is assumed to be at the top
of the PBL, and the height of the cumulus cloud is given by the highest
level for which the moist static energy is less than that at the PBL top
(see Planetary Boundary Layer). It is assumed that all the humidity entering
each cloudy layer since the last call of the convective scheme (30 minutes
prior) is pumped into this cloud. The environmental humidity is reduced
accordingly, while the environmental temperature is taken as the grid-scale
value; the cloud temperature and humidity profiles are defined to be those
of a moist adiabat. The fractional area of the convective cloud is obtained
from a suitably normalized, mass-weighted vertical integral (from cloud
bottom to top) of differences between the humidities and temperatures of
the cloud vs those of the environment. As a result of mixing, the environmental
(grid-scale) temperature and humidity profiles evolve to the moist adiabatic
values in proportion to this fractional cloud area, while the excess of
moisture precipitates (see Precipitation). Mixing of momentum also occurs.
There is no explicit simulation of shallow convection, but the moist
convective adjustment produces similar effects in the moisture field (cf.
Le Treut and Li 1991) . See also Cloud Formation.
Cloud cover is prognostically determined, as described by Le Treut and
Li (1991) . Time-dependent cloud liquid water content (LWC) follows a conservation
equation involving rates of water vapor condensation, evaporation of cloud
droplets, and the transformation of small droplets to large precipitating
drops (see Precipitation).
Cloud water content and cloud fraction are interactively used in the
radiative code through the calculation of cloud optical thickness and cloud
emissivity (Le Treut 1994).
The fraction of convective cloud in a grid box is unity if moist convective
adjustment is invoked; otherwise, it is given by the surface fraction of
the active cumulus cloud obtained from the Kuo (1965) scheme (see Convection).
Cloud forms in those layers where there is a decrease in water vapor from
one call of the convective scheme to the next (every 30 minutes), and the
cloud LWC is redistributed in these layers proportional to this decrease.
The fraction of stratiform cloud in any layer is determined from the
probability that the total cloud water (liquid plus vapor) is above the
saturated value. (A uniform probability distribution is assumed with a
prescribed standard deviation--cloud typically begins to form when the
relative humidity exceeds 83 percent of saturation.) This stochastic approach
also crudely simulates the effects of evaporation of cloud droplets. Cf.
Le Treut and Li (1991) for further details. See also Precipitation.
Both convective and large-scale precipitation are linked to cloud LWC (see
Cloud Formation). For warm clouds, the precipitation is parameterized using
the relationship proposed by Sundqvist (1981) in which Cl = 2 x 104
kg per kg and Ct = 5.5 10-4 s-1. For cold clouds,
we use a relationship that takes into account the terminal falling speed
of the crystals as described in Starr and Cox (1985).
Planetary Boundary Layer
The PBL is represented by the first 4 levels above the surface (at sigma
= 0.979, 0.941, 0.873, and 0.770). The PBL top is prescribed to be at the
sigma = 0.770 level; here vertical turbulent eddy fluxes of momentum, heat,
and moisture are assumed to vanish. See also Diffusion, Surface Fluxes,
and Surface Characteristics.
Raw orography obtained at 10 x 10-minute resolution from the U.S. Navy
dataset (cf. Joseph 1980) is area-averaged over the model grid boxes. The
orographic variance about the mean value for each grid box is also computed
from the same dataset for use in the gravity-wave drag parameterization
(see Gravity-wave Drag). The rugosity is also linked to orographic variance.
For LGM, we add the anomaly of Peltier (Peltier(21k) -Peltier(0k)) to
the control run.
For control: SSTs and sea ice extent are specified from the climatolgy
of the 79-88 Reynold?s data (1988) used in AMIP (10 years mean).
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. When points remain ocean, 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. When
points are ocean in summer and ice in winter (for Nothern Hemisphere),
a "trapezoidal" function is used with an ad-hoc duration of the sea-ice
deduced from present-day temperatures equivalent.
For computed SSTs experiments, the AGCM was coupled to a mixed layer
ocean (50m) (Le Treut et al. 1994). The ocean heat transport is daily prescribed
using the present-day diagnosed ocean heat transport accounting for the
sea level drop (105m.) as described in the PMIP Newsletter (Webb et al.
For 0fix and 6fix, sea ice extents are prescribed monthly. The surface
temperature of the ice is predicted from the balance of energy fluxes (see
Surface Fluxes) that includes conduction heating from the ocean below.
This conduction flux is proportional to the difference between the surface
temperature and that of melting ice (271.2 K), and is inversely proportional
to the ice thickness (prescribed to be a uniform 3 m). Snow that accumulates
on sea ice modifies its albedo and thermal properties. See also Snow Cover
and Surface Characteristics.
For 21fix, sea ice edge of CLIMAP data is used for February and August.There
is not any percent of sea ice for each grid box. When the CLIMAP 21k sea-ice
extent is interpolated over the model grid 48x36 if the result is higher
than 50%, we consider that the grid box is sea ice. When the result is
lower than 50%, the whole box is free of ice.
For computed SSTs experiments, the sea ice appears when the temperature
is lower than -2°C and disappear when it is higher.
If the air temperature at the first level above the surface (at sigma =
0.979) is <0 degrees C, precipitation falls as snow. Prognostic snow
mass is determined from a budget equation, with accumulation and melting
over both land and sea ice. Snow cover affects the surface albedo and the
heat capacity of the surface. Sublimation of snow is calculated as part
of the surface evaporative flux, and snowmelt contributes to soil moisture.
See also Surface Characteristics, Surface Fluxes, and Land Surface Processes.
The surface roughness lengths over the continents are prescribed as a function
of orography and vegetation from data of Baumgartner et al. (1977) , and
their seasonal modulation is inferred following Dorman and Sellers (1989)
. Roughness lengths over ice surfaces are a uniform 1 x 10-2
m. Over ocean, the surface drag/transfer coefficients (see Surface Fluxes)
are determined without reference to a roughness length.
Surface albedos for oceans and snow-free sea ice are prescribed from
monthly data of Bartman (1980) , and for snow-free continents from monthly
data of Dorman and Sellers (1989) . When there is snow cover, the surface
albedo is modified according to the parameterization of Chalita and Le
Treut (1994) , which takes account of snow age, the eight designated land
surface types, and spectral range (in visible and near-infrared subintervals).
The longwave emissivity is prescribed as unity (blackbody emission)
for all surfaces.
The surface solar absorption is determined from surface albedos, and longwave
emission from the Planck equation with prescribed emissivity of 0.96 (see
In the lowest atmospheric layer, surface turbulent eddy fluxes of momentum,
heat, and moisture are expressed as bulk formulae multiplied by drag/transfer
coefficients that are functions of wind speed, stability, and (except over
ocean) roughness length (see Surface Characteristics). The transfer coefficient
for the surface moisture flux also depends on the vertical humidity gradient.
Over the oceans, the neutral surface drag/transfer is corrected according
to the local condition of surface winds. For strong surface winds, the
drag/transfer coefficients are determined (without reference to a roughness
length) as functions of surface wind speed and temperature difference between
the ocean and the surface air, following Bunker (1976) . For conditions
of light surface winds over the oceans, functions of Golitzyn and Grachov
(1986) that depend on the surface temperature and humidity gradients are
utilized. In the transition region between these wind regimes, surface
drag/transfer coefficients are calculated as exponential functions of the
surface wind speed.
In addition, the momentum flux is proportional to the wind vector extrapolated
to the surface. The sensible heat flux is proportional to the difference
between the potential temperature at the ground and that extrapolated from
the atmosphere to the surface. The surface moisture flux is proportional
to the potential evaporation (the difference between the saturated specific
humidity at the surface and the extrapolated atmospheric humidity) multiplied
by an evapotranspiration efficiency beta. Over oceans, snow, and ice, beta
is set to unity, while over land it is a function of soil moisture (see
Land Surface Processes).
Above the surface layer, but only within the PBL, turbulent eddy fluxes
are represented as diffusive processes (see Diffusion and Planetary Boundary
Land Surface Processes
Ground temperature and bulk heat capacity (with differentiation for bare
soil, snow, and ice) are defined as mean quantities over a single layer
of thickness about 0.15 m (over which there is significant diurnal variation
of temperature). The temperature prediction equation, which follows Corby
et al. (1976) , includes as forcing the surface heat fluxes (see Surface
Fluxes) and the heat of fusion of snow and ice.
Prognostic soil moisture is represented by a single-layer "bucket" model
after Budyko (1956) , with uniform field capacity 0.15 m. Soil moisture
is increased by both precipitation and snowmelt, and is decreased by surface
evaporation, which is determined from the 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 maximum of unity or twice the ratio of local
soil moisture to the constant field capacity. Runoff occurs implicitly
if the soil moisture exceeds the field capacity. Cf. Laval et al. (1981)
for further details.
Last update November 9, 1998. For further information, contact: Céline