PMIP Documentation for UKMO

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

PMIP documentation for UKMO

United Kingdom Meteorological Office: Model UKMO HADAM2 (2.5x3.75 L19) 1997

PMIP Representative(s)

Chris Hewitt, Hadley Centre for Climate Prediction and Research, United Kingdom Meteorological Office,London Road,London Road, Bracknell, Berkshire RG12 2SY, United Kingdom; Phone: +44 1344 85 4520; Fax:+44 1344 85 4898; e-mail: cdhewitt@meto.gov.uk ;

Model Designation

for fixed experiments:

UKMO HADAM2 (2.5x3.75 L19) 1997 or UKMO - Unified Model vn3.2 AGCM - 2x4L19

for computed experiments:

UKMO HADAM2b (2.5x3.75 L19) 1997 or UKMO - Unified Model vn3.4 Slab GCM - 2x4L19

Model Identification for PMIP

UKMO

PMIP run(s)

0fix, 6fix, 0cal, 21cal

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

Model Lineage

The UKMO HADAM2 model part of the line of Unified Models (UM) intended to provide a common framework for forecasting and climate applications (cf. Cullen 1993 ). The dynamical formulations are those described by Bell and Dickinson (1987) ; the physical parameterizations are substantially modified from those of an earlier UKMO model documented by Slingo (1985) .

The UKMO HADAM2b model is a coupled atmosphere-mixed layer ocean model.

Model Documentation

The atmospheric model is described in :

Hewitt, C.D., and Mitchell, J.F.B., 1996: GCM simulations of the climate of 6 kyr BP : Mean changes and interdecadal variability. J. Climate, 9, 3505-3529.

Johns,T.C., R.E. Carnell, J.F. Crossley, J.M. Gregory, J.F.B. Mitchell, C.A. Senior, S.F.B. Tett, and R.A. Wood, 1997. The second Hadley Centre coupled ocean-atmosphere GCM : Model Description, Spinup and Validation. Climate Dyn., 13, 103-134.

Computed experiments are described briefly in:

Hewitt, C.D., and Mitchell, J.F.B., 1997 : Radiative forcing and response of a GCM to ice age boundary conditions : cloud feedback and climate sensitivity. Climate Dynamics, 13, 821-834.

There is a Climate Research Technical Note by Tim Johns (1996), number 71 which describes the UKMO HADAM2 model in detail.

The differences between HADAM2b and HADAM2 are briefly described in the Hewitt and Mitchell (1997) paper. The PMIP documentation (below) refers to HADAM2.

Numerical/Computational Properties

Horizontal Representation

Fourth-order finite differences on a B-grid (cf. Arakawa and Lamb 1977 , Bell and Dickinson 1987 ) in spherical polar coordinates. Mass-weighted linear quantities are conserved, and second moments of advected quantities are conserved under nondivergent flow.

Horizontal Resolution

2.5 x 3.75-degree latitude-longitude grid.

dim_longitude*dim_latitude: 96*73

Vertical Domain

Surface to about 5 hPa; for a surface pressure of 1000 hPa, the lowest atmospheric level is at about 997 hPa.

Vertical Representation

Finite differences in hybrid sigma-pressure coordinates after Simmons and Strüfing (1981) . Mass and mass-weighted potential temperature and moisture are conserved. See also Horizontal Representation.

Vertical Resolution

The model uses a hydrid cor-ordinate in the vertical (the top 3 levels are pressure layers, and the bottom 4 are terrain-following sigma layers, with a smooth transition in between) with the following levels: 0.005, 0.015, 0.030, 0.057, 0.100, 0.149, 0.200, 0.250, 0.300, 0.355, 0.422, 0.505, 0.599, 0.699, 0.793, 0.870, 0.930, 0.975, 0.997.

There are 19 unevenly spaced hybrid levels. For a surface pressure of 1000 hPa, 4 levels are below 800 hPa and 7 levels are above 200 hPa.

Computer/Operating System

The PMIP simulation was run on a Cray C90 computer in a UNICOS environment.

Computational Performance

For the PMIP experiment, about 2 minutes Cray C90 computation time per simulated day (about half this time being associated with output postprocessing).

Initialization

