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


PMIP Documentation for GISS-IIP

Goddard Institute for Space Studies: Model GISS Model IIp (4x5 L9) 1997



PMIP Representative(s)

Dr. Robert S. Webb, NOAA Paleoclimatology Program, NGDC, 325 Broadway E/GC, Boulder, CO 80303, USA; Phone: +1-303 497 6967; Fax: +1-303 497 6513; e-mail:. rwebb@ngdc.noaa.gov

and:

Mr. Richard J. Healy, Marine Chemistry & Geochemistry, MS 25, Woods Hole Oceanographic Institution, Woods Hole, MA 02543, USA; Phone: +1-508 289 3514; Fax: +1-508 457 2193; e-mail:. rhealy@whoi.edu

World Wide Web URL: http://www.giss.nasa.gov

Model Designation

GISS Model IIp MY01 - (4x5L9) 1997 for 0fix

GISS Model IIp MY02 - (4x5L9) 1997 for 6fix

GISS Model IIp MY06M9 - (4x5L9) 1998 for 21fix
 

Model Identification for PMIP

GISS-IIP

PMIP run(s)

0fix, 6fix, 21fix

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

Differences relatives to AMIP and Hansen et al 1997: corrected calculation of insolation as a function of orbital forcing.
 

Model Lineage

 The GISS model used for the PMIP experiment is a modified version of Model II that is described by Hansen et al. (1983) . The current model (described by Hansen et al., 1997) differs from this predecessor principally in numerics and in the treatment of convection, planetary boundary layer (PBL), large-scale clouds, and ground hydrology. Table: Primary Differences between NASA-GISS Model II and Model II?
Process
Model II
Model II?
Impact
Boundary 
Layer		 infinite depth modified-Ekman layer, crude bulk formula for surface fluxes & exchange coeffients finite modified-Ekman layer drag & mixing coefficents based on similarity theory alters fluxes under unstable conditions; affects latitudinal precipitation distribution
Convection 50 percent of grid box mass elevated when moist static energy profile unstable variable mass flux proportional to instability, entraining & non-entraining plume, down draft mass 1/3 updrafts limits vertical mass exchange
Large Scale Clouds diagnostic cloud parameterization as the saturated fraction of the grid box, optical thickness based on type, altitude, thickness liquid water and ice prognostic variables, microphysical sources/sinks of cloud water, interactive optical thickness function of cloud particle size atmospheric dynamics associated with cloud radiative forcing, positive feedback on low latitude warming
Land Surface variable size two-layer bucket model specificed as a function of vegetation cover physically-based calculation of evapotranspiration, evaporation of intercepted precipitation, dew, and over bare soil, infiltration, soil water movement, & runoff more realistic partitioning of latent & sensible heating, stomatal control of transpiration
Heat and Moisture Advection second-order-centered differencing scheme quadratic upstream scheme preserve strong gradients, alter eddy available potential energy & moisture/ radiative forcing 
Momentum Advection second-order B grid scheme fourth-order B grid scheme propagate waves faster affecting Intrahemispheric transport & high-frequency synoptic variations
Surface Radiative fluxes radiative flux calculation for an average albedo of surface types in grid box individual radiative flux calculations for albedo of each surface type in grid box  

Model Documentation

Basic documentation of the model is provided by Hansen et al. (1983) , with subsequent changes summarized by Hansen et al. (1997). Changes in the convective scheme described by Del Genio and Yao (1993) , large-scale cloud parameterization by Del Genio et al. (1996), ground hydrology scheme by Rosenzweig and Abramopoulos et al. (1997), planetary boundary layer by Hartke and Rind (1997)

Numerical/Computational Properties

Horizontal Representation

Atmospheric mass and zonal and meridional velocity components are represented on a B-grid (cf. Arakawa 1972) , which conserves mass, kinetic energy (but not angular momentum) under advection, and enstrophy in the nondivergent limit. The prognostic variable for mass (see Atmospheric Dynamics) represents a mean value over a grid box. For heat and moisture prognostics, however, the linear gradients and second-order moments in three dimensions, and the three cross-term second-order moments are included in addition to the mean quantity for each grid box. (Potential enthalpy and water vapor are advected via a stable and accurate quadratic upstream scheme that utilizes these first- and second-order moments.)

Horizontal Resolution

4 x 5-degree latitude-longitude grid.

Dim_longitude*dim_latitude : 72*46

Vertical Domain

Surface to 10 hPa. For a surface pressure of 1000 hPa, the first vertical level above the surface is at about 975 hPa.

Vertical Representation

Finite-difference sigma coordinates up to 10 hPa, the model top for dynamics. Above 10 hPa the atmosphere interacts only radiatively with lower levels (i.e., its temperature profile here is determined solely by radiation).

Vertical Resolution

The model uses a sigma-coordinate in the vertical with 9 levels: 0.016897, 0.094938, 0.195759, 0.318899, 0.470418, 0.640124, 0.796957, 0.907372, 0.974264

For a surface pressure of 1000 hPa, 2 levels are below 800 hPa and 3 levels are above 200 hPa.

Computer/Operating System

The PMIP simulation was run on an SGI Origin 2000 running IRIX 6.4 operating system and FORTRAN-77 7.1.

The 21fix simulation was run on a Silicon Graphics multiprocessor R10000 cpu / R10010 fpu.

Computational Performance

For the PMIP experiment, about 3 minutes of SGI Origin 2000 computer time per simulated day or 18 hours per year.

Initialization

For the PMIP experiment, the model atmospheric state, soil moisture, and snow cover/depth are initialized for 1 January 1979 from a previous model simulation of December.

Time Integration Scheme(s)

Time integration is by a leapfrog scheme that is initiated each hour with an Euler-backward step. The dynamical time step is 7.5 minutes, while the physical source terms (except radiation) are updated hourly. (The hourly update of the PBL source term employs 4 successive 15-minute time steps.) Full radiation calculations are performed once every 5 hours.

Smoothing/Filling

The raw orography data is area-averaged (see Orography). An eighth-order Shapiro (1970) filter is used to smooth the surface pressure. Because the quadratic upstream scheme (see Horizontal Representation) does not guarantee positive-definite atmospheric moisture, the first- and second-order moments of water vapor mixing ratios are minimally reduced, if necessary, at each advective step to maintain non-negative values.

Sampling Frequency

For the PMIP simulation, Monthly averages of selected variables are saved as model history.

Dynamical/Physical Properties

Atmospheric Dynamics

Primitive-equation dynamics are expressed in terms of wind velocity, potential temperature, water-vapor mixing ratio, geopotential, pressure, and atmospheric mass (or PS-PT, where PS is the surface pressure and PT is a constant 10 hPa at the dynamical top--see Vertical Representation).

Diffusion

Neither horizontal nor vertical diffusion is explicitly modeled. However, there is horizontal transport of momentum associated with the convective mass flux (see Convection), and the quadratic upstream advective scheme (see Horizontal Representation) is weakly diffusive as well.

Gravity-wave Drag

A momentum drag proportional to air density and the square of the velocity is introduced in the top layer of the model (cf. Hansen et al. 1983) .

Solar Constant/Cycles

The solar constant is the prescribed value of 1367 W/(m2). The orbital parameters and seasonal insolation distribution are calculated after PMIP recommendations. A seasonal and diurnal cycle in solar forcing is simulated.

Chemistry

Carbon dioxide concentration is prescribed value of 345, 280 and 200 ppm for control, 6ka and 21ka, respectively. Ozone concentrations are prescribed as a function of season, height, and geographic location, with column abundances obtained from the climatology of London et al. (1976) . Radiative effects of water vapor, and of methane, nitrous oxide, nitric oxide, oxygen, and aerosols also are included (see Radiation).

CH4 and N2O concentrations are of 1643 ppb and 306 ppb for control, 650 ppb and 280 ppb for 6k, 350 ppb and 190 ppb for 21k.

Radiation

A correlated k-distribution method that is a generalization of the approach of Lacis and Hansen (1974) is used for both shortwave and longwave radiative calculations. The k-distribution for a given gas and frequency interval is obtained by least-squares fitting to calculations of line-by-line absorption for a range of temperatures and pressures by McClatchey et al. (1973) and Rothman (1981) .

The shortwave radiative fluxes are calculated after Lacis and Hansen (1974) , but with modifications to obtain accurate results at all solar zenith angles and optical thicknesses. Gaseous absorbers include water vapor, carbon dioxide, ozone, oxygen, nitrous oxide, and nitric oxide. Multiple scattering computations are made of 12 k-profiles, with strong line (exponential) absorption of the direct solar beam computed separately for water vapor, carbon dioxide, and oxygen. Absorption and scattering by aerosols are also included using radiative properties obtained from Mie calculations for the global aerosol climatology of Toon and Pollack (1976) . The spectral dependence of Mie parameters for clouds, aerosols, and Rayleigh scattering is specified in 6 intervals; these are superimposed on the 12 k-profiles to account for overlapping absorption.

In the longwave, the k-distribution method is used to model absorption by water vapor, carbon dioxide, ozone, nitrous oxide, nitric oxide, and methane (but with scattering effects neglected). A single k-distribution is specified for each gas, with 11 k-intervals for water vapor, 10 for carbon dioxide, and 4 for ozone.

Shortwave optical thickness of large-scale cloud is based on the prognostic cloud water path (see Cloud Formation). The required droplet effective radius is diagnosed from the cloud water content by assuming constant number concentration with different values for land/ocean cloud and liquid/ice cloud (cf. Del Genio et al. 1993) . The shortwave optical thickness of a convective cloud is proportional to its pressure depth. Cloud particle phase function and single-scattering albedo are functions of spectral interval, based on Mie computations for cloud droplet data of Squires (1958) and Hansen and Pollack (1970) . For purposes of the radiation calculations, partial cloud cover of a grid box is represented as full cloud cover that occurs for a percentage of the time (implemented via a random number generator--cf. Hansen et al. 1983) . Longwave effects of clouds are treated by an emissivity formulation, where the longwave cloud properties are self-consistent with the shortwave properties as a result of the application of Mie theory.

Convection

Convection is initiated only when the grid box is convectively unstable in the mean. Dry convection can occur below the condensation level if the moist static energy of a parcel exceeds that of the layer above. Moist convection is triggered if the moist static energy of a parcel exceeds the saturation energy value in the layer above and if, in addition, the implied lifting produces saturation; this level defines the cloud base. The convective cloud top is defined as the upper boundary of the highest vertical layer for which the cloud parcel is buoyant.

The amount of convective mass flux is obtained from a closure assumption that the cloud base is restored to neutral buoyancy relative to the next higher layer. The mass of the rising convective plume is changed by the entrainment of drier environmental air, with associated decreases in buoyancy. The entrainment rate is prescribed for an ensemble of two convective cloud types (entraining and non-entraining). Heating/cooling of the environment occurs through compensating environmental subsidence, detrainment of cloud air at cloud top, a convective-scale downdraft whose mass flux detrains into the cloud base layer, and evaporation of falling condensate (see Precipitation). (Latent heat release serves only to maintain cloud buoyancy.) The convective plume and subsiding environmental air transport gridscale horizontal momentum, under the assumption that exchanged air parcels carry with them the momentum of the layer of origin see Del Genio and Yao (1993).
 

Cloud Formation

Clouds result from either large-scale or convective condensation (see Convection). Condensation of cloud droplets is assumed to occur with respect to the saturated vapor pressure of water if the local temperature is >-35 degrees C, and with respect to the mixed-phase pseudo-adiabatic process of Sassen and Dodd (1989) at lower temperatures. Liquid droplets form if the cloud temperature is >-4 degrees C (-10 degrees C) over ocean (land), and ice forms at temperatures <-40 degrees C. At intermediate temperatures either phase may exist, with increasing probability of ice as temperature decreases.

The local convective cloud fraction is given by the ratio of convective mass flux to the total atmospheric mass of the grid box (see Convection and Atmospheric Dynamics). Large-scale clouds are predicted from a prognostic cloud water budget equation, where fractional cloudiness is an increasing function of relative humidity above a 60-percent threshold. (Upper-tropospheric convective condensate is also detrained into large-scale anvil cloud.) Cloud top entrainment instability is accounted for using a restrictive instability criterion. see Del Genio et al. (1996) for further details.

Precipitation

Precipitation can result either from large-scale or convective condensation (see Cloud Formation and Convection); the amount of condensate is computed from an iterative solution of the Clausius-Clapeyron equation. For large-scale condensation, the prognostic cloud water is converted to rainwater by autoconversion and accretion, and into snow by a seeder-feeder parameterization.

Evaporation of both the large-scale and convective condensate is modeled, with the residual falling to the surface as rain or snow (see Snow Cover). Evaporation of large-scale precipitation occurs in all unsaturated layers below its origin. The fraction of convective condensate that evaporates below cloud base is equated to the ratio of the convective mass flux to the total air mass of the grid box, while the fraction that evaporates above the cloud base is taken as half this value.

Planetary Boundary Layer

The PBL top is defined as the height of dry convection (see Convection). The surface wind velocity is determined using parameterizations of a drag law for the surface layer and a spiral layer above. Similarity theory is applied to calculate the turbulent transport coefficients for the spiral layer and to specify instability functions that determine the surface drag coefficient. The surface atmospheric temperature and water-vapor mixing ratio are obtained by equating surface fluxes computed by the bulk aerodynamic method to a diffusive flux from the surface layer into the layer above; the latter flux depends on a stability-dependent vertical diffusion coefficient and the depth of the PBL. See Hartke and Rind (1997) and also Surface Fluxes.
 

Orography

Orography is specified by area-averaging the 5-minute resolution topographic height data of National Geophysical Data Center (1988) on the 4 x 5-degree model grid (see Horizontal Resolution).

The topography used was based on Peltier?s modern and 21k topographic data. The procedure used to generate it took the difference between 21k and 0k and applied that to giss?s 0k which is based on the ETOPO5 data (5 minutes resolution). 106 meters are subtracted for the change in sea level between modern and 21k.

Ocean

For control: SSTs and sea ice extent are specified from the AMIP 10 year.

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

For 21 fix : The change in SST (LGM minus present-day) given by CLIMAP (1981) available at NGCD rather than the LGM absolute values in order to avoid differences due to differences in present day climatologies, was used. To obtain seasonally varying SSTs and sea ice edge from data for February and August, a simple sinusoidal variations, with extrems in February and August, is used.

The original 4x5 GISS CLIMAP LGM SSTs for Model II were generated from the CLIMAP LGM SSTs at 1x1 degree resolution from a tape provided by LDGO. In September 1995 we regridded the Model II data for the Model II? grid. The difference between the Model II and the Model II? 4x5 grid is a 2.5 shift degrees in longitude.
 

Sea Ice

For 0fix and 6fix, sea ice extents are prescribed monthly. The ice consists of two layers: an upper layer of constant thickness (0.1 m) and a lower layer whose thickness depends on the fractional ice coverage for each grid box and the number of months with some ice present. The total ice thickness also may be augmented by snow accumulation (see Snow Cover). The temperature profile in each ice layer is assumed to be a quadratic function of depth, with coefficients that are solved subject to six constraints. (These include consistency conditions on the temperatures and heat fluxes at the interfaces between the atmosphere and the ice, between the ice layers, and between the ice and the ocean; in addition, the mean temperature of each ice layer is set equal to its value from the previous hour.) After the ice temperature profiles are determined, the heat fluxes at the atmosphere-ice, ice-ice, and ice-ocean interfaces are updated. See also Surface Fluxes.

For 21fix, the sea-ice edge of CLIMAP data is used for February and August.

Snow Cover

Precipitation falls as snow if the temperature of the first vertical layer above the surface is <0 degrees C. Snow may accumulate on land, on continental ice, and on sea ice up to 0.1 m equivalent water, after which it augments the ice thickness (see Sea Ice). Snow depth is computed prognostically as the balance of snowfall, snowmelt, and sublimation (which contributes to the surface evaporative flux--see Surface Fluxes). Net heating of the snow surface raises the ground temperature to 0 degrees C; additional heating produces snowmelt, which contributes to soil moisture and affects the thermal properties of the soil (see Land Surface Processes). Snow also modifies the surface albedo (see Surface Characteristics). In some areas (e.g. the Tibetian Plateau) snow is retained year round, likely related to the fact that the snow albedo only changes from fresh snow when the entire snow pack and upper top layer of the soil reachs 0 deg C, hence a much larger heat capacity change than is realistic which seems to delay the changing the snow albedo due to snow aging/melting.

Surface Characteristics

Each grid box is assigned appropriate fractions of land and ocean, and part of the ocean fraction may be covered by ice (see Sea Ice). In addition, permanent ice sheets are specified for Antarctica, Greenland, and some Arctic islands. Vegetation type is a composite over each grid-box from 32 classifications distinguished in the 1 x 1-degree data of Matthews (1983 , 1984 ).

Over land, surface roughness is a fit to the data of Fiedler and Panofsky (1972) as a function of the standard deviation of the orography (see Orography). The maximum of this roughness and that of the local vegetation (including a "zero plane displacement" value for tall vegetation types--cf. Monteith 1973 determines the roughness over land. Over sea ice, the roughness length is a uniform 4.3 x 10-4 m, after Doronin (1969) . Over ocean, the surface roughness is a function of the momentum flux and is used to compute the neutral drag coefficient as well as the Stanton and Dalton numbers (see Surface Fluxes).

The surface albedo is prescribed for ocean, land ice, sea ice, and for eight different land surface types after the data of Matthews (1983 , 1984 ). The visible and near infrared albedos are distinguished, and seasonal variations in the albedo for vegetated surfaces are included. The albedo of snow-covered ground is the snow-free value modified by factors depending on snow depth, snow age, and vegetation masking depth, which varies with vegetation type. Ocean albedo is a function of surface wind speed and solar zenith angle after Cox and Munk (1956) .

The spectral dependence of longwave emissivity for deserts is included from data of Hovis and Callahan (1966) , and for snow and ice from data of Wiscombe and Warren (1980) . The emissivity of the ocean is a function of the surface wind speed and of the albedo.

Soil properties are specified from Webb et al. (1993)

Surface Fluxes

The absorbed surface solar flux is determined from albedos, and the surface longwave emission from the Planck function with spatially variable emissivities (see Surface Characteristics).

The surface wind stress is expressed as a product of air density, a drag coefficient, and the surface wind speed and velocity (see Planetary Boundary Layer). The surface drag coefficient is a function of both roughness length (see Surface Characteristics) and vertical stability.

Over land, the latent heat flux is computed separately for bare and vegetated surfaces (see Land Surface Processes), following Penman (1948) and Monteith (1981) . Over ocean and ice surfaces, the latent heat flux is expressed by a bulk formula that includes the product of the surface air density, a transfer coefficient, the surface wind speed, and the difference between the saturation value of mixing ratio at the ground temperature and the surface atmospheric value (see Planetary Boundary Layer).

Over land, the sensible heat flux is determined as the residual of the total net heat flux computed by the surface model (see Land Surface Processes) minus the latent heat flux. (The partitioning of the latent and sensible heat fluxes, or Bowen ratio, implies the surface temperature--see Land Surface Processes.) Over ocean and ice surfaces, the sensible heat flux is calculated from a bulk aerodynamic formula as a product of the surface air density and heat capacity, a transfer coefficient, a surface wind speed, and the difference between the skin temperature and the surface air temperature (see Planetary Boundary Layer). (The transfer coefficient is a stability-dependent function of the drag coefficient and is different from that used for the latent heat flux over ocean and ice.).

Land Surface Processes

Soil temperature is computed by solving a heat diffusion equation in six layers. The thickness of the top layer varies, but is approximately 0.1 m. The thicknesses of deeper layers increase geometrically, with the bottom boundary of the soil column at a nominal bedrock depth of 3.444 m. The upper boundary condition is the balance of surface energy fluxes (see Surface Fluxes); at the bottom boundary, zero net heat flux is specified. The thermal conductivity and heat capacity of the ground vary with snow cover, as well as soil moisture amount and phase.

Land-surface hydrology is treated after the physically based model of Abramopoulos et al. (1988) . The scheme includes a vegetation canopy, a composite over each grid box from the vegetation types of Matthews (1983 , 1984 ), that intercepts precipitation and dew. Evaporation from the wet canopy and from bare soil is treated, as well as soil-moisture loss from transpiration according to moisture availability and variable vegetation resistance and root density. Diffusion of moisture is predicted in the six soil layers, accounting for spatially variable composite conductivities and matric potentials that depend on soil type and moisture content. Infiltration of precipitation and snowmelt is explicitly calculated, with surface runoff occurring when the uppermost soil layer is saturated; underground runoff that depends on topographic slope is also included. See Rosenzweig and Abramopoulos et al. (1997) for additional details.

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