c_ ---------------------------------------------------------------------
c_ RCS lines preceded by "c_ "
c_ ---------------------------------------------------------------------
c_
c_ $Source: /www/html/ipsl/OCMIP/phase2/simulations/Abiotic/Cchem/RCS/co2calc.f,v $ 
c_ $Revision: 1.8 $
c_ $Date: 1999/07/16 11:40:33 $   ;  $State: Exp $
c_ $Author: orr $ ;  $Locker:  $
c_
c_ ---------------------------------------------------------------------
c_ $Log: co2calc.f,v $
c_ Revision 1.8  1999/07/16 11:40:33  orr
c_ Modifications by Keith Lindsay to fix inconsistency with common block
c_ "species" (not the same in "ta_iter_1.f").
c_ Also comment lines changed/added by J. Orr.
c_
c_ Revision 1.7  1999/04/26  13:04:54  orr
c_ Modified USAGE comment, to include new arguments
c_
c_ Revision 1.6  1999/04/14 12:55:52  orr
c_ Changed units for input arguments for tracers:
c_ formerly in mol/metric ton (T); now in mol/m^3.
c_ Used 1024.5 kg/m^3 as a constant conversion factor.
c_ Modelers can now pass tracers on a per volume basis, as carried in models.
c_
c_ Revision 1.5  1999/04/06 16:57:37  orr
c_ Changed calc of dpCO2 to account for diff atm pressure
c_
c_ Revision 1.4  1999/04/06 13:17:58  orr
c_ Added 2 output arguments: pCO2surf and dpCO2
c_ (see section 4 of Biotic HOWTO)
c_
c_ Revision 1.3  1999/04/05 15:59:11  orr
c_ Cleaned up comments regarding units
c_
c_ Revision 1.2  1999/04/04 02:35:00  orr
c_ Changed units for input and output tracers:
c_ previously in umol/kg; now in mol/T
c_ Can also pass input and output in mol/m^3 with little error.
c_
c_ Revision 1.1  1999/04/03  21:59:56  orr
c_ Initial revision
c_
c_ ---------------------------------------------------------------------
c_ 
      subroutine co2calc(t,s,dic_in,ta_in,pt_in,sit_in
     &                  ,phlo,phhi,ph,xco2_in,atmpres
     &                  ,co2star,dco2star,pCO2surf,dpco2)
C
C-------------------------------------------------------------------------
C SUBROUTINE CO2CALC
C
C PURPOSE
C	Calculate delta co2* from total alkalinity and total CO2 at
C temperature (t), salinity (s) and "atmpres" atmosphere total pressure. 
C
C USAGE
C       call co2calc(t,s,dic_in,ta_in,pt_in,sit_in
C    &                  ,phlo,phhi,ph,xco2_in,atmpres
C    &                  ,co2star,dco2star,pCO2surf,dpco2)
C
C INPUT
C	dic_in = total inorganic carbon (mol/m^3) 
C                where 1 T = 1 metric ton = 1000 kg
C	ta_in  = total alkalinity (eq/m^3) 
C	pt_in  = inorganic phosphate (mol/m^3) 
C	sit_in = inorganic silicate (mol/m^3) 
C	t      = temperature (degrees C)
C	s      = salinity (PSU)
C	phlo   = lower limit of pH range
C	phhi   = upper limit of pH range
C	xco2_in=atmospheric mole fraction CO2 in dry air (ppmv) 
C	atmpres= atmospheric pressure in atmospheres (1 atm==1013.25mbar)
C
C       Note: arguments dic_in, ta_in, pt_in, sit_in, and xco2_in are 
C             used to initialize variables dic, ta, pt, sit, and xco2.
C             * Variables dic, ta, pt, and sit are in the common block 
C               "species".
C             * Variable xco2 is a local variable.
C             * Variables with "_in" suffix have different units 
C               than those without.

C OUTPUT
C	co2star  = CO2*water (mol/m^3)
C	dco2star = delta CO2 (mol/m^3)
c       pco2surf = oceanic pCO2 (ppmv)
c       dpco2    = Delta pCO2, i.e, pCO2ocn - pCO2atm (ppmv)
C
C IMPORTANT: Some words about units - (JCO, 4/4/1999)
c     - Models carry tracers in mol/m^3 (on a per volume basis)
c     - Conversely, this routine, which was written by observationalists 
c       (C. Sabine and R. Key), passes input arguments in umol/kg  
c       (i.e., on a per mass basis)
c     - I have changed things slightly so that input arguments are in mol/m^3,
c     - Thus, all input concentrations (dic_in, ta_in, pt_in, and st_in) 
c       should be given in mol/m^3; output arguments "co2star" and "dco2star"  
c       are likewise in mol/m^3.

