program test
#
# input and output values for dic and ta are umol/kg, permil for pco2
#
        real k0,k1,k2,kw,kb,ks,kf,k1p,k2p,k3p,ksi
	common /const/k0,k1,k2,kw,kb,ks,kf,k1p,k2p,k3p,ksi,ff,htotal
	common /species/bt,st,ft,sit,pt,dic,ta
phlo=6.
phhi=9.
open(11,file="test.out",form="FORMATTED")
for(i=1;i<=6;i=i+1)
	{
	t=float(i)*5.0
	for(j=0;j<=8;j=j+1)
		{
		s=float(j)*.5+34.
		write(11,100)t,s
		for(k=1;k<=8;k=k+1)
			{
			for(m=1;m<=8;m=m+1)
				{
#
# Initialize silicate, phosphate, xco2 and atmospheric pressure
#  for this example. GCM code will presumably be carrrying these
#  variables around in some manner or another. Units must be correct
#  as described in co2calc.f. For this example these values must be set
# HERE since they are changed inside co2calc and not reset. If this is
# cumbersome then add code to reset the units inside co2calc.f
#
				sit=0.0
				pt=1.0
				xco2=364.0
				atmpres=1.0
				ta=2000+float(m-1)*25.
				dic=1900+float(k-1)*25.
				call co2calc(t,s,
					phlo,phhi,ph,xco2,atmpres,
					co2star,dco2star)
				fco2=co2star*10**6/k0
#
# convert units and print
#
				write(11,200)dic,ta,fco2,
					co2star*10**6,dco2star*10**6,ph				}
			}
		}
	}
close(11)
100	format(2f8.2)
200	format(6f8.2)
stop
end