The model evolution has been documented in a number of papers which include publications by Sadourny and Laval (1984), Laval et al. (1981), Le Treut and Li (1991), Ducoudre et al. (1993), Le Treut et al. (1994), and Grenier (1997). Computational aspects are described by Butel (1991). A description of the version 5.2 is also available in the AMIP documentation. The version 5.3 is different from this AMIP version mainly because of the treatment of the surface conditions. It is described by Harzallah and Sadourny (1995) wo also describe the simulated climatology. The standard resolution of the model is also higher. The differences with the version of the model used at UCL also concern the treatment of the turbulent fluxes over mixed surface conditions.The two model versions are otherwise very similar.

Finite differences on a uniform-area, staggered C-grid (e.g., Arakawa and Lamb, 1977), with points equally spaced in sine of latitude and in longitude.

The most usual grid used for the resolution of the equations includes 72 grid points equally spaced in the sine of latitude and 96 points equally spaced in longitude. The lower resolution version of the model (50 x 64 grid-points) is retained for some applications and developments (in particular at LMCE)

Surface to about 3 hPa. But the resolution within the stratosphere is poor ( 4-5 levels out of 15, or 3-4 levels out of 11, depending on model resolution).

Finite-difference sigma coordinate.

There are 15 unevenly spaced sigma levels. The use of a lower vertical resolution - 11 levels - is also retained for some applications.

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 (Sadourny, 1975a,b). Total energy is also conserved for irrotational flow (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.

IPSL-at: Diffusion

Linear horizontal diffusion is applied on constant-pressure surfaces to potential enthalpy, divergence, and rotational wind via a biharmonic operator del(del*del*)del, where del denotes a first-order difference on the model grid, while del* 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, no further horizontal diffusion of specific humidity is included. See 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 following Smagorinsky et al. (1965). Estimation of TKE involves calculation of a countergradient term after Deardorff (1966) and comparison of the bulk Richardson number with a critical value. See Sadourny and Laval (1984) for further details. See also (h).

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-Vaisala frequency, a launching height determined from the orographic variance over the grid box, 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 height.

The carbon dioxide concentration is taken equal to 345 ppmv. The ozone concentration is prescribed as a function of pressure, latitude, and time. A new release of the Morcrette's radiation code, used in the coupled simulations, makes it now possible to prescribe aerosol concentration and other trace gases concentration.

Shortwave radiation is modelled after 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 (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.25-0.68 micron and 0.68-4.0 microns) using more exact gas transmittances (Rothman, 1981) and with adjustments made for the presence of clouds. 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). The diurnal cycle of insolation is not resolved, and the solar constant is taken as 1367 W m^-2.

Longwave radiation is modelled in six spectral intervals between wavenumbers 0 and 2.82 x 10⁵ m^-1 after 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 grey bodies 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 modelled 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, except in the case of contiguous cloud layers, where maximum overlap is preferred.

The effective radius of cloud droplets used to determine the cloud optical depth is prescribed as 10 x 10^-6 m for warm clouds and as 40 x 10^-6 m for cold clouds. The absorption coefficient used in the cloud-emissivity formulation is also different form warm and cold clouds: 130 m² kg^-1 and 50 m² kg^-1, respectively. In each grid box, warm clouds constitute a fraction X and cold clouds a fraction 1-X of the total cloud cover. X is taken equal to 0 for temperatures below 258.15 K and 1 for temperatures above 273.15K, and we assume it varies linearly between these values at intermediate temperatures. (These coefficients have been slightly modified between the CMIP1 and CMIP2 experiments)

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 rate, and the excess of moisture contributes to the liquid water content (LWC) of clouds (see (f)).

If the temperature lapse rate is conditionally unstable but the air is unsaturated, condensation also occurs following Kuo's (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. It is assumed that all the humidity entering each cloudy layer since the last call of the convective scheme 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 adiabatic. 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 contributes to the cloud LWC (see (f)).

There is no explicit simulation of shallow convection, but the moist convective adjustment produces similar effects in the moisture field (Le Treut and Li, 1991). See also (f).

