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**

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).

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**.

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

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

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 (

**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.

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

**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.**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.**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

**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:

- 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. - 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
- 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.