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

PMIP Documentation for CNRM-2

Centre National de Recherches Météorologiques: Model CNRM ARPEGE Cy14c (T31 L19) 1996

PMIP Representative(s)

Mr. Michel Deque, CNRM METEO-FRANCE, 42 avenue coriolis, 31057 TOULOUSE cedex, FRANCE; Phone: 33 561 07 93 82; Fax:33 561 07 96 10; e-mail:


Mr. Jean-Francois Royer, CNRM METEO-FRANCE, 42 avenue coriolis, 31057 TOULOUSE cedex, FRANCE; Phone: 33 561 07 93 77; Fax:33 561 07 96 10; e-mail:

World Wide Web URL:

Model Designation

CNRM ARPEGE Cy14c (T31 L19) 1996

Model Identification for PMIP


PMIP run(s)

0fix, 6fix.

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

but the 1rst, 5th and 9th year are leap years (29 days in february)

Model Lineage

The model is a version (Cycle 14) of the ARPEGE climate model, which is designed for use by the French climate community. (ARPEGE is the acronym for Action de Recherche Petite Echelle Grande Echelle: Research Project on Small and Large Scales.) ARPEGE is the successor to the EMERAUDE baseline model, and differs substantially in its numerics (especially the horizontal representation), as well as in the formulation of horizontal diffusion, gravity-wave drag, cloud formation, and land surface schemes. Some modifications that have been introduced in this version 2 of Arpege-climat compared to version 1 (cycle 11).

Model Documentation

main reference:

There is yet no paper based on the PMIP version or version 2 of ARPEGE Climat.

As a reference paper you can quote a published paper describing the results of version 0:

Deque,M; Dreveton,C; Braun,A; Cariolle,D (1994): The ARPEGE/IFS atmosphere model: A contribution to the French community climate modelling. Climate Dyn. 10(4/5, Sep), 249-266.

secondary reference(s):

If required for detailed description of the parameterizations of version 2 an unpublished technical document: "Community modelling, 1996: ARPEGE-CLIMAT - Version 2- I. Algorithmic documentation." is available from:

Dandin, Ph; J.J. Morcrette,(1996): The ECMWF FMR scheme in the Meteo-France Climate model Arpege-Climat. Note de travail du Groupe de meteorologie de grande Echelle et Climat (GMGEC), no 50.( Available from: CNRM, 42 Av G Coriolis, F-31057 Toulouse Cedex, France).

Deque,M; Piedelievre,JP (1995): High resolution climate simulation over Europe. Climate Dyn. 11(6, Aug), 321-339.

Douville, H.; J.F. Royer (1995): A new snow parameterization for the Meteo-France climate model. Part I: validation in stand-alone experiments. Climate Dynamics, 12, 21-35

Douville, H.; J.F. Royer (1995): A new snow parameterization for the Meteo-France climate model. Part II: validation in a 3-D GCM experiment. Climate Dynamics, 12, 37-52

Geleyn,JF; Bazile,E; Bougeault,P; Deque,M; Ivanovici,V; Joly,A; Labbe,L; Piedelievre,JP; Piriou,JM; Royer,JF (1995): Atmospheric parameterization schemes in Meteo-France's Arpege NWP model. ECMWF Seminar Proc. Parametrization of sub-grid scale physical processes(5-9 September 1994), 385-402.

Land surface processes are represented by the scheme of Noihlan and Planton (1989) which was implemented in the ARPEGE model by JF Mahfouf:

Mahfouf,JF; Manzi,AO; Noilhan,J; Giordani,H; Deque,M (1995): The land surface scheme ISBA within the Meteo-France climate model ARPEGE. Part I: Implementation and preliminary results. J. Climate 8(8, Aug), 2039-2057.

Manzi, AO and S Planton (1994): Implementation of the ISBA parameterization scheme for land surface processes in a GCM: an annual cycle experiment. J. Hydrology, 155, 355-389 Noilhan, J and JF Mahfouf (1996): The ISBA land surface parameterization scheme. Glob. Planet. Change, 13, 145-159

Sensitivity experiments are described in:

Douville H and JF Royer (1996): Sensitivity of the Asian summer monsoon to an anomalous Eurasian snow cover within the Meteo-France GCM. Climate Dynamics, 12, 449-466