For the PMIP experiment, the model atmosphere, soil moisture, and snow cover/depth were initialized for 1 June from a previous model simulation. Snow mass for areas of permanent land ice was initially set to 5 x 10⁴ kg/(m²). The model was then integrated forward for months before starting the fix experiments.

Time Integration Scheme(s)

Time integration proceeds mainly by a split-explicit scheme, where the solution procedure is split into "adjustment" and "advection" phases. In the adjustment phase, a forward-backward scheme that is second-order accurate in space and time is applied. The pressure, temperature, and wind fields are updated using the pressure gradient, the main part of the Coriolis terms, and the vertical advection of potential temperature. In the advective phase, a two-step Heun scheme is applied. A time step of 30 minutes (including a 10-minute adjustment step) is used for integration of dynamics and physics, except for full calculation of shortwave/longwave radiation once every 3 hours. In addition, an implicit scheme is used to compute turbulent vertical fluxes of momentum, heat, and moisture in the planetary boundary layer (PBL). Cf. Cullen et al. (1991) for further details. See also Diffusion, Planetary Boundary Layer, and Surface Fluxes.

Smoothing/Filling

To prevent numerical instability, the orography is smoothed in high latitudes (see Orography), and Fourier filtering is applied to mass-weighted velocity and to increments of potential temperature and total moisture. Negative values of atmospheric moisture are removed by summing the mass-weighted positive values in each horizontal layer, and rescaling them to ensure global moisture conservation after the negative values are reset to zero.

Sampling Frequency

For the PMIP simulation, the model history is written once every 6 hours. (All average quantities in the PMIP monthly-mean standard output data are computed from samples taken at every 30-minute time step.)

Dynamical/Physical Properties

Atmospheric Dynamics

Primitive-equation dynamics, formulated to ensure approximate energy conservation, are expressed in terms of u and v winds, liquid/ice water potential temperature, total water, and surface pressure (cf. White and Bromley 1988 ).

Diffusion

Linear conservative horizontal diffusion is applied at fourth-order (Ñ⁴) to moisture, and at sixth-order (Ñ⁶) to winds and to liquid water potential temperature (cf. Cullen et al. 1991).

Stability-dependent, second-order vertical diffusion of momentum and of conserved cloud thermodynamic and water content variables (to include the effects of cloud-water phase changes on turbulent mixing), operates only in the PBL. The diffusion coefficients are functions of the vertical wind shear (following mixing-length theory), as well as surface roughness length and a bulk Richardson number that includes buoyancy parameters for the cloud-conserved quantities (cf. Smith 1990b ). See also Cloud Formation, Planetary Boundary Layer, and Surface Fluxes.

Gravity-wave Drag

The parameterization of orographic gravity-wave drag follows a modified Palmer et al. (1986) scheme, as described by Wilson and Swinbank (1989) . The drag is given by the vertical divergence of the wave stress. Near the surface, the stress is equal to the product of a representative mountain wave number, the square of the wave amplitude (taken to be the subgrid-scale orographic variance--see Orography), and the density, wind, and Brunt-Vaisalla frequency evaluated in near-surface layers. At higher levels, the stress is given by this surface value weighted by the projection of the local wind on the surface wind. If this projection goes to zero, the stress is also zero; otherwise, if the minimum Richardson number falls below 0.25, the gravity wave is assumed to break. Above this critical level, the wave is maintained at marginal stability, and a corresponding saturation amplitude is used to compute the stress.

Solar Constant/Cycles

For 0fix, 6fix, 0cal and 21cal, the solar constant is the AMIP-prescribed value of 1365 W/(m²). The orbital parameters and seasonal insolation distribution are calculated after PMIP recommendations. Both seasonal and diurnal cycles in solar forcing are simulated.The seasonal cycle of solar insolation is based on a 360-day year (each month 30 days in length), with the date of perihelion adjusted to minimize discrepancies (cf. Ingram 1993).

Chemistry

The carbon dioxide concentration are prescribed value of 323, 262 and 229 ppm for control, 6ka and 21ka, respectively. Zonally averaged monthly ozone profiles are specified from the climatology of Keating et al. (1987) above hybrid level 0.0225 (see Vertical Representation), and from satellite data of McPeters et al. (1984) below this level. Radiative effects of water vapor and clouds, but not those of aerosol, are also included (see Radiation).

