2. Model runs

2.1 Conservation equations

For the inorganic carbon and radiocarbon, both passive tracers, the conservation equations carried in the model are

(1a)   d[DIC]/dt = L([DIC]) + Jv + J

and

(1b)   d[DIC14]/dt = L([DIC14]) - Lambda * [DIC14] + Jv14 + J14

where

• [DIC] is the model's concentration (moles/m^3) of total dissolved inorganic carbon;
• [DIC14] is the model's DIC-normalized concentration (also in moles/m^3) of total dissolved inorganic C-14 (see below);
• L is the 3-D transport operator, which represents effects due to advection, diffusion, and convection;
• Lambda is the radioactive decay constant for C-14 (ln(2) / 5730 year = 1.2097e-04 year^-1), converted to s^-1 using the number of seconds/year in your particular model;
• Jv is the "virtual" source-sink term representing the changes in surface [DIC] due to evaporation and precipitation, which must be accounted for because of the relatively high background concentration of [DIC];
• Jv14 is the "virtual" source-sink term for changes in surface [DIC14]due to evaporation and precipitation (E-P changes in background concentrations are of the same order as observed variability);
• J is the the source-sink term due to air-sea exchange of CO2; and
• J14 is the source-sink term due to air-sea exchange of 14CO2.

The source-sink terms Jv, Jv14, J, and J14 are added only as surface boundary conditions. That is they are equal to zero in all subsurface layers. These source-sink terms are equivalent to the fluxes, described below, divided by the surface layer thickness dz1.

    Jv = Fv/dz1

    Jv14 = Fv14/dz1

    J = F/dz1

    J14 = F14/dz1

2.2 Virtual flux (Fv)

In models where surface salinity is restored to observed values, this results in a surface flux of salt, not a surface flux of water as in the real world. Such surface salt fluxes are typically found in models with a rigid lid, and even in some models with a free surface (e.g., the OGCM from Louvain-la-Neuve). For simplicity, we categorize both classes of models as "rigid-lid-like". Conversely, non-rigid-lid-like models have a free surface and restore surface salinity by an equivalent flux of water leading to dilution or concentration (e.g., the MPI LSG model). Salinity in the latter type of free-surface model is conserved; E-P fluxes are taken into account by the velocity fields and thus do not need to be explicitly formulated in the transport model.

Yet for all rigid-lid-like models, we must explicitly take into account the concentration-dilution effect of E-P (Evaporation minus Precipitation), which changes surface [DIC] and [Alk]. Thus we add the virtual flux to the surface layer, each time step according to

(2a)   Fv = DICg * (E-P)

(2b)   Fv14 = DIC14g * (E-P)

where DICg and DIC14g are the model's globally averaged surface concentrations of DIC and DIC14, respectively. Both global averages must be computed at least once per year. For rigid-lid-like models with only salinity restoring, we suggest that (P - E) be computed as