C FILES and PROGRAMS NEEDED
C	drtsafe
C	ta_iter_1
C
C--------------------------------------------------------------------------
C
	real invtk,is,is2
        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
	external ta_iter_1
C

c       ---------------------------------------------------------------------
C       Change units from the input of mol/m^3 -> mol/kg:
c       (1 mol/m^3)  x (1 m^3/1024.5 kg)
c       where the ocean's mean surface density is 1024.5 kg/m^3
c       Note: mol/kg are actually what the body of this routine uses 
c       for calculations.  
c       ---------------------------------------------------------------------
	permil = 1.0 / 1024.5
c       To convert input in mol/m^3 -> mol/kg 
	pt=pt_in*permil
	sit=sit_in*permil
	ta=ta_in*permil
	dic=dic_in*permil

c       ---------------------------------------------------------------------
C       Change units from uatm to atm. That is, atm is what the body of 
c       this routine uses for calculations.
c       ---------------------------------------------------------------------
	permeg=1.e-6
c       To convert input in uatm -> atm
	xco2=xco2_in*permeg
c       ---------------------------------------------------------------------
C
C*************************************************************************
C Calculate all constants needed to convert between various measured
C carbon species. References for each equation are noted in the code. 
C Once calculated, the constants are
C stored and passed in the common block "const". The original version of this
C code was based on the code by Dickson in Version 2 of "Handbook of Methods
C for the Analysis of the Various Parameters of the Carbon Dioxide System
C in Seawater", DOE, 1994 (SOP No. 3, p25-26). 
C
C Derive simple terms used more than once
C
	tk = 273.15 + t
	tk100 = tk/100.0
	tk1002=tk100*tk100
	invtk=1.0/tk
	dlogtk=log(tk)
	is=19.924*s/(1000.-1.005*s)
	is2=is*is
	sqrtis=sqrt(is)
	s2=s*s
	sqrts=sqrt(s)
	s15=s**1.5
	scl=s/1.80655
C
C f = k0(1-pH2O)*correction term for non-ideality
C
C Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)
C
	ff = exp(-162.8301 + 218.2968/tk100  +
     &		90.9241*log(tk100) - 1.47696*tk1002 +
     &		s * (.025695 - .025225*tk100 + 
     &		0.0049867*tk1002))
C
C K0 from Weiss 1974
C
	k0 = exp(93.4517/tk100 - 60.2409 + 23.3585 * log(tk100) +
     &		s * (.023517 - 0.023656 * tk100 + 0.0047036 * tk1002))

C
C k1 = [H][HCO3]/[H2CO3]
C k2 = [H][CO3]/[HCO3]
C
C Millero p.664 (1995) using Mehrbach et al. data on seawater scale 
C
	k1=10**(-1*(3670.7*invtk - 62.008 + 9.7944*dlogtk -
     &		0.0118 * s + 0.000116*s2))
C
	k2=10**(-1*(1394.7*invtk + 4.777 - 
     &		0.0184*s + 0.000118*s2))
C
C kb = [H][BO2]/[HBO2]
C
C Millero p.669 (1995) using data from Dickson (1990)
C
	kb=exp((-8966.90 - 2890.53*sqrts - 77.942*s +
     &		1.728*s15 - 0.0996*s2)*invtk +
     &		(148.0248 + 137.1942*sqrts + 1.62142*s) +
     &		(-24.4344 - 25.085*sqrts - 0.2474*s) *
     &		dlogtk + 0.053105*sqrts*tk)
C
C k1p = [H][H2PO4]/[H3PO4]
C
C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
C
	k1p = exp(-4576.752*invtk + 115.525 - 18.453 * dlogtk +
     &		(-106.736*invtk + 0.69171) * sqrts +
     &		(-0.65643*invtk - 0.01844) * s)
C
C k2p = [H][HPO4]/[H2PO4]
C
C DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))
C
	k2p = exp(-8814.715*invtk + 172.0883 - 27.927 * dlogtk +
     &		(-160.340*invtk + 1.3566) * sqrts +
     &		(0.37335*invtk - 0.05778) * s)
C
C------------------------------------------------------------------------
C k3p = [H][PO4]/[HPO4]
C
C DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)
C
	k3p = exp(-3070.75*invtk - 18.141 + 
     &		(17.27039*invtk + 2.81197) *
     &		sqrts + (-44.99486*invtk - 0.09984) * s)