Radiation

The radiation schemes are as described by Ingram (1993)). Incoming insolation is based on a 360-day year (see Solar Constant/Cycles). Shortwave computations follow Slingo (1985)) , extended to include 4 spectral intervals (with boundaries at 0.25, 0.69, 1.19, 2.38, and 4.0 microns) and the interactive cloud optical properties of Slingo (1989) . Rayleigh scattering is represented by reflection of 3 percent of the incoming insolation before any interaction with the atmosphere. Shortwave absorption (by ozone, carbon dioxide, and water vapor) is treated by use of look-up tables. The effects of pressure broadening on carbon dioxide and water vapor absorption bands are included; the overlap of these bands is assumed to be random within each spectral interval.

Calculation of longwave fluxes follows the method of Slingo and Wilderspin (1986) . Absorption by water vapor, carbon dioxide, and ozone is treated in 6 spectral bands (with boundaries at 0.0, 4.0x10⁴, 5.6x10⁴, 8.0x10⁴, 9.0x10⁴, 1.1x10⁵, and 1.2x10⁵ m^-1), using exponentials for the water vapor continuum and look-up tables otherwise. The temperature dependence of e-type absorption in the water vapor continuum is after Roberts et al. (1976) . Pathlengths of gaseous absorbers are scaled to account for pressure-broadening effects(and for carbon dioxide, also for temperature-broadening effects).

The convective cloud in each grid box is confined to a single tower (with full vertical overlap). All absorption of shortwave radiation by the convective cloud occurs in its top layer, with the remainder of the beam passing unimpeded through lower layers. In each vertical column, the shortwave radiation interacts only with 3 stratiform clouds defined by vertical domain: low (levels 1 to 5), middle (levels 6 to 9), and high (levels 10 to 18), all randomly overlapped. The cloud with the greatest area in each domain interacts with the shortwave radiation, and it is assigned the cloud water content of all the layer clouds in the domain. (Neither the shortwave nor longwave scheme allows cloud in the top layer, however.) In the longwave, similar to the formulation of Geleyn and Hollingsworth (1979) , clouds in different layers are treated as fully overlapped if there is cloud in all intervening layers, while clouds separated by cloud-free layers are treated as randomly overlapped.

Sunlight reflected from a cloud is assumed to pass directly to space without further cloud interactions, but with full gaseous absorption; the surface albedo (see Surface Characteristics) is adjusted to account for the increased absorption due to multiple reflections between clouds and the surface. Cloud shortwave optical properties (optical depth, single-scattering albedo, and asymmetry factor) are calculated from the cloud water path (CWP), the effective radius of the drop-size distribution (7 microns for water and 30 microns for ice), and the sun angle, following the Practical Improved Flux Method of Zdunkowski et al. (1980) . Cloud longwave emissivity is a negative exponential function of CWP, with absorption coefficient 65 m²/kg for ice clouds and 130 m²/kg for water clouds. See also Cloud Formation and Precipitation.

Convection

Moist and dry convection are both simulated by the mass-flux scheme of Gregory and Rowntree (1990)) that is based on the bulk cloud model of Yanai et al. (1973) . Convection is initiated if a parcel in vertical layer k has a minimum excess buoyancy beta that is retained in the next higher level k + 1 when entrainment effects and latent heating are included. The convective mass flux at cloud base is taken as proportional to the excess buoyancy; the mass flux increases in the vertical for a buoyant parcel, which entrains environmental air and detrains cloud air as it rises. Both updrafts and downdrafts are represented, the latter by an inverted entraining plume with initial mass flux related to that of the updraft, and with detrainment occurring over the lowest 100 hPa of the model atmosphere.

When the parcel is no longer buoyant after being lifted from layer m to layer m + 1, it is assumed that a portion of the convective plumes has detrained in layer m so that the parcel in layer m + 1 has minimum buoyancy b. Ascent continues until a layer n is reached at which an undiluted (without entrainment) parcel originating from the lowest convectively active layer k would have zero buoyancy, or until the convective mass flux falls below a minimum value. Cf. Gregory (1990) and Gregory and Rowntree (1990) for further details.

Cloud Formation

