References:

The model has been documented in Johns (1996) and in Johns et al (1997)

The AGCM ARPEGE-Climat is an adaptation of the ARPEGE model (an.acronym for "Action de Recherche Petite Echelle Grande Echelle", i.e."Research Project on Small and Large Scales) for climate studies. The ARPEGE model was developed at Meteo-France in collaboration with ECMWF where the model is called IFS (Integrated Forecast System) Physical parameterizations of ARPEGE and IFS are generally different. The version detailed below is the version 2 of the ARPEGE-Climat model, released to the scientific community at the beginning of 1996 in the context of a community climate model project. The OGCM was developed by the Laboratoire d'Oceanographie Dynamique et de Climatologie (LODYC, Paris). The two models are coupled through the software OASIS developed by the Centre Europeen de Recherche et de Formation Avancee en Calcul Scientifique (CERFACS, Toulouse).

A complete description of the physics and dynamics of the model is included in a technical report entitled "Community modelling, 1996: ARPEGE-Climat, Version 2, I. Algorithmic documentation". In addition to this documentation, three complementary reports describe how to use the model (II. User Guide), the main characteristics of the code (III. Code Documentation) and the climatology of the model (IV. Climate Validation). All these reports are available under request from Michel Deque (CNRM, Toulouse, E-mail: Michel.Deque@meteo.fr). Several papers cover the description of the different components of the model physics. Some key publications on the model design are Deque et al (1994), Deque and Piedelievre (1995), Geleyn et al (1995).

The horizontal representation is spectral (spherical harmonic basis functions) with transformation to a Gaussian grid for calculation of nonlinear quantities and some physics. The Gaussian grid is reduced longitudinally near the poles so that its horizontal resolution is everywhere approximately the same. Through a conform grid transform, the effective resolution of the model may be varied, locating the centre of highest resolution at any geographical point (Courtier and Geleyn, 1988).

A classical resolution for coupled experiments is the spectral triangular 31 (T31) with a reduced Gaussian grid implying a quasi uniform horizontal resolution of 400 Km (48 points in latitude, 96 points in longitude).

Finite-difference with hybrid sigma-pressure coordinate.

With the 19 level version, for a surface pressure of 1000 hPa, 5 levels are below 800 hPa and 6 levels are above 200 hPa.

The model uses a spherical representation of the primitive equations expressed in terms of vorticity and divergence, temperature and specific humidity, and the natural logarithm of surface pressure (or surface pressure itself as an optional choice) (Simmons and Burridge, 1981). The possibility exists of including arbitrary scalar variables which obey equations similar to the specific humidity equation. In the present model, the ozone mixing ratio is such a variable. The mass of non-precipitating liquid and solid water and of ozone are neglected. The mass flux of precipitation/evaporation may be taken into account in the atmosphere's total mass budget as an optional choice. Also in option is the treatment of the advection through a semi-lagrangian scheme.

The atmosphere is assumed to be a mixture of two perfect gases, water vapor and dry air. This implies a combination of two gases constants and two specific heat. The latent heat of vaporization (sublimation) and the saturated partial pressure of liquid water (ice) are computed by the Clausius Clapeyron relations, implying a dependence of the latent heat of vaporization (sublimation) with the temperature.

A Linear del⁶ horizontal diffusion is applied on constant hybrid sigma-pressure surfaces to vorticity, divergence, temperature, specific humidity and ozone mixing ratio. This horizontal diffusion is calculated in the spectral space. The diffusivity is constant below 100 Hpa. Above 100 Hpa, the diffusivity increases as the inverse of pressure, yielding very strong diffusion in the model stratosphere.

A mesospheric drag consisting in a relaxation to zero for the wind and to a standard value for temperature, is added for levels above 1 Hpa.

The vertical exchange of static energy, momentum and humidity, depends on exchange coefficients estimated as functions of the vertical stability, surface roughness, wind shear and a mixing length (Louis, 1979; Louis et al, 1982). The mixing length is a function of the asymptotic lengths and of a typical boundary layer depth. This scheme is applied to the whole atmospheric depth, but the exchange coefficients vanish well above the boundary layer.

The gravity-wave drag scheme accounts for the effect of the subgrid-scale orography. The surface momentum flux is proportional to the Brunt-Vaissala frequency, to an effective surface wind and to the root mean square of the unresolved orography in the direction of this wind. The effective surface wind accounts for the depth of the boundary layer generated by the unresolved orography. The angle between direction of the surface momentum flux and the effective wind depends on the anisotropy of this orography. The magnitude of the momentum flux is always maximum at the surface and its direction is invariant. Three effects are taken into account : the dissipation, the resonance and the reflexion. Following the instability criterion of Lindzen , the dissipative flux is proportional to the square of the Froude number and the momentum flux thus vanishes above the first level where the wind becomes perpendicular to the surface wind direction. The resonant term follows the experimental results of Clark and Peltier (1984). The surface momentum flux may increase and the momentum flux profile is then taken linear from the surface to a calculated critical level. In case of attenuation, the momentum flux is truncated to the reduced surface value. The reflection occurs at the first level where the Brunt-Vaissala frequency vanishes. The resulting contribution to the momentum flux is linear from the surface to this level in order to keep unchanged the surface value and to repeal the flux at the reflection level.

Following the recent work of Lott and Miller (1995), an additional drag is introduced in the model layers corresponding to the mountain boundary layer, to better parameterize the blocking of the flow by the unresolved orography. As the first component of the drag (see above), the corresponding surface momentum flux accounts for the anisotropy of the unresolved orography.

The parameterization of convective wave drag assumes that the source of perturbation is proportional to the convective precipitation rate. The waves are assumed to move upwards and in the direction of the wind at the top of the cloud. The horizontal phase velocity of the waves is constant in the vertical and corresponds to the projection of the wind at the cloud bottom onto the direction of propagation. As for the orography waves, above the cloud top, the Lindzen criterion is used to determine the wave dissipation. Inside the cloud, to account for an exchange of momentum between the cloud layers, the momentum flux profile is taken linear from zero at the cloud bottom to its calculated value at the cloud top.

The CO² concentration is supposed constant whatever the location and the altitude.

A specific parameterization calculates ozone fluxes due to photochemical sources and sinks. It consists in a linear scheme (Cariolle and Deque, 1986), adjusted from a 2D statospheric model including complex photochemical reactions.

The profile of aerosol is taken identical at each location. The optical depth, single scattering albedo and assymmetry factor, result from a combination of the optical properties of five aerosol types proposed by Tanre et al (1983).

Two different schemes may be used. The FMR scheme (Morcrette, 1990), is similar to the one described in the UCL's coupled model documentation. The second scheme is derived from the methods described in Geleyn and Hollingsworth (1979) with optical properties characterized by a few coefficients, following Ritter and Geleyn (1992). It is simplified, so that it can be called at each time step. We give some more characteristics of this scheme in the following.

Each vertical layer is supposed to be horizontally homogeneous and mean values of cloudiness and temperature are calculated for each layer. There are only two spectral intervals, one in the solar visible range and one in the infrared one. Following the delta two-stream calculation, which supposes that the diffusion is isotropic, only two diffusive fluxes, one upwards and the other one downwards, are considered. The optical depth for the solar radiation is taken exactly along the parallel downward beam, and estimated for the diffuse beams. For thermal radiation, ``cooling to space'' calculations are handled exactly, while exchange with the surface and between layers are intentionally underestimated through an approximation, so that the scheme remains stable.

Cloud radiative properties are specified after Stephens (1979) for water clouds with eight different droplet size distributions that are related to diagnostic cloud liquid water content following Betts and Harshvardhan (1987). Optical properties of ice clouds are also included. A random overlap of adjacent cloud layers is assumed, which consists in calculating the radiative fluxes as a linear combination (which depends on the cloudiness) of the fluxes obtained at each level for clear sky and cloudy sky conditions. A maximum overlap version is also coded in the model but is more expensive, since it involves solving the complete system of flux calculation once for clear sky and once for cloudy sky conditions. It is better suited for high vertical resolution versions.

The deep convection scheme used has been described by Bougeault (1985). The deep convection occurs under two conditions: a convergence of humidity at low layers is required and the vertical temperature profile must be unstable. Within the grid mesh, the total area of convective clouds is neglected and the tendencies of grid-scale temperature and specific humidity are thus taken equal to those of the cloud environment variables. The scheme uses the mass-flux concept where the vertical ascent in the cloud is compensated by a grid-scale subsidence which warms and dries the environment. To account for the detrainment of cloud air and of the subgrid-scale vertical transport, the grid-scale temperature and specific humidity are relaxed towards the profiles of a specific cloud type. These profiles follow a moist adiabatic but account for the entrainment of environmental air. The mass flux at each level is proportional to the square root of the difference of moist static energy between the cloud and the grid-scale. The constant of proportionality is calculated by means of a Kuo- type closure hypothesis. The partition of the available moisture between precipitation and moistening of the environment is determined from the conservation of the integrated moist static energy. The evaporation below the cloud is parameterized. A wind tendency due to convection is calculated as the sum of a subsidence term and a mixing term representing the subgrid-scale vertical transport.

The shallow convection is taken into account through a modification of the Richardson number (Geleyn, 1986).

Cloud cover is diagnostically determined. The stratiform cloudiness are calculated as a function of the humidity and a critical humidity profiles. The total convective cloudiness is a function of the convective precipitation. It is equally distributed amongst the layers where convection occurs. The cloud cover at each level are calculated by combining the stratiform and convective fractions with random overlap or with maximum overlap according to a switch value.

The condensated water is diagnostically determined as a function proportional to the square of the pressure, and to the derivative of saturated specific humidity with respect to pressure along the moist adiabatic. This condensated water is partitioned into solid and liquid water as a function of temperature.

The precipitation scheme assumes that the atmosphere is adiabatic and not allowed to become oversaturated. All the condensed water is supposed to precipitate in one time-step and precipitation is at the same temperature as the environment. An estimate of the proportion of falling ice phase is evaluated as a function of the temperature. The rainfall rate is computed from a Kessler-type formula. Further details on this parameterization may be found in Matveev (1982), Rockel et al (1991) and Clough and Francks (1991).

The PBL is generally represented by he first 6 levels of the model, which, for a surface pressure of 1000 Hpa, correspond to pressures of 995, 981, 954, 912, 858, 793 Hpa (see also (a), (j) and (k)).

The snow parameterization, combined with the ISBA soil- vegetation scheme (l), is described in Douville et al (1995).

The proportion between liquid and snow convective precipitation is determined as a function of the temperature and the pressure of the level closest to the surface. The snow stratiform precipitation are calculated as a function of the temperatures of all the precipitating levels. The total snow precipitation feeds a reservoir of snow content. Melting of snow is calculated when the snow temperature, a combination of the surface and the first soil layer temperature (see (l)) according to the fraction of vegetation covered by snow, is greater than the melting point temperature. Sublimation is calculated as part of the evaporation flux and snow melt is added to the soil moisture content (see (l)).

The fraction of bare soil covered by snow depends on snow amount and on the roughness length of the orography. The fraction of vegetation covered by snow is estimated as a function of the snow depth and the height of the vegetation, following the formula proposed by Pitman et al (1991). A value of the snow density is computed at each time step from the density of new snow and the one of old snow determined through an exponential increase with time (Verseghy, 1991) The snow amount, snow cover and snow density affect the surface characteristics (see (j)).

The ISBA land surface scheme (see (l)) requires the definition of parameters describing the land/cover at each continental grid point. The original datasets consist of 1deg.X1deg. fields of vegetation (primary and secondary) and soil types (Wilson and Henderson Sellers, 1987), and the percentage of clay and sand in the soil (Webb et al, 1991). The inferred soil/vegetation parameters are aggregated at the model grid scale following a procedure similar to the one described in Manzi and Planton (1994).

The roughness length linked to the subgrid orography is calculated by means of an orography data set on a 50min.X50min. grid (from NASA). The roughness length linked to the vegetation cover, results from the aggregation procedure and incorporates a seasonal cycle for deciduous and cultivated vegetation. The resultant roughness length combines the two previous and is also modulated by the snow cover. A thermal roughness length is also computed.

The surface albedo are function of the solar zenith angle. At each grid point, the albedo value for diffuse radiation is a surface weighted combination of bare soil, vegetation and snow albedo. The bare soil albedo is inferred from satellite observations and the vegetation data. The vegetation albedo is calculated by the aggregation procedure. During the melting period, the snow albedo is calculated through a prognostic equation following Verseghy (1991), while a weak linear decrease is imposed otherwise.

The emissivity varies geographically (CLIMAP dataset, 1981) and is combined with a snow value according to the snow cover fraction.

The vegetation parameters (fraction of vegetation, Leaf Area Index (LAI), minimum stomatal resistance) are aggregated at the model grid scale. For the LAI a seasonal cycle for deciduous and cultivated vegetation, is also included.

The depth of the active soil layer required by the land surface scheme, is assigned as the larger of the vegetation root depth and a soil depth determined from the soil type.

The thermo-hydric properties of the soil are function of the soil texture and the soil water content as described in Noilhan and Planton (1989). The soil thermal inertia also depends on the vegetation cover, and on the snow cover and density.

Net solar flux is computed from the incoming direct and diffuse solar fluxes and the corresponding surface albedos. The net infrared radiation results from the incoming infrared radiation, the emitted infrared flux from the Planck equation with a geographically varying emissivity (see (d) and ((j))

At the surface, the vertical diffusive fluxes are calculated through drag coefficients proposed by Louis et al. (1979). They are functions of the Karman constant, the roughness length and the Richardson number. The value of the roughness length as well as the similarity functions are different for the velocity and for the thermodynamical variables.

The total evaporation flux is the sum of three components: evaporation from bare ground, foliage transpiration and interception loss. The evaporation from bare ground is function of the surface humidity determined in terms of surface soil water content. The transpiration is modulated by a resistance function of the LAI, the minimum stomatal resistance and several limiting factors (Photosynthetically Active Radiation, soil water, vapor pressure deficit, temperature). The sublimation and the evaporation of intercepted water on the foliage are at the potential rate. The surface of foliage covered by intercepted water is determined as a function of the content of the corresponding reservoir, following Deardorff (1978).

The surface precipitation flux (see (g)) fills in part the reservoir of rain intercepted by the foliage, which collecting surface depends on the fraction of vegetation and on the amount of convective rain. The remaining part reaches the soil surface.

The ISBA (Interface Sol Biosphere Atmosphere) land-surface scheme, developed by Noilhan and Planton (1989) and improved by Mahfouf et al (1995).

The temperature is calculated at four levels in the soil including the surface, and without climatological restore term of the deeper layer. The basic equations extend the original formalism of Bhumralkar (1975) and Blackadar (1976).

The soil hydrology is characterized by the surface volumetric water content, the mean water content of the active soil layer and the reservoir of rain intercepted by the canopy. The evolution of the two first reservoirs is computed with a scheme based on the force-restore method for water transfer (Deardorff, 1978). Run-off occurs when either the surface or the mean reservoir reach the maximum volumetric water content depending on the soil texture. The Run-off of the rain intercepted by the canopy is calculated in a similar way, and is added to the soil surface precipitation. The maximum value of the corresponding reservoir is proportional to the LAI and the fraction of vegetation. A deep drainage by gravity (Goutorbe et al, 1989) is included by a restore term towards the field capacity when the mean water content exceeds it.

Detailed information about the oceanic model can be found in the description of this part of the IPSL coupled-model.

Overall documentation of the LMD5 model is provided by Polcher et al. (1991).

Other key model documents include publications by Sadourny and Laval (1984), Laval et al. (1981), Le Treut and Li (1991), Ducoudré et al. (1993), Le Treut et al. (1994), and Grenier (1997). Computational aspects are described by Butel (1991).

Horizontal representation

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.

There are 50 grid points equally spaced in the sine of latitude and 64 points equally spaced in longitude. (The mesh size is 225 km north-south and 625 km east-west at the equator, and is about 400 x 400 km at 50 degrees latitude.)

Surface to about 4 hPa. For a surface pressure of 1000 hPa, the lowest level is at 979 hPa.

Finite-difference sigma coordinate.

There are 11 unevenly spaced sigma levels. For a surface pressure of 1000 hPa, 3 levels are below 800 hPa and 2 levels are above 200 hPa.

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.

UCL-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. The radiative effects of aerosols are not included.

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^-²

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.

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.15 K, and we assume it varies linearly between these values at intermediate temperatures.

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 (= 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' = 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.

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. In order to account for the large difference in the turbulent heat fluxes over sea ice and leads or polynyas at high latitudes and to ensure conservation of heat and freshwater while transferring fluxes between the AGCM and the OGCM, each surface grid box is divided into appropriate fractions of land, ocean, and eventually sea ice, and the vertical diffusion processes occurring in the planetary boundary layer are resolved separately over the continental, oceanic, and sea-ice portions (Grenier, 1995). Mixing between the subgrid columns of air is accomplished by lateral diffusion, with the diffusivity increasing with altitude. See also (a), (j), (k), and (l).

If the air temperature at the first level above the surface (at sigma = 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, two different formulations are available. Following Bunker (1976), transfer coefficients are calculated using an empirical formula based on surface wind speed and on the temperature difference between the ocean and the surface air. The other possibility is to use Louis' (1979) parameterization, which 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 Ducoudré 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 Ducoudré 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 Ducoudré 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.

Detailed information about the model physics and numerics can be found in Deleersnijder and Campin (1995), Goosse et al. (1997a,b), and Campin (1997). Fichefet and Morales Maqueda (1997) give an exhaustive description of the sea-ice component.

Horizontal representation

Finite differences on a staggered B-grid (e.g., Arakawa and Lamb, 1977).

The horizontal resolution is of 3 degrees by 3 degrees. To cope with the singularity at the North Pole associated with the geographical spherical coordinates, two spherical grids connected in the equatorial Atlantic are used: a grid with its poles located on the geographical equator for the North Atlantic and the Arctic, and a classical latitude-longitude grid for the rest of the World Ocean (Deleersnijder et al., 1993).

Surface to a maximum depth of 5500 m.

Finite-difference z coordinate.

There are 20 levels. The bases of the corresponding layers lie at the following depths: 10, 22, 36, 54, 76, 104, 139, 187, 253, 346, 484, 694, 1007, 1443, 1992, 2622, 3304, 4018, 4752, and 5500 m.

The UCL OGCM is a primitive-equation, free-surface model that rests on the usual set of approximations, i.e., the hydrostatic equilibrium and the Boussinesq approximation. Its dependent variables are the three components of the velocity, the sea-surface elevation, the potential temperature, and the salinity. The equation of state is that of Eckart (1958). See Deleersnijder and Campin 1995 for further details.

The horizontal eddy viscosity and diffusivity are taken as 10⁵ m² s^-1 and 150 m² s^-1, respectively. No isopycnal mixing scheme is used.

The Mellor and Yamada's (1974, 1982) level 2.5 turbulence closure scheme as modified by Kantha and Clayson (1994) is employed to compute the vertical eddy viscosity and diffusivity. The stability functions are derived from the quasi-equilibrium model of Galperin et al. (1998). The turbulent kinetic energy is computed from a prognostic equation. The mixing length is taken as a function of vertical stability and turbulent kinetic energy in stratified regions (e.g., Kantha and Clayson, 1994) and as a function of depth in unstratified regions. The vertical eddy viscosity is not allowed to drop below 10^-4 m² s^-1 throughout the whole water column. A minimum value is also prescribed for the vertical eddy diffusivity. This value is of 10^-5 m² s^-1 in the first 400 m and then increases with depth to 1.1 x 10^-4 m² s^-1 at the bottom in accordance with Bryan and Lewis (1979).

The parameterization of Paulson and Simpson (1977) is used to compute the differential absorption of shortwave radiation in the ocean interior. Regarding the ocean turbidity, the geographical distribution of Jerlov's (1976) water types proposed by Simonot and Le Treut (1986) is utilized.

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

When static instability occurs, the vertical eddy diffusivity is increased to 100 m² s^-1 to simulate convection.

See above (diffusion).

The bottom stress is determined from a quadratic drag law using a drag coefficient of 1.5 x 10^-3.

The parameterization developed by Campin (1997) for the export of dense waters formed over continental shelves is used. In this formulation, the amount of water exported along the continental slope is a function of the density difference between the shelf and the open ocean. See Campin (1997) for further details.

The non-overlapping of the two grids in the area of the Bering Strait prevents an explicit calculation of the water flow through the strait. So, this transport needs to be parameterized. In accordance with the geostrophic control theory, the throughflow is taken to be proportional to the cross-strait sea-level difference. See Goosse et al. (1997a) for further details.

The sea-ice model (Fichefet and Morales Maqueda, 1997) incorporates parameterizations of the most relevant thermodynamic and dynamic sea-ice processes.

The model component that determines the vertical growth and decay of the ice due to thermodynamic processes is essentially an improved version of Semtner's (1976) three-layer model. Within the ice-covered portion of each grid cell, sea ice is supposed to be a horizontally homogeneous slab of ice (divided into two layers of equal thickness) on which snow may accumulate when the surface temperature is below the melting point. Internal temperatures of snow and ice are governed by a one-dimensional heat diffusion equation. Effective snow and ice thermal conductivities are used to account for the fact that the unresolved ice floes of varying thickness contribute differently to the average heat conduction. These conductivities are determined by assuming that the snow and ice thicknesses are uniformly distributed between zero and twice their mean value over the ice-covered part of the grid cell. The mechanism of "brine damping'' is represented following the approach of Semtner (1976), wherein the solar energy absorbed inside the ice is stored in a heat reservoir that represents internal meltwater. Energy from this reservoir is utilized to keep the temperature of the upper ice layer from dropping below the freezing point in autumn, thereby simulating the release of latent heat through refreezing of the internal brine pockets. When the load of snow is large enough to depress the lower boundary of the snow layer under the water level, seawater is assumed to infiltrate the entirety of the submerged snow and to freeze there, forming a snow-ice cap.

To take into consideration the existence of leads and polynyas within the ice pack, we introduce the concentration variable A, which is defined as the fraction of the grid cell area covered by ice. Variations of A depend on the heat budget of the open water area. See Fichefet and Morales Maqueda (1996) for further details about the lateral growth and decay of the ice.

For the momentum balance, the ice is considered as a two-dimensional continuum in dynamical interaction with the atmosphere and ocean. The advection of momentum is neglected. Following Hibler (1989), sea ice is assumed to have a viscous-plastic constitutive law.

The physical fields that are advected are the ice concentration, the snow volume per unit area, the ice volume per unit area, the snow enthalpy per unit area, the ice enthalpy per unit area, and the latent heat contained in the brine reservoir per unit area.

The sea-ice model supplies heat flux, salt flux, and momentum-exchange boundary conditions for the top of the ocean. The ocean model, in turn, supplies current and heat-exchange information to the ice model. It should be noted that the freezing point of seawater is taken as a function of salinity.