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

PMIP Documentation for LLN_NH_1

Institut d?astronomie et de geophysique G. Lemaître - Universite catholique de Louvain : Model LLN 2-D NH (7 sectors) 1995

PMIP Representative(s)

Dr Marie-France Loutre, Institut d'astronomie et de geophysique G. Lemaitre, Universite catholique de Louvain, Chemin du cyclotron, 2 B-1348 Louvain-la-Neuve, BELGIUM; Phone : +32-10 473299; Fax : +32-10 474722; e-mail :

World Wide Web URL:

Model Designation


LLN 2-D NH (7 sectors) 1995.

Model Identification for PMIP


PMIP run(s)

0cal, 6cal, 21cal

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

Model Lineage

The LLN 2-D NH climate model is the model designed in Louvain-la-Neuve to test the Milankovitch theory

and used for long-term paleoclimate simulations.

Boundary conditions were adapted to fill the PMIP requirement.

Model Documentation

Gallée H., van Ypersele J.P., Fichefet Th., Tricot Ch. and Berger, A., 1991 : Simulation of the last glacial cycle by a coupled, sectorially averaged climate - ices-sheet model. I. The climate model. J. Geophys. Res., 96(D7): 13,139--13,161.

Gallée H., van Ypersele J.P., Fichefet Th., Marsiat I., Tricot Ch. and Berger A., 1992 : Simulation of the last glacial cycle by a coupled, sectorially averaged climate - ice-sheet model. II. Response to insolation and CO2 variation. J. Geophys. Res., 97(D14), 15,713-15,740.

Numerical/Computational Properties

Horizontal Representation

Regular grid in latitude (5deg.) ; sectorially averaged in longitude ; finite differences.

Horizontal Resolution

Resolution in latitude : 5 deg.

Resolution in longitude : 7 'sectors' (i.e. types of surfaces) (sea ice, ice free ocean, snow free land, snow field, Greenland ice sheet, North american ice sheet, Eurasian ice sheet).

Vertical Domain

From 0 to 1000hPa.

Vertical Representation

In pressure coordinates.

Vertical Resolution

The atmospheric dynamics is the classical two-level quasi-geostrophic (QG) system written in pressure coordinates and zonally averaged (Sela and Wiin-Nielsen, 1971; Ohring and Adler, 1978). The QG potential vorticities are computed at level 250 and 750 hPa.

The atmospheric temperature is calculated at 500 hPa and the surface temperature is set equal to the ground temperature.

The eastward components of the wind are computed at 250 and 750, and extrapolated for 1000 hPa.