Douville H and JF Royer (1997): Influence of the temperate and boreal forests on the Northern hemisphere climate in the Meteo-france climate model. Climate Dynamics, 13, 57-74

Numerical/Computational Properties

Horizontal Representation

Formulation of horizontal variation of model variables spectral representation on spherical harmonics basis functions with transformation to a Gaussian grid for calculation of nonlinear quantities and physics. All prognostic variables including water vapor specific humidity are spectrally transformed.

Horizontal Resolution

triangular truncation T31 resolution of the Gaussian grid: about 3.75 degrees In this simulation with Cy14c (as in all simulations from ARPEGE Cy11 onwards ) the Gaussian grid is reduced longitudinally near the poles so that its horizontal resolution is everywhere approximately the same (cf. Hortal and Simmons, 1991; Coutier and Naughton, 1994).

dim_longitude*dim_latitude: 96*48.

Vertical Domain

The upper layer is from 0 to 20 hPa with the midlevel corresponding to 7.4 hPa

For a surface pressure of 1000 hPa, the bottom layer is from 990 to 1000 hPa with the midlevel at 995 hPa.

Vertical Representation

Finite difference scheme in hybrid vertical coordinates according to Simmons and Burridge (1981)

Relevant references.

Simmons A and D Burridge (1981): An energy and angular momentum conserving vertical finite difference scheme and hybrid vertical coordinates. Mon. Wea. Rev, 109, 758-766

Vertical Resolution

There are 19 layers unevenly spaced hybrid sigma-pressure levels.

For a surface pressure of 1000 hPa, there are 5 layers below 800 hPa 6 layers above 200 hPa

The 19 vertical levels in hybrid coordinate are defined by: Pressure p(L) at half levels L (i.e. at the interface between atmospheric layers) are computed from the surface pressure p_s from the formula:

p(L) = A(L)+B(L)* ps

where A(L) and B(L) are coefficients which in this T31-L19 simulation

have been given the values:

A:2000.000,4000.000,6491.873,10000.000,13466.847,15602.481,16578.893,16568.072,15742.010,14272.696,12332.122, 10092.277,7725.153,5402.739,3297.026,1580.005,423.666,0.000,0.000


Computer/Operating System

the PMIP simulation was run on Cray J916 in UNICOS environment.

Computational Performance

For the PMIP experiment, approximately 5 minutes monoprocessor J916 computation time per simulated day.


The model was started from the 31 May 1988 situation extracted from a previous T42 AMIP simulation after interpolation on a T31 grid.

Time Integration Scheme(s)

semi-implicit centered (leap-frog) scheme for the dynamics Asselin time filter with a coefficient of 0.05 Restart of an integration with a direct uncentered time step every 5 day in order to remove the computational mode The time step for dynamics and physics is 30 mn Radiative fluxes are recomputed every 3 hours and kept constant in the intermediate time steps

Relevant references.

Asselin (1972): Frequency filter for time integrations. Mon. Weath. Rev., 100, 487-490


The horizontal representation of the ARPEGE model makes use of a vertical moisture borrowing scheme to correct spurious negative humidities produced by the Gibbs effect of the spectral truncation

Sampling Frequency

For the PMIP simulation, the history of selected variables is written every 6 hours, with a full model history saved at 5-day intervals.

Dynamical/Physical Properties

Atmospheric Dynamics

As in the baseline model, primitive-equation dynamics are expressed in terms of vorticity and divergence, temperature, and specific humidity, and the natural logarithm of surface pressure (or, on option, the surface pressure itself). Conservation of mass is enforced now for dry air (instead of total mass formerly) by a uniform correction of surface pressure at each restart of the model (usually every 5 day)


A linear Ñ6 formulation of horizontal diffusion of vorticity, divergence, temperature and specific humidity on constant hybrid sigma-pressure vertical surfaces replaces the Ñ4 scheme of the baseline model. Below 100 hPa, the diffusivity K increases as (n/N)6, where n and N are the meridional and truncation wavenumbers, respectively (i.e., n <= N = 31 at T31 resolution). Above 100 hPa, K increases as the inverse of the pressure, yielding very strong diffusion in the model stratosphere. The mesospheric drag is not applied in this version of Arpege-climat.

Gravity-wave Drag

The parameterization of gravity-wave drag is different from that of the baseline model.

