Dear OCMIPpers,

Attached you will find the source code for the carbon chemistry calculations
as per Ray's recent message(s). The carbon constants are Mehrbach as refit by
Dickson and Millero (1987) and taken from Millero (1995, GCA, 59) on 
the seawater pH scale (see also CO2 handbook). The subroutines are all
documented within the source files. The files are:

test.r		test program showing units conversions, etc.
co2calc.f	main driver subroutine
ta_iter_1.f	subroutine which defines alk in terms of tco2 and htotal
drtsafe.f	"safe" root finding function

NOTES:
In co2calc the initial pH range is set to 6 to 9 which should cover anything
that your models generate (in terms of combined t,s,tco2,alk). With this
range, it takes an average of 12 iterations for the solution to converge.
The number of iterations can be dropped to around 4 IF these bounds can
be narrowed to a span of about 1 to 1.5 pH units. Without knowing exactly
how your individual programs work, we didn't try to narrow the bounds, but
doing so should speed up the calculations. One example of how you might do 
this would be to use the pH field from the previous time step. If at the 
last time step the pH (called htotal in carbon.inc) was 7.8, you could 
probably set the search bounds to 7.8+/-0.5 (values for x1 and x2 in 
co2calc). All error trapping has been disabled in this code, so do not set 
x1 and x2 so narrow that the calculated htotal (with the new ta and tco2) 
is outside the bounds of x1 and x2. We pass ph from co2calc.f so it can be
saved if wanted. We also calculate K0 in case you want to calculate the fCO2
of the water (FCO2=co2star*10**6/k0).

We have included a test program that runs through a range of values for error
checking. Note that all of co2calc.f expects units of umol/kg. These values 
are converted to moles/kg for the calculations. If you enter the wrong units 
it will effect your answers.

Based on our discussions at the US-OCMIP meeting in Boulder, we have designed
the program to output delta co2*=(co2*air - co2*water). 
co2*air is a function of xco2 (mole fraction in dry air) and barometric
pressure. The net flux is simply 

Flux = kw * delta CO2*
	where kw is the piston velocity 
	and delta CO2* is the air-sea disequilibrium (mole/kg).


OTHER:
There are a few type specifications, otherwise the variable names are
normal fortran (a-h,o-z = real).
The code can be very easily modified to work with array or matrix input (except
for the iteration step).
Precision is more than sufficient even on a SUN.
Program test.r was not translated, but an example output file (test.out)
is included in the mailing.
If t, and s are changing slowly enough, you may want to bypass recalculation
of the constants at each cell or time step.
Try to devise a rational way to limit the range of x1 and x2 (pH) for a really
significant speed increase.

3/9/99
Code modified to allow for atmospheric pressures different from 1 atmosphere
	via the variable atmpres (in test.r and co2calc.f).
Code modified so that input values for pt, sit, dic, ta are returned from
	co2calc (via common block) unchanged.