A Three-layer model is assumed for the solar radiation scheme ((1) from the surface to the cloud bottom, (2) the cloud layer, (3) from the cloud summit to the top of the atmosphere.

For the long-wave radiation the atmosphere is divided into 10 to 15 layers depending on the surface pressure over each surface type.

Computer/Operating System

Machine : HP-Exemplar

OS : SPP/UX - 1 processor.

Computational Performance

20s (CPU) per year.


The model is run until a quasi-equilibrium is established, i.e. after 50 years of integration.

Time Integration Scheme(s)

The QG potential vorticity is solved by the way of a Cranck-Nicholson (semi-implicit) scheme. The tridiagonal method of Richtmayer (1957) is used both in the integration of the QG prognostic equations and in the solution of the thermal stream function. The time step used for QG dynamics is 1 day. The surface fluxes and the diabatic heating are reevaluated every three days.


Input values are linearly interpolated.

Sampling Frequency

Results are those of the last year of integration.

Dynamical/Physical Properties

Atmospheric Dynamics

2 level quasi-geostrophic model.

The dependent variables are the QG potential vorticity at 250 and 750 hPa.

The atmospheric temperature is calculated at 500 hPa and the surface temperature is set equal to the ground temperature.

The eastward components of the wind are computed at 250 and 750, and extrapolated for 1000 hPa.


The transport of QG potential vorticity is computed following a diffusion law. The coefficient is a function of the latitude and is different for the two levels. In the equatorial regions, the parameterisation adopted is the one suggested by Branscome (1983). Outside the tropics, the diagnostic values of Fuenzalida (1973) are used. The diffusion coefficients are re-evaluated every day.

Gravity-wave Drag

Not represented.

Solar Constant/Cycles

Solar constant according to PMIP specification : 1365 W/m**2.

The orbital parameters and seasonal insolation distribution are calculated after PMIP recommendations. Seasonal cycle is simulated.


The total ozone amount is parameterised as a function of latitude and time following Van Heuklon (1979) and the vertical profiles are given by the formulation of Lacis and Hansen (1974).

Three kinds of aerosols are considered : continental, maritime and unperturbed stratospheric. Their vertical profiles and radiative parameters are taken from World Meteorological Organization (WMO, 1986).

Carbon dioxide concentration (PMIP-prescribed values) : 280ppmv (0 and 6 kyrBP), 200ppmv (21 kyrBP)

Radiatively active gases :



The solar radiation scheme is that of Tricot and Berger (1988) with slight modifications. The following processes are taken into account : absorption by H2O, CO2 and O3, Rayleigh scattering, absorption and scattering by cloud and aerosol layers, and reflection by the surface. A three-layer model of the atmosphere is assumed. The first layer extends from the surface to the pressure level of the cloud bottom, the second layer is filled up with a zonally averaged effective cloud and the third layer extends from the cloud summit to the top of the atmosphere.

Three similar layers are defined for the clear fraction of the sky. O3 is assumed to be present only in the upper layer, while H2O, CO2 and aerosols are distributed within each layer.

One single spectral interval (0.25 - 4 micron) is considered to compute solar radiative fluxes. We use the parameterisations given by Fouquart (1986) for the spectrally averaged optical thickness for Rayleigh scattering, and for the spectrally averaged gaseous transmittances which are used for each layer as a function of pressure depths and gaseous amounts. Unlike the method discussed by Tricot and Berger (1988) the transmittance and reflectance for each homogeneous layer are obtained with the delta-Eddington approximation for the cloud layer only, whereas first-order expansions of the analytical solutions of the two-stream formulation are used in the present work for other layers. Then the three layers are combined using the ascending method presented by Fouquart and Bonnel (1980) and extended by Fouquart (1986). In the present scheme, downward radiation in each layer is treated as a collimated radiation with the same zenith angle as the incident radiation at the top of the atmosphere, in the case of a clear sky. For a cloudy sky, the same assumption is made for the upper layer, while the downward radiation into and below the cloud is assumed to be diffuse, the same being assumed for the upward radiation in both cloudy and clear conditions. To take into account partial cloudiness, the fluxes are first calculated separately for a clear and a cloudy sky and then weighted linearly by the clear and cloudy fractions.

The longwave radiation scheme follows a wide-band formulation of the radiative transfer equation which was designed for use in GCMs (Morcrette, 1984; Morcrette et al., 1986). The upward flux at the top of the atmosphere and the downward flux at the surface are computed by dividing the model atmosphere into 10--15 layers, the exact number depending on the surface pressure over each surface type. A balance between flux accuracy and computer burden is achieved by dividing the longwave spectrum into four wide wavenumber intervals : [0--50000m-1 + 125000--282000m-1], [50000--80000m-1], [80000--125000m-1] and [97000--111000m-1]. (Morcrette, 1984) The line absorptions by H2O, CO2 and O3 are separately treated for each gas. The absorption by the H2O continuum is parameterised following the results of Clough et al. (1980) and its dependency on temperature in the two wide bands between 50000 and 125000 m-1 is taken into account according to the formulation of Roberts et al. (1976).

Morcrette (1984) has presented fast parameterisations of the spectrally averaged transmissivity functions for each absorber based on results from a detailed narrow-band model. In each spectral band, the effects of overlapping between gaseous absorbers is treated by making the usual assumption that the transmissivities may be multiplied together to give the overlapped values. The modification of the clear-sky fluxes due to the cloud layer is introduced following the method discussed by Washington and Williamson (1977) and the cloud cover is supposed, for simplicity, to behave as a blackbody, as is the surface of the Earth. The first step is to calculate the fluxes corresponding to a clear-sky atmosphere between the different levels along the vertical. The fluxes are then re-evaluated assuming an overcast cloud layer of unity emissivity. Finally, the actual fluxes are derived from a linear combination of the fluxes obtained for the clear and overcast atmospheres.


Cloud Formation

An effective single cloud is prescribed in each latitude belt, with annual mean cloud amounts deduced from London's (1957) data (Ohring and Adler, 1978) and supposed to be equal over each surface type. Monthly variations of the zonal cloudiness are introduced by superimposing a seasonal cycle on the annual mean values following the monthly mean data of Berlyand and Strokina (1980). The base and top altitudes of each cloud layer and its optical thickness are kept fixed throughout the year following Chou et al. (1981).


Precipitation is calculated from the observed seasonal cycle of the zonally averaged precipitation of Jaeger (1976). Over the ice sheets the observations are first modified in order to take into account the sectorially averaged surface slope, the sectorially averaged elevation and the sectorially averaged continentality of the ice sheet in each latitudinal belt (Oerlemans, 1982). Then a second correction is achieved for each surface type by the constraint for the total latent heat release due to water condensation to be equal to the total surface water vapor production (from evaporation or sublimation processes) at each time step.

Planetary Boundary Layer

Not resolved.

Orography/Land-Sea Mask

Topography and ice cover from Peltier data, averaged over the different sector and over 5 deg. zone in latitude.


Computed SSTs and sea ice : parametrisation of the zonal equator-to-pole (meridional) oceanic heat transport (OHT) (Sellers, 1973). The upper-ocean model is represented by the variable depth and temperature mixed-layer model of Gaspar (1988).

Ocean salinity is not represented in this model.

Sea Ice

Computed SSTs and sea ice. The upper-ocean model is represented by the variable depth and temperature mixed-layer model of Gaspar (1988). The sea ice model is a thermodynamic model which calculates ice thickness and leads percentage. Ice dynamics is ignored. A simple parameterisation of sea ice formation and decay in warm water (above the freezing temperature) is introduced.

Snow Cover

In the present model a snow mass budget is computed for land and ice sheet in each latitudinal band and is a function of snowfall and snowmelt. Over the ice sheet, the snow depth is assumed to be uniform and snow is melted until it disappears. Then the excess of energy is used to melt continental ice. Over land, the melting is partitioned between a decrease in snow extent and a decrease in snow depth, so that snow area gradually decreases during the snow melting period. The fractional area of land over which precipitation falls as snow is parameterised as in Harvey (1988). The proportion of precipitation over an ice sheet that falls as snow rather than rain is parameterised according to Ledley (1983).

Surface Characteristics

In each latitudinal belt the surface is divided into at most seven oceanic or continental surface types (sea ice, ice free ocean, snow free land, snow field, Greenland ice sheet, North American ice sheet, Eurasian ice sheet). Land surface albedo is chosen to be a function of water availability, following Saltzman and Ashe (1976). The albedo of sea water is function of sea ice surface temperature and of solar zenith angle (Briegleb and Ramanathan; 1982). Over land areas not covered by trees and over the ice sheet, the snow surface albedo is a function of snow precipitation frequency and snow surface temperature (Danard et al.; 1984). Corrections for snow depth and solar zenith angle are also applied. The shifts of the taiga/tundra boundary modify the surface albedo, because of a larger snow albedo over tundra than over taiga. These shifts are forced by a variation of surface temperature (Otterman et al., 1984; Harvey, 1988).

Surface Fluxes

The sensible heat flux is parameterised following Saltzman-Ashe (1976). It takes into account synoptic and diurnal temperature variances and the wind speed at the surface.

The latent heat flux, which corresponds to the evaporation of surface water, is linearly related to the sensible heat flux and is multiplied by a water availability factor (Saltzman 1980). Over snow or continental ice, this factor is a function of the surface temperature. Over snow-free land, it is function of the mean moisture.

Radiative fluxes are taken into account as described in the related sections.

Land Surface Processes

The variation in the mean moisture in the surface layer is a balance between the precipitation rate and the rate of evaporation from land surface. A minimum and a maximum value of the mean moisture are assumed (i.e. runoff is implicitly taken into account).

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