The gravity-wave stress is assumed to be maximum at the surface, and in a direction that depends on the (2 x 2) covariance matrix of the unresolved oragraphy (a measure of the anisotropy of the mountains) as well as on the mean wind in the planetary boundary layer. The magnitude of the surface stress is proportional to the product of the air density, the Brunt-Vaisalla frequency at the surface, the wind speed at the lowest vertical level, and the root-mean square of the unresolved orography in the direction of this wind (calculated from the subgrid-scale orographic variance).

At levels above the surface, the gravity-wave stress is assumed to be in the same direction as the surface stress. The vertically propagating gravity waves do not interact with the mean flow below a critical level of resonance. Above this level, their resonant amplification follows the experimental results of Clark and Peltier (1984), while dissipation proportional to the square of the Froude number also operates. The gravity waves are trapped and reflected with dissipation at the vertical level where the Brunt-Vaisalla frequency becomes zero. Cf. Deque et al. (1994))for further details.

A modification of the parameterization of momentum deposition by gravity-waves in order to increase momentum deposition in the lower atmospheric layers has been introduced by Geleyn et al, 1995. (A parameterization of gravity-wave source due to convection has been introduced in version 2, but is not used in this PMIP simulation).

Radiative Boundary Conditions

Solar constant is the PMIP-prescribed value of 1365 W/m2. 1365. The orbital parameters and seasonal insolation distribution are calculated according to PMIP recommendations.

Seasonal and diurnal cycles are included



Update of the coefficients for the ozone parameterization by recomputing them from a more recent bidimensional photochemical model. Ozone is not treated as a prognostic variable. Enumeration of radiatively active gases and concentration for control, 6ka and 21ka (CO2, ozone, water vapor, other trace gases) and aerosols.


CO2 : 353 ppmv

CH4: 1.72 ppmv

N2O: 310 ppbv

CFC-11: 280 pptv

CFC-12: 484 pptv


CO2 : 286 ppm

CH4: 1.72 ppmv

N2O: 310 ppbv

CFC-11: 280 pptv

CFC-12: 484 pptv

There are a globally uniform distribution for the trace gases (CO2,CH4,N2O, CFC) and a globally uniform vertical profile for ozone imposed geographical distribution for 5 types of aerosols water vapour is treated spectrally by an Eulerian transport scheme.


This PMIP version makes use of the Fouquart and Morcrette Radiation scheme (used at ECMWF) which has been implemented in version 2 of ARPEGE-Climate by P. Dandin (Dandin and Morcrette, 1996). The radiative computations are performed at 3 hour intervals, and the radiative fluxes are kept constant in the intermediate time-steps. Two-stream formulation with photon path distribution method (Fouquart and Bonnel, 1980) in two spectral intervals (0.25-0.68 and 0.68-4.0 microns). Broadband flux emissivity methods, with temperature and pressure dependence of absorption (Morcrette et al, 1986). Absorption coefficients are fitted from AGFL 1982)

6 spectral bands covering 0-282000 m-1:

0- 35000 m-1 + 145000-188000 m-1

50000- 80000 m-1

80000- 97000 m-1 + 111000-125000 m-1

97000-111000 m-1

35000- 50000 m-1

125000-145000 m-1 + 18800-282000 m-1

Droplet absorption and scattering are taken into account and are determined from Liquid and Ice water paths computed from a diagnostic liquid and ice concentration in clouds. A Delta-Eddington method is used for the shortwave and an emissivity formulation for the longwave. A maximum overlap between cloud layers has been assumed

Aerosol effects are taken into account for 5 types of aerosols (Tanré et al, 1984) Four types (desertic, land sea and urban) have a specified geographical distribution, while a well mixed distribution is assumed for the background aerosol. In the shortwave, aerosol scattering and absorption are computed from their Mie parameters. In the longwave, absorption effects of aerosols are based on an emissivity formulation.


The mass flux scheme of Bougeault (1985) used in the basline model has been modified by the Introduction of a vertical variation of the entrainment coefficient used to compute the convective cloud profile The entrainment has been increased in the lower layers. Evaporation of convective precipitation has been parameterized. A distinction between ice and liquid phases for falling precipitation has been introduced (Geleyn et al, 1995)

Shallow convection is parameterized as part of the stability dependent computation of the turbulent exchange coefficients by a modification of the Richardson number using the vertical gradients of specific humidity (Geleyn, 1987)