The convection scheme (see Convection) determines the vertical extent of subgrid-scale convective cloud, which is treated as a single tower in each grid box. The convective cloud base is taken as the lower boundary of the first model layer at which saturation occurs, and the cloud top as the upper boundary of the last buoyant layer. The fractional coverage of each vertical column by convective cloud is a logarithmic function of the mass of liquid water condensed per unit area between cloud bottom and top (cf. Gregory 1990).

Large-scale (stratiform) cloud is prognostically determined in a similar fashion to that of Smith (1990a) . Cloud amount and water content are calculated from the total moisture (vapor plus cloud water/ice) and the liquid/frozen water temperature, which are conserved during changes of state of cloud water (i.e., cloud condensation is reversible). In each grid box, these cloud-conserved quantities are assumed to vary (because of unresolved atmospheric fluctuations) according to a top-hat statistical distribution, with specified standard deviation. The mean local cloud fraction is given by the part of the grid box where the total moisture exceeds the saturation specific humidity (defined over ice if the local temperature is < 273.15 K, and over liquid water otherwise). See also Radiation for treatment of cloud-radiative interactions.

Precipitation

Large-scale precipitation forms in association with stratiform cloud (see Cloud Formation). For purposes of precipitation formation and the radiation calculations (see Radiation), the condensate is assumed to be liquid above 0 degrees C, and to be ice below -15 degrees C, with a liquid/ice fraction obtained by quadratic-spline interpolation for intermediate temperatures. The rate of conversion of cloud water into liquid precipitation is a nonlinear function of the large-scale cloud fraction and the cloud-mean liquid water content, following Sundqvist (1978 , 1981 ) and Golding (1986) . The precipitation of ice is a nonlinear function of the cloud-mean ice content, as deduced by Heymsfield (1977) . Liquid and frozen precipitation also form in subgrid-scale convection (see Convection).

Evaporation/sublimation of falling liquid/frozen precipitation are modeled after Kessler (1969) and Lin et al. (1983) . Frozen precipitation that falls to the surface defines the snowfall rate (see Snow Cover). For purposes of land hydrology, surface precipitation is assumed to be exponentially distributed over each land grid box, with fractional coverage of 0.5 for large-scale precipitation and 0.1 for convective precipitation. Cf. Smith and Gregory (1990) and Dolman and Gregory (1992) for further details.

Planetary Boundary Layer

Conditions within the PBL are typically represented by the first 5 levels above the surface (centered at about 997, 975, 930, 869, and 787 hPa for a surface pressure of 1000 hPa), where turbulent diffusion of momentum and cloud-conserved thermodynamic and moisture variables may occur (see Cloud Formation and Diffusion). The PBL top is defined either by the highest of these layers, or by the layer in which a modified bulk Richardson number (that incorporates buoyancy parameters for the cloud-conserved variables) exceeds a critical value of unity. Nonlocal mixing terms are included for heat and moisture. See also Surface Characteristics, Surface Fluxes, and Land Surface Processes.

Orography

For 0fix, 0cal and 6fix : Orography obtained from the U.S. Navy 10-minute resolution dataset (cf. Joseph 1980 ) is grid-box averaged, and is further smoothed with a 1-2-1 filter at latitudes poleward of 60 degrees. The orographic variances required by the gravity-wave drag parameterization are obtained from the same dataset (see Gravity-wave Drag).

For 21cal we use Peltier (1994).

Ocean

For fixed SSTs experiments, the model used is HADAM2.

For 0fix: The SST and sea ice datasets come from the Met. Office observational climatology MOHSST3 for 1951-1980.

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

For computed SSTs experiments, the model was coupled to an ocean model based on the Bryan-Cox primitive equation model (Bryan 1969a, Cox 1984), with 20 depht levels in the vertical. (HADAM2b)

For 0cal: SSTs and sea ice were calculated using the ocean model.

For 21 cal: The 21 kyr BP experiments with computed SSTs and sea ice should be performed in the same way as the CO2 experiments. Typically this means that a coupled atmosphere-ocean model will be used with prescribed present day ocean heat flux. This PMIP run uses Peltier ice sheets.

Sea Ice