Cloud cover is prognostically determined, as described by Le Treut and Li (1991). Time-dependent cloud 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. The LWC also determines cloud cover and cloud optical properties (see (d)).

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 Kuo's (1965) scheme (see (e)). Cloud forms in those layers where there is a decrease in water vapor, and the cloud LWC is redistributed in these layers proportionally 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 (varying along the vertical) ; cloud typically begins to form when the relative humidity exceeds 83% of saturation. This stochastic approach also crudely simulates the effects of evaporation of cloud droplets. See Le Treut and Li (1991) for further details. See also (g).

For warm clouds, the precipitation rate is parmeterized after Sundqvist (1981) as the product of a characteristic precipitation time scale, T (D 5.5 x 10^-4 s^-1), the prognostic cloud liquid mixing ratio, M, and an exponential function of (M/C)², where C is a prescribed precipitation threshold value taken equal to 2 x 10^-4 kg kg^-1. For cold clouds, the precipitation rate is determined by the ratio of M to a different time scale T' D Z/V, where Z is the depth of vertical layer and V is the terminal velocity of the water droplets, which is determined as an empirical function of M after Heymsfield and Donner (1990). Evaporation of falling convective and large-scale precipitation is not explicitly modelled, but evaporation of small stratiform cloud droplets making up the LWC is simulated stochastically.

Vertical mixing can occur throughout the atmosphere, but is of course more active in the unstable lower layers of the model. The surface drag is parameterized according to Louis (1979) (see below). The mixing coefficients within the atmosphere are computed as the product of a mixing length (which is slowly decreasing with altitude) and a velocity lenght diagnosed as a function of a Richardson number (using a diagnostic equation to determine a turbulent kinetic energy) (see Laval and Sadourny, 1984) Mixed sea-ice and free ocean conditions are possible within a single grid. Surface fluxes are computed saparately over the free ocean and the sea-ice. In the coupled model no attempt is made to use a moist Richardson number.

If the air temperature at the first level above the surface (at sigma D 0.979) is lower than 273.15 K, 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 (l)).

For each continental grid box, eight coexisting land surface types are specified from aggregation of the data of Matthews (1983, 1984): bare soil, tundra, grassland, grassland with shrub cover, grassland with tree cover, deciduous forest, evergreen forest, and rainforest. A fractional area of the grid box, as well as some prescribed canopy characteristics, are associated with each of these eight land surface types.

The surface roughness lengths over the continent 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.1 x 10^-3 m. Over ocean, roughness lengths are evaluated as a function of wind speed according to Charnock (1955).

Over the land free of snow and ice, the surface albedos are derived from the monthly data of Dorman and Sellers (1989). Over ocean, the formulation of Larson and Barkstrom (1977) is used. The albedo of snow-covered surfaces is computed following Chalita and Le Treut (1994) as a function of the fractional area of snow cover, the snow age, the vegetation type, and the spectral range.

The longwave emissivity is set equal to 0.96 for all surface types.

The net radiative flux at the surface and the turbulent fluxes of momentum, sensible heat, and moisture are calculated separately for the continental, oceanic, and sea-ice portions of each grid box (see (h)). It is worth pointing out that bare soil and vegetation are treated as a single medium for calculations of the surface radiative budget and the turbulent flux of momentum and sensible heat over the land area.

The surface turbulent fluxes of momentum, sensible heat, and moisture are computed by bulk formulas using drag/transfer coefficients that are functions of wind speed, stability, and roughness length. The drag/transfer coefficient for the surface moisture flux also depends on the vertical humidity gradient. (Over land, the formulation of the moisture flux is somewhat more elaborated; see below.) Over ocean, the Louis' (1979) parameterization is used, It corrects for stability a neutral coefficient related to the roughness length.

The evaporative flux is calculated separately for each of the eight coexisting land surface types (see (j)). The total evaporative flux from land is then computed as an area-weighted average of the individual fluxes. The total flux includes sublimation from snow, evaporation from bare soil and from moisture intercepted by the canopy of each vegetation class, and transpiration from the dry foliage of each class. Sublimation and evaporation from intercepted canopy moisture occur at the potential rate, while canopy transpiration and evaporation from bare soil depend on the surface relative humidity which is parameterized in terms of soil moisture. Evaporation from sub-canopy soil is neglected. The surface moisture flux is computed by a bulk method and depends on the moisture gradient between the surface and the overlying air and on resistances of different kinds (aerodynamic, soil, architectural, and canopy) that vary according to surface type and/or the nature of the moisture flux (sublimation, evaporation, transpiration). See Ducoudre et al. (1993) for further details. See also (l).