C
C------------------------------------------------------------------------
C ksi = [H][SiO(OH)3]/[Si(OH)4]
C
C Millero p.671 (1995) using data from Yao and Millero (1995)
C
	ksi = exp(-8904.2*invtk + 117.385 - 19.334 * dlogtk +
     &		(-458.79*invtk + 3.5913) * sqrtis +
     &		(188.74*invtk - 1.5998) * is +
     &		(-12.1652*invtk + 0.07871) * is2 +
     &		log(1.0-0.001005*s))
C
C------------------------------------------------------------------------
C kw = [H][OH]
C
C Millero p.670 (1995) using composite data
C
	kw = exp(-13847.26*invtk + 148.9652 - 23.6521 * dlogtk +
     &		(118.67*invtk - 5.977 + 1.0495 * dlogtk) *
     &		sqrts - 0.01615 * s)
C
C------------------------------------------------------------------------
C ks = [H][SO4]/[HSO4]
C
C Dickson (1990, J. chem. Thermodynamics 22, 113)
C
	ks=exp(-4276.1*invtk + 141.328 - 23.093*dlogtk +
     &		(-13856*invtk + 324.57 - 47.986*dlogtk) * sqrtis +
     &		(35474*invtk - 771.54 + 114.723*dlogtk) * is -
     &		2698*invtk*is**1.5 + 1776*invtk*is2 +
     &		log(1.0 - 0.001005*s))
C
C------------------------------------------------------------------------
C kf = [H][F]/[HF]
C
C Dickson and Riley (1979) -- change pH scale to total
C
	kf=exp(1590.2*invtk - 12.641 + 1.525*sqrtis +
     &		log(1.0 - 0.001005*s) + 
     &		log(1.0 + (0.1400/96.062)*(scl)/ks))
C
C------------------------------------------------------------------------
C Calculate concentrations for borate, sulfate, and fluoride
C
C Uppstrom (1974)
	bt = 0.000232 * scl/10.811
C Morris & Riley (1966)
	st = 0.14 * scl/96.062
C Riley (1965)
	ft = 0.000067 * scl/18.9984
C*************************************************************************
C
C Calculate [H+] total when DIC and TA are known at T, S and 1 atm.
C The solution converges to err of xacc. The solution must be within
C the range x1 to x2.
C
C If DIC and TA are known then either a root finding or iterative method
C must be used to calculate htotal. In this case we use the Newton-Raphson
C "safe" method taken from "Numerical Recipes" (function "rtsafe.f" with
C error trapping removed).
C
C As currently set, this procedure iterates about 12 times. The x1 and x2
C values set below will accomodate ANY oceanographic values. If an initial
C guess of the pH is known, then the number of iterations can be reduced to
C about 5 by narrowing the gap between x1 and x2. It is recommended that
C the first few time steps be run with x1 and x2 set as below. After that,
C set x1 and x2 to the previous value of the pH +/- ~0.5. The current
C setting of xacc will result in co2star accurate to 3 significant figures
C (xx.y). Making xacc bigger will result in faster convergence also, but this
C is not recommended (xacc of 10**-9 drops precision to 2 significant figures).
C
C Parentheses added around negative exponents (Keith Lindsay)
C
	x1 = 10.0**(-phhi)
	x2 = 10.0**(-phlo)
	xacc = 1.e-10
	htotal = drtsafe(ta_iter_1,x1,x2,xacc)
C
C Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2, 
C ORNL/CDIAC-74, Dickson and Goyet, eds. (Ch 2 p 10, Eq A.49)
C
	htotal2=htotal*htotal
	co2star=dic*htotal2/(htotal2 + k1*htotal + k1*k2)
	co2starair=xco2*ff*atmpres
	dco2star=co2starair-co2star
	ph=-log10(htotal)

c
c       ---------------------------------------------------------------
cc      Add two output arguments for storing pCO2surf
cc      Should we be using K0 or ff for the solubility here?
c       ---------------------------------------------------------------
        pCO2surf = co2star / ff
        dpCO2    = pCO2surf - xco2*atmpres
C
C  Convert units of output arguments
c      Note: co2star and dco2star are calculated in mol/kg within this routine 
c      Thus Convert now from mol/kg -> mol/m^3
       co2star  = co2star / permil
       dco2star = dco2star / permil

c      Note: pCO2surf and dpCO2 are calculated in atm above. 
c      Thus convert now to uatm
       pCO2surf = pCO2surf / permeg
       dpCO2    = dpCO2 / permeg
C
	return
	end