For 0fix and 6fix, sea ice extents are prescribed monthly. The ice may occupy only a fraction of a grid box, and the effects of the remaining ice leads are accounted for in the surface roughness length, shortwave albedo and longwave emission, and turbulent eddy fluxes (see Surface Characteristics and Surface Fluxes). The spatially variable sea ice thickness is prescribed from climatological data. Snow falling on sea ice affects the surface albedo (see Surface Characteristics), but not the ice thickness or thermodynamic properties. Ice temperature is prognostically determined from a surface energy balance (see Surface Fluxes) that includes a conduction heat flux from the ocean below. Following Semtner (1976) , the conduction flux is proportional to the difference between the surface temperature of the ice and the subsurface ocean temperature (assumed to be fixed at the melting temperature of sea ice, or -1.8 degrees C), and the conduction flux is inversely proportional to the prescribed ice thickness. For computed SSTs experiments, the model was coupled with a sea ice model. The thermodynamics of the model is based on the zero-layer model of Stemtner (1976), and a simple parameterization of sea ice dynamics based on Bryan (1969) is included using the top layer ocean currents prescribed from values diagnosed in a full coupled ocean-atmosphere model.

Snow Cover

Surface snowfall is determined from the rate of frozen large-scale and convective precipitation in the lowest vertical layer. (Snowfall, like surface rainfall, is assumed to be distributed exponentially over each land grid box--see Precipitation.) On land only, prognostic snow mass is determined from a budget equation that accounts for accumulation, melting, and sublimation. Snow cover affects the roughness and heat conduction of the land surface, and it also alters the albedo of both land and sea ice (see Surface Characteristics). Snow melts when the temperature of the top soil/snow layer is > 0 degrees C, the snowmelt being limited by the total heat content of this layer. Snowmelt augments soil moisture, and sublimation of snow contributes to the surface evaporative flux over land. See also Surface Fluxes and Land Surface Processes.

Surface Characteristics

Surface types include land, ocean, sea ice, and permanent land ice. Sea ice may occupy a fraction of a grid box (see Sea Ice). On land, 24 different soil/vegetation types are specified from the 1 x 1-degree data of Wilson and Henderson-Sellers (1985) . The effects of these surface types on surface albedo (see below) and on surface thermodynamics and moisture (see Land Surface Processes) are treated via parameters derived by Buckley and Warrilow (1988) .

On each surface, roughness lengths are specified for momentum, for heat and moisture, and for free convective turbulence (which applies in cases of very light surface winds under unstable conditions). Over oceans, the roughness length for momentum is a function of surface wind stress (cf. Charnock 1955 ), but is constrained to be at least 10^-4 m; the roughness length for heat and moisture is a constant 10^-4 m, and it is 1.3 x 10^-3 m for free convective turbulence. Over sea ice, the roughness length is a constant 3 x 10^-3 m, but it is 0.10 m for that fraction of the grid box with ice leads (see Sea Ice). Over land, the roughness length is a function of vegetation and small surface irregularities; it is decreased as a linear function of snow cover, but is at least 5 x 10^-4 m. Cf. Smith (1990b) for further details.

The surface albedo of open ocean is a function of solar zenith angle. The albedo of sea ice varies between 0.60 and 0.85 as a linear function of the ice temperature above -5 degrees C, and it is also modified by snow cover. Where there is partial coverage of a grid box by sea ice, the surface albedo is given by the fractionally weighted albedos of sea ice and open ocean. Surface albedos of snow-free land are specified according to climatological vegetation and land use. Snow cover modifies the albedo of the land surface depending (exponentially) on the depth (following Hansen et al. 1983 ) and (linearly) on the snow temperature above -2 degrees C (following Warrilow et al. 1990 ), as well as on the background vegetation (following Buckley and Warrilow 1988 ). Cf. Ingram (1993) for further details.

Longwave emissivity is unity (blackbody emission) for all surfaces. Thermal emission from grid boxes with partial coverage by sea ice (see Sea Ice) is calculated from the different surface temperatures of ice and the open-ocean leads, weighted by the fractional coverage of each.

Surface Fluxes

Surface solar absorption is determined from the surface albedos, and longwave emission from the Planck equation with constant emissivity of 1.0 (see Surface Characteristics).