The surface temperature is estimated from a surface energy balance that takes into account the radiative budget, the sensible heat flux, the latent heat flux, the heat of fusion of snow or ice, and, for sea ice only, the heat conduction from below the surface. The bulk heat is that of a layer whose depth is chosen so that vertical variation of temperature at a one-day time scale is significant.

Soil hydrology is simulated using the land-surface scheme SECHIBA of Ducoudre et al. (1993). The total depth of the soil column (corresponding to the vegetation root zone) is a constant 1.0 m. Soil moisture is computed in two layers, the upper layer being the most reactive: when precipitation exceeds evapotranspiration, the upper layer fills first; when the reverse is true, it empties first. Runoff occurs whenever the soil column is completely saturated (water depth 0.15 m). The remaining prescribed parameters for bare soil are a constant evaporative resistance and an empirical constant used to compute surface relative humidity for calculation of evaporation.

In SECHIBA, each of the seven prescribed vegetation classes interacts individually with the soil hydrology and contributes individually to the surface moisture flux. All the vegetation is assumed to have a single-story canopy that transpires or intercepts precipitation, but the canopy moisture capacity varies with the leaf area index, which is prescribed differently for each vegetation class. Different architectural and canopy resistances for evaporation/transpiration are alos prescribed for each vegetation class. See Ducoudre et al. (1993) for further details. See also (j) and (k).

The runoff water is instantaneously transferred to the ocean. 26 drainage basins are considered.

IPSL: The oceanic model OPA

The OPA OGCM has been developed at the Laboratoire d'Oceanographie Dynamique et de Climatolgie (LODYC) to study the large-scale ocean circulation and its interactions with the atmosphere and the sea ice. The general philosophy consists in solving the primitive equations with appropriate parameterizations and algorithms. This model was first used to study the tropical Atlantic. In 1988, it was rewritten, with some improvements on the physics and on the numerics (vectorization, multitasking). Version 7, which is currently used for studying the coupled atmosphere-ocean system, has been developed in 1992, with many improvements on physical parametrizations and facilities (OPA version 7, Ocean General Circulation Model, Reference Manual, September 1993). The code was adapted to the global ocean by Madec and Imbard (1996). A new version, OPA version 8, has been recently written, which is more adapted for the coupling of the ocean (OPA version 8.0, Ocean General Circulation Model, Reference Manual, Technical Report, April 1997). Presently, the version used onto a global domain is OPA7 (OPAG2 configuration and OPAICE, which includes a sea-ice model).

The use of vectorial operators relying on tensorial formalism ensures second-order accuracy on any curvilinear orthogonal grid (Marti et al., 1992). Equations are discretized on a staggered Arakawa C-grid. The discretization of vorticity in the momentum equation ensures the potential enstrophy conservation for any horizontal non-divergent flow. Time stepping for the advection, Coriolis, and pressure terms is achieved by a basic leap-frog scheme associated with an Asselin filter applied at each time step in order to avoid time splitting.

The horizontal resolution is 182 points by 152 points. A distinctive feature of the model grid is that, in the Southern Hemisphere, it is regular, while in the Northern Hemisphere, it is stretched with the pole centered on Asia to overcome the North Pole singularity (Madec and Imbard, 1996). So it does not have a geographical configuration. Typically, the horizontal resolution gives a mean size for the mesh of 2 degrees for longitude and 1.5 degree for latitude at mid-latitudes, and 0.5 degree at the equator. The bathymetry is derived from ETOPO5, a 5 minutes by 5 minutes gridded global bathymetry file provided by the Marine Geology and Geophysics Division of the National Geophysical Data Center. The model has 11 islands (Antartic, America, New Zealand, Australia, Madagascar, New Guinea, Borneo, Philippine, Cuba, Iceland, and Spitsbergen). The Bering Strait is open.

Surface to a maximum depth of 5000 m.

Finite-difference z coordinate.

There are 31 vertical levels, with 10 levels in the top 100 m. The vertical mesh is deduced from a mathematical function of z, not regular with depth. The ocean surface corresponds to the w level (kD1), and there is always a T level in the ground (kD31). The depths of the vertical levels are the following: 0.00, 10.00, 20.00, 30.00, 40.01, 50.02, 60.04, 70.09, 80.18, 90.35, 100.69, 111.36, 122.65, 135.16, 150.03, 169.42, 197.37, 241.13, 312.74, 429.72, 611.89, 872.87, 1211.59, 1612.98, 2057.13, 2527.22, 3011.90, 3504.46, 4001.16, 4500.02, and 5000.00 m.