Cloud Formation

Cloud formation is by the same diagnostic method as is used in the baseline model. The stratiform cloud cover C in each vertical layer is a quadratic function of the relative humidity H = q/qsat in excess of a given critical humidity threshold Hc. C = (max (0, H - Hc)/(1 -Hc))**2.

The critical humidity profile Hc is specified as a function of the vertical hybrid sigma-pressure coordinate at full levels by the following formula: Hc = 1 - HUCOE * z * (1 - z) * (1 + sqrt(HUTIL) * (z - 0.5)).

The 2 empirical coefficients HUCOE and HUTIL of the critical profile are tuned to accomplish several objectives: to produce generally larger cloud fractions than in the baseline model; to yield a planetary albedo of about 0.30; to result in approximate balance of global annual-mean radiation at the top of the atmosphere. Cf. Deque and Piedelievre (1995)and Deque et al. (1994) for further details. In the T31 version used for this PMIP simulation we have taken HUCOE = 1.2 and HUTIL = 3. A convective cloud cover is computed from the total convective precipitation and combined with the stratiform cloud cover


Diagnostic method based on large-scale saturation: for each vertical layer at each grid point the water in excess of the saturation specific humidity is removed and precipitated in the underlying layer, where it can evaporate if the layer in not saturated. Both liquid and solid precipitation is taken into accountre

potential instability on the vertical and positive water vapor convergence over the depth of the convective cloud as computed by the Bougeault (1985) parameterization. Large scale precipitation is evaporated and melted in the underlying layers a subgridscale distribution is assumed only for the interception of convective precipitation by the vegetation canopy

Planetary Boundary Layer

Stability dependent vertical turbulent exchange coefficient specified as a function of the Richardson number according to the formulation of Louis et al (1981) with some modifications, and a specified variation of the mixing length with height. No diagnosis of PBL top.


AMIP data set (U.S. Navy 10' x 10' data set) is used. Spectral truncation with a constraint to reduce Gibbs oscillations over the oceans is the method of tranforming raw data to model grid, including smoothing/filtering procedures. Global-average value of model orography is 175 m (538 m over the continents)

A fractional land cover for each grid point on the Gaussian grid was computed from the U.S. Navy data set. Points with a fractional land cover above 50 % were specified as land.


For control: PMIP datasets will be used for SSTs and sea-ice have been prepared at PCMDI and were calculated by averaging the 10-year AMIP datasets (1979-1988).

For 6 fix: SSTs and sea-ice prescribed at their present day value, as in the control run.

Sea Ice

PMIP monthly sea ice extents are prescribed based on National Oceanic and Atmospheric Administration Ice Joint Center. The surface temperature of the ice is determined from a balance of energy fluxes (see Surface Fluxes) that includes conduction from the ocean below. The conduction flux is obtained by the Deardorff (1978) force-restore method, where the restore temperature is the ice melting point and the thermal inertia is modified from that used over land surfaces. Accumulated snow modifies the albedo, but not the thermal properties of the ice. See also Snow Cover and Surface Characteristics.

Snow Cover

The introduction of a more physically based snow cover parameterization in which the albedo and density of the snow cover is a function of time since the last snowfall is described in Douville and Royer, 1995.

The land surface scheme includes the evolution of a single snow layer represented by the snow mass Wn The prognostic equation for Wn takes into account snow precipitation Pn, snow sublimation En, and snowmelt Fn:

d(Wn)/dt = Pn - En - Fn.

Snowmelt occurs when the surface temperature is larger than 273.16 K. The snow density increases exponentially with an e-folding time of about 4 day since the last snowfall from a value of 100 kg/m3 for fresh snow up to a maximum value of 300 kg/m3 for old snow. Fractional snow cover for radiative computations Pnr is computed as a function of the snow mass Wn and a critical value Wnc by the following formula:

Pn = Wn / (Wn + Wnc * (z0 / z0cp))

This fractional snow cover is used for computing surface albedo and emissivity

surface roughness is modified by the presence of snow.

z0 = z0f * (1-Pnz) + zon * Pnz

where z0f is surface roughness in the absence of snow z0n is roughness length of snow (set to 0.001 m) Pnz is a fractional cover for roughness computed as

Pnz = Wn / (Wn + Wnc * (1+ z0f/z0cz) )

where z0cz is a characteristic roughness length (0.025 m), Wn a characteristic snow mass ( 10 km/m2).

Surface Characteristics

Surface characteristics of land are different from those of the baseline model (cf. Mahfouf et al. (1995).

Primary and secondary vegetation cover and type specified in each grid box from 13 classes determined by the Manzi and Planton (1994) simplification of the Wilson and Henderson-Sellers (1985) dataset. Roughness length, leaf area index (LAI), and minimum stomatal resistance required by the land surface scheme are specified according to the vegetation class and soil characteristics of the grid box by blending values associated with the primary and secondary vegetation weighted in a 3 to 1 ratio, respectively. The coverage and roughness length of deciduous and cultivated vegetation also undergo a seasonal cycle. The roughness length includes a contribution from local subgrid-scale orography as well.

It does not really include different soil types but rather use a specification of soil properties such as albedo and percentage of silt, clay and sand from 3 categories of soil color, soil texture and drainage as explained in Mahfouf et al (1995). Soil color, column depth (derived from drainage data of Wilson and Henderson-Sellers (1985)), and texture (fraction of sand and clay required by the land land surface scheme obtained from Webb et al. (1991) data) also are specified for each grid box. The depth of the active soil layer required by the land surface scheme is assigned according to the larger of the vegetation root depth vs the bare-soil depth.

As in the baseline model, surface albedos are a function of solar zenith angle, but their values are assigned differently. Albedos over land are prescribed according to vegetation cover (blended as described above for primary and secondary vegetation types) and soil color. Observed albedo from Meteosat and NOAA satellites are used to correct the bareground albedo when the vegetation cover is less than 50% Surface longwave emissivities are prescribed in the same manner as in the baseline model.

Surface Fluxes

The model follows the Louis et al. (1981) formulation of turbulent surface fluxes, as in the baseline model, but with roughness lengths for calculation of the momentum flux over land set 10 times larger than those for the heat flux.

Over land, the surface moisture flux is made up of evaporation from bare ground, snow sublimation, and from moisture intercepted by the vegetation canopy, as well as from transpiration by the foliage according to formulations of the land surface scheme.

Land Surface Processes

Land surface processes are simulated by the Interactions between Soil-Biosphere-Atmosphere (ISBA) scheme of Noilhan and Planton (1989) as implemented in the ARPEGE model by Mahfouf et al. (1995). Use of ISBA results in less extreme ground temperatures over the summer continents than in the baseline model.

An improvement of the ISBA soil-vegetation scheme has been achieved by the introduction of a deep drainage term, interception of convective precipitation by the vegetation canopy, and use of a thermal roughness length smaller than the dynamic roughness length. A new climatology for surface albedo and thermo-hydric properties of bare ground has been used. (Mahfouf et al, 1995)

The recent version of the ISBA scheme is described in Noilhan and Mahfouf (1996). It includes 5 prognostic variables: surface temperature, mean surface temperature, surface volumetric water content, mean volumetric water content, and the water amount intercepted by the vegetation canopy. The time dependence of the prognostic variables are formulated as force-restore equations after Deardorff (1978).

ISBA also requires 7 parameters that are prescribed or derived from other surface characteristics: the vegetation cover, leaf area index (LAI), minimum stomatal resistance, surface shortwave albedo, longwave emissivity, active soil depth, and surface roughness length (see Surface Characteristics). In addition, climatoligical/equilibrium temperatures and volumetric water contents, the maximum moisture capacity of the vegetation canopy, as well as transfer coefficients and restoring time constants are specified in the prognostic equations.

8 land surface characteristics:

- soil depth including the root zone: d_veg

- albedo: \alpha

- roughness length: Z_ov

- minimal vegetation density: veg_min (minimum during the seasonal cycle)

- maximal vegetation density: veg_max (maximum "")

- minimal leaf aera index : LAI_min

- maximal Leaf Area Index : LAI_max

- minimum surface resitance to transpiration: Rs_min

The prescribed values for the 13 different vegetation types are given in table 2 of Mahfouf et al (1995)

Elimination of the relaxation of deep soil temperature toward climatological values, by using a 4-layer scheme for solving the heat layer equation in the ground. Introduction of a limitation on the maximum value of the Richardson number in order to maintain a minimum turbulent transfer between the atmosphere and the surface in very stable stratifications.

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