Turbulent eddy fluxes are formulated as bulk formulae in a constant-flux surface layer, following Monin-Obukhov similarity theory. The momentum flux is expressed in terms of a surface stress, and the heat and moisture fluxes in terms of cloud-conserved quantities (liquid/frozen water temperature and total water content) to account for effects of phase changes on turbulent exchanges (see Cloud Formation). These surface fluxes are solved by an implicit numerical method.

The surface atmospheric variables required for the bulk formulae are taken to be at the first level above the surface (at 997 hPa for a 1000 hPa surface pressure). (For diagnostic purposes, temperature and humidity at 1.5 m and the wind at 10 m are also estimated from the constant-flux assumption.) Following Louis (1979) , the drag/transfer coefficients in the bulk formulae are functions of stability (expressed as a bulk Richardson number) and roughness length, and the same transfer coefficient is used for heat and moisture. In grid boxes with fractional sea ice, surface fluxes are computed separately for the ice and lead fractions, but using mean drag and transfer coefficients obtained from linearly weighting the coefficients for ice and lead fractions (see Sea Ice).

The surface moisture flux also is a fraction beta of the local potential evaporation for a saturated surface. Over oceans, snow and ice, and where there is dew formation over land, beta is set to unity; otherwise, beta is a function of soil moisture and vegetation (see Land Surface Processes).

Above the surface layer, momentum and cloud-conserved variables are mixed vertically within the PBL by stability-dependent diffusion. For unstable conditions, cloud-conserved temperature and water variables are transported by a combination of nonlocal and local mixing. Cf. Smith (1990b, 1993 ) for further details. See also Diffusion and Planetary Boundary Layer.

Land Surface Processes

A vegetation canopy model (cf. Warrilow et al. 1986 and Shuttleworth 1988 ) includes effects of moisture condensation, precipitation interception, direct evaporation from wet leaves and from surface ponding, and evapotranspiration via root uptake of soil moisture. The fractional coverage and water-storage capacity of the canopy vary spatially by vegetation type (see Surface Characteristics), with a small storage added to represent surface ponding. The canopy intercepts a portion of the precipitation, which is exponentially distributed over each grid box (see Precipitation). Throughfall of canopy condensate and intercepted precipitation occurs in proportion to the degree of fullness of the canopy; evaporation from the wet canopy occurs in the same proportion, as a fraction of the local potential evaporation. Cf. Gregory and Smith (1990) for further details.

Soil moisture is predicted from a single-layer model with spatially nonuniform water-holding capacity; it is augmented by snowmelt, precipitation, and the throughfall of canopy condensate. This moisture infiltrates the soil at a rate depending on saturated soil hydraulic conductivity enhanced by effects of root systems that vary spatially by vegetation type (see Surface Characteristics). The noninfiltrated moisture is treated as surface runoff. Subsurface runoff from gravitational drainage is also parameterized as a function of spatially varying saturated soil hydraulic conductivity and of the ratio of soil moisture to its saturated value (cf. Eagleson 1978 ). Soil moisture is depleted by evaporation at a fraction beta of the local potential rate (see Surface Fluxes). The value of beta is obtained by combining the canopy evaporation efficiency (see above) with a moisture availability parameter that depends on the ratio of soil moisture to a spatially varying critical value, and of the ratio of the stomatal resistance to aerodynamic resistance (cf. Monteith 1965 ). The critical soil moisture and stomatal resistance vary spatially by soil/vegetation type (see Surface Characteristics). Cf. Gregory and Smith (1990) and Smith (1990b) for further details.

Soil temperature is predicted after Warrilow et al. (1986) from heat conduction in four layers. The depth of the topmost soil layer is given by the penetration of the diurnal wave, which depends on the spatially varying soil heat conductivity/capacity (see Surface Characteristics). The lower soil layers are, respectively, about 3.9, 14.1, and 44.7 times the depth of this top layer. The top boundary condition for heat conduction is the net downward surface energy balance (see Surface Fluxes), including the latent heat of fusion for snowmelt (see Snow Cover); the bottom boundary condition is zero heat flux. Heat insulation by snow (see Snow Cover) is modeled by reducing the thermal conductivity between the top two soil layers; however, subsurface moisture does not affect the thermodynamic properties of the soil. Cf. Smith (1990b) for further details.

Last update November 9, 1998. For further information, contact: Céline Bonfils (pmipweb@lsce.ipsl.fr )