The primitive equations are expressed in terms of potential temperature, salinity, u and v current components and surface pressure. These equations are the Navier-Stokes equations plus the classical hypotheses: turbulent closure hypothesis, spherical Earth approximation, thin-shell approximation, quasi-Boussinesq hypothesis, hydrostatic and incompressibility hypotheses. The non-linear equation of state (UNESCO, 1981) is used. It is very accurate but expensive in CPU time (it takes about 15% of the total CPU time). The rigid lid approximation is made at the surface for filtering the external gravity waves. With this approximation, the horizontal velocity can be splitted into barotropic and baroclinic components. The barotropic velocity is obtained by solving an elliptic equation for the barotropic streamfunction (a preconditioned conjugate gradient solver is used).

IPSL-oc: Diffusion

The lateral turbulent fluxes are assumed to depend linearly on the lateral gradients of large-scale quantities. The resulting lateral diffusive and dissipative operators are of second order. With the global ocean grid, operators acting on momentum, temperature, and salinity are of Laplacian type.

The horizontal eddy viscosity varies spatially. Outside the equatorial regions it takes a value of 40000 m² s-1. Near the equator it is reduced to 2000 m2 s-1, except near the western coasts where it is kept fixed at 40000 m2 s-1. A forward scheme is used.

For active tracers, the diffusion is along the isopycnal surfaces (more precisely, neutral surfaces). The eddy diffusivities are taken as 2000 m2 s-1.

The vertical eddy viscosity and diffusivities are computed by a 1.5 order Turbulent Kinetic Energy (TKE) closure scheme (Blanke and Delecluse, 1993) which allows the formulation of the mixed layer as well as a minimum diffusion in the thermocline. An implicit scheme is used for vertical diffusive processes.

The solar radiation penetrates the top meters of the ocean (Blanke and Delecluse, 1993). The absorption is very selective in the first meters of the ocean with a preference for the longer and the shorter wavelengths. The downward irradiance I(z) is formulated with two extinction coefficients (Paulson and Simpson, 1977), whose values correspond to a Type I water in Jerlov's classification (i.e., the most transparent water).

The longwave radiation is assumed to be absorbed exclusively in the top layer of the ocean.

The static instabilities are resolved in the turbulent closure scheme.

See (a).

Zero fluxes of heat and salt and no-slip conditions are applied at solid boundaries.

+ OPAG2 configuration

This version does not include a sea-ice component but a simple pseudo-parametrization that only involves a test on the sea-surface temperature. Ice appears when the ocean surface temperature becomes less than the freezing temperature. A constant heat flux is applied to the ocean under the ice (-2 W m-2 in the Arctic Ocean and -4 W m-2 in the Southern Ocean), and the AGCM computes its own surface heat flux.

+ OPAICE configuration

A realistic representation of the ice distribution, including open-water and several ice types, is performed. The different ice types behave in different ways: for example, thin ice grows and decays more rapidly that thick ice. The existence of a variable fraction of leads in the ice pack is also simulated. The model takes into account the presence of snow on top of sea ice, the storage of latent heat, and the transformation of snow into ice due to snow flooding. The snow albedo evolves with the snow age. The model structure is initially based on Hibler (1980). Particularly, it uses a distribution function associated with the different ice types. Only the thermodynamic aspects are included in the version used for coupled studies (the dynamical part is under development). The most important difference with Hibler's model is an improvement of the ice distribution description. In Hibler (1980), this function only depends on ice thickness; in our model, it depends on two variables : the thickness and the interface temperature. These variables define axes in a phase space in which the evolution of the distribution function is written. We use a multi-level model for the prediction of temperature. The number of levels can be chosen arbitrarily. In experiments with the coupled model, a one-level version is used.

In absence of ice, a flux boundary condition is used in the temperature vertical diffusion equation. In presence of ice, the sea ice is assumed to be totally included in the first level of the ocean. The ocean surface temperature is fixed to the freezing point. The ocean supplies heat flux to the ice. This flux is used by the ice to determine its growth rate. This rate is used to determine the water flux at the top of the ocean.