(3)   P - E = (S - S')/ Sg * dz1 /Tau

where S' is the observed local salinity to which modeled local salinity S is being restored, Sg is the model's globally averaged surface salinity, dz1 is the top layer thickness, and Tau is the restoring time scale for salinity. For rigid-lid-like models which in addition include explicit P - E water fluxes, that term must of course also be added to eq (3).

2.3 Air-sea gas exchange fluxes (F and F14)

For simulations of DIC and DIC14, OCMIP-2 simulations will directly model the finite air-sea fluxes F and F14, respectively. Modelers must use the formulation for the standard OCMIP-2 air-to-sea flux,

(4a)   F = Kw (Csat - Csurf)

(4b)   F14 = Kw (14Csat - 14Csurf)

with

(5a)   Csat = alphaC * pCO2atm * P/Po

(5b)   14Csat = Csat * Ratm

where

• Kw is the CO2 gas transfer (piston) velocity [m/s] ;
• Csurf is the surface aqueous [CO2] concentration [mol/m^3], which is computed from the model's surface [DIC], T, S, and [Alk] (see section 2.5);
• 14Csurf is the surface ocean [14CO2] (see section 2.5);
• alphaC is the C solubility for water-vapor saturated air [mol/(m^3 * uatm)];
• pCO2atm is the partial pressure of CO2 in dry air at one atmosphere total pressure [in microatm], which is the same as the dry air mixing ratio of CO2 multiplied by 10^6 ;
• P is the total air pressure at sea level [atm], locally;
• Po is 1 atm; and
• Ratm is the normalized atmospheric ratio of C-14/C-12, which for our purposes we divide by the analogous ratio for the standard Rstd

(6)   Ratm = (1 + D14Catm/1000)

where D14Catm is the atmospheric Delta C-14, the fractionation corrected ratio of C-14/C-12, given in permil (see below).

Those familiar with C-14, may be surprised that in equation (6) we define Ratm, without multiplying the right hand term by Rstd (1.176e-12). Instead, we prefer to be able to compare [DIC14] to [DIC], directly, in order to simplify early interpretation and code verification. With the above formulation for the OCMIP equilibrium runs (where pCO2atm=278 ppm and D14Catm=0 permil), if both tracers are initialized identically, the only difference between units for the [DIC] and [DIC14] tracers will be due to radioactive decay. For the anthropogenic runs, there will also be contributions due to differences between atmospheric records for pCO2atm and D14Catm.

2.4 The Piston Velocity Kw

For simulations of DIC and DIC14, modelers must use the standard OCMIP-2 formulation for the piston velocity Kw for CO2. The monthly climatology of Kw, to be interpolated linearly in time by each modeling group, is computed with the following equation adapted from Wanninkhof (1992, eq. 3):

(7)    Kw = (1 - Fice) [Xconv * a *(u2 + v)] (Sc/660)**-1/2

where

• Fice is the fraction of the sea surface covered with ice, which varies from 0.0 to 1.0, and is given as monthly averages from the Walsh (1978) and Zwally et al. (1983) climatology (OCMIP-2 modelers must reset Fice values less than 0.2 to zero, after interpolation to their model grid)
• u2 is the instantaneous SSMI wind speed, averaged for each month, then squared, and subsequently averaged over th e same month of all years to give the monthly climatology. (see the OCMIP-1 README.satdat for further details);
• v is the variance of the instantaneous SSMI wind speed computed over one month temporal resolution And 2.5 degree spatial resolution, and subsequently averaged over the same month of all years to give the monthly climatology. Again, see the OCMIP-1 README.satdat for further details.
• a is the coefficient of 0.337, consistent with a piston velocity in cm/hr. We adjusted the coefficient a for OCMIP-2, in order to obtain Broecker et al.'s (1986) radiocarbon-calibrated, global CO2 gas exchange of 0.061 mol CO2 /(m^2 * yr * uatm), when using the satellite SSMI wind information (u2 + v) from Boutin and Etcheto (pers. comm.). Our computed value for a is similar to that determined by Wanninkhof (a = 0.31), who used a different wind speed data set and assumptions about wind speed variance; we use the observed variance.
• Xconv = 1/3.6e+05, is a constant factor to convert the piston velocity from [cm/hr] to [m/s]. This conversion factor is already included in the forcing field xKw, provided below.
• Sc is the Schmidt number which is to be computed using modeled SST, using the formulation from Wanninkhof (1992). The function scco2.f computes the Sc (unit-less) for CO2.

Practically speaking, to use equation (2) each group will interpolate the OCMIP-2 standard information to their own model grid. The standard information is provided by IPSL/LSCE as a monthly climatology on the 1 x 1 degree grid of Levitus (1982) in netCDF format (in file gasx_ocmip2.nc). Gridded variables in that file include

• the variable Fice,
• the second term, [Xconv * a * (u2 + v)], denoted as xKw [m/s]
• the mask Tmask (1 if ocean; 0 if land),
• the total atmospheric pressure at sea level P [atm]
• the longitude Lon at the center of each 1 x 1 degree grid box,
• the latitude Lat at the center of each 1 x 1 degree grid box.

For the variables Fice and xKw, continents on the 1 x 1 degree standard grid have been flooded with adjacent ocean values. Such an approach avoids discontinuities at land-sea boundaries during interpolation. See the Fortran program rgasx_ocmip2.f for an example of how to read the information in gasx_ocmip2.nc.gz into your interpolation routines. After compilation, to link and use rgasx_ocmip2.f, one must have already installed netCDF.

http://www.unidata.ucar.edu/packages/netcdf/

The file gasx_ocmip2.nc may also be inspected with software that uses netCDF format, such as ncdump or Ferret. Ferret will be used for some of the analysis during OCMIP-2. We encourage participants to become familiar with Ferret now

http://ferret.wrc.noaa.gov/Ferret/

After installation, one can visualize maps of the standard information in gasx_ocmip2.nc, by using the Ferret script vgasx_ocmip2.jnl.

After launching Ferret, simply issue the following command (at Ferret's "yes?" prompt)

yes? go vgasx_ocmip2.jnl

2.5 Oceanic and Atmospheric Components

Apart from Kw, there are a total of four other terms in equation (4a) and (4b) which require further development.

Ocean

The oceanic terms Csurf and 14Csurf [in mol/m^3] are not carried as tracers, so they must be computed each timestep to determine gas exchange

Csurf is the surface [CO2] concentration [mol/m^3], which is computed from the model's surface [DIC], T, S, and [Alk] through the equations and constants found in the subroutine co2calc.f. As input, we must provide alkalinity, which we determine as a normalized linear function of salinity.

(8)    [Alk] = Alkbar * S/Sbar

where [Alkbar] is 2310 microeq/kg and Sbar is the model's annual mean surface salinity, integrated globally (horizontally). Two other input arguments, both nutrient concentrations, are needed as input. Although accounting for both of their equilibria makes a difference, neither nutrient is included in the solubility pump run. Hence we take concentrations of both as being constant, equal to the global mean of surface observations: 0.5 micromol/kg for phosphate and 7.5 micromol/kg for silicate. Note that for the later OCMIP-2 run which includes the biological pump, we will use observed seasonal distributions of surface phosphate.

IMPORTANT: The carbonate chemistry subroutine co2calc.f was originally designed to require tracer input ([DIC], [Alk], [PO4], and [SiO2]) on a per mass basis (umol/kg); however, for OCMIP-2 co2calc.f has been modified to pass tracer concentrations on a per volume basis (mol/m^3), as carried in ocean models. To do so, we use the mean surface density of the ocean (1024.5 kg/m^3) as a constant conversion factor; we do NOT use model-predicted densities. For example, OCMIP-2 modelers should used SiO2 = 7.7e-03 mol/m^3 and PO4 = 5.1e-04 mol/m^3 as input arguments; again both are constant for the abiotic simulation. The output arguments co2star (Csurf) and dco2star (Csat - Csurf) are also returned in mol/m^3.

14Csurf is the surface ocean [14CO2], defined as

(9)    14Csurf = Csurf * Rocn,

where

(10)    Rocn = [DIC14]/[DIC].

Furthermore, for comparison to ocean measurements, we compute

(11)    D14Cocn = 1000*(Rocn - 1).

Following equation (4), we do not include Rstd when calculating D14Cocn in the model.

Atmosphere

The atmospheric components Csat and 14Csat in equations (4a) and (4b) are specified a priori via four remaining terms:

1. alphaC: The CO2 solubility alphaC is to be computed using modeled SST and SSS, both of which vary in time at each grid point. For OCMIP-2 we use the solubility formulation of Weiss (1974), corrected for the contribution of water vapor to the total pressure (Weiss and Price, 1980, Table IV for solubility in [mol/(l * atm)]). The solubility alphaC is calculated within the routine co2calc.f.
2. pCO2atm: For the Equilibrium run, pCO2atm is held constant at 278 ppm. For the anthropogenic perturbation, we define the equilibrium state as year 1765.0. Then for the Historical run, the model must be integrated until the end of 1999, following the observed record until 1990.5 ( splco2.dat) and IPCC scenario S650 ( stab.dat) until 2000.0 (Enting, 1994). That same scenario, Future run S650, will be continued from 2000.0 to 2300.0. Similarly, a second Future run CIS92A (see also cis92a.dat) will be run from 1990.0 to 2100.0, after initializing with model output from the Historical run in 1990.0. Additionally a Pulse run will be made, where preindustrial atmospheric CO2 is doubled and allowed to decline for 1000 years. Finally to eliminate effects due to model drift, we will make essentially two Control runs: (1) the first will be held to the same atmospheric boundary conditions as the Equilibrium run , carrying both DIC and DIC14 during 1765-2000 but only DIC from 2000-2300; (2) the second will be analogous to the Pulse run, made in forward mode for 1000 years, except that atmospheric CO2 will not be doubled on the first time step.
3. D14Catm: is atmospheric Delta C-14 [in permil]. For the Equilibrium run, D14Catm is held constant at 0 permil. For the Historical run, we define the equilibrium state as year 1765.0. Then the model must be integrated until the end of year 1999 following the observed record (Enting, 1994). The observed atmospheric C-14 record is given for three latitudinal bands:
   • 90S-20S
   • 20S-20N
   • 20N-90N
   There will be NO future or pulse simulations for C-14.
4. P: Is the total atmospheric pressure [atm] from the monthly mean climatology of Esbensen and Kushnir (1981). The latter, given originally on a 4 x 5 degree grid (latitude x longitude) in bars, is converted to atm by multiplying by (1/1.101325). Land and sea ice values in the original data set were filled with average values from adjacent ocean points. These monthly mean arrays were then linearly interpolated to the 1 x 1 degree grid of Levitus (see netCDF file gasx_ocmip2.nc).

Technical notes:

1. The ASCII file splco2.dat provides values of atmospheric pCO2 [in microatm], every half year, for the period from 1765.0 to 1990.5. Thereafter, there are two files used for future scenarios: for scenario S650, the ASCII file stab.dat provides half-year values of atmospheric pCO2 [in microatm] for the period from 1990.5 to 2300.5; for scenario CIS92A, the ASCII file cis92a.dat provides yearly values of atmospheric pCO2 for the period from 1990.5 to 2100.5. The subroutine read_co2atm.f reads atmospheric CO2 information from all three files.
2. The ASCII files c14nth.dat, c14equ.dat, and c14sth.dat provide mid-year values of atmospheric D14Catm [in permil] for the period from 1764.5 to 2000.0. See the subroutine read_c14atm.f
3. The Fortran subroutine c_interp.f temporally interpolates (linearly) both pCO2atm and D14Catm at a given timestep. That routine is called by the demonstration program try_c_interp.f, which spatially assigns D14Catm to the three latitudinal bands for C-14 (see above). Thus both routines together effect (1) temporal interpolation for both pCO2atm and D14Catm and (2) spatial "interpolation" for D14Catm as a function of latitude.