! source file: /Users/oschlies/UVIC/master/testcase/updates/co2calc.F
! individual files downloaded from:
! http://www.ipsl.jussieu.fr/OCMIP/phase2/simulations/Abiotic/Cchem/

!_ ---------------------------------------------------------------------
!_ RCS lines preceded by "c_ "
!_ ---------------------------------------------------------------------
!_
!_ $Source: /www/html/ipsl/OCMIP/phase2/simulations/Abiotic/Cchem/RCS/
!_          co2calc.f,v $
!_ $Revision: 1.8 $
!_ $Date: 1999/07/16 11:40:33 $   ;  $State: Exp $
!_ $Author: orr $ ;  $Locker:  $
!_
!_ ---------------------------------------------------------------------
!_ $Log: co2calc.f,v $
!_ Revision 1.8  1999/07/16 11:40:33  orr
!_ Modifications by Keith Lindsay to fix inconsistency with common block
!_ "species" (not the same in "ta_iter_1.f").
!_ Also comment lines changed/added by J. Orr.
!_
!_ Revision 1.7  1999/04/26  13:04:54  orr
!_ Modified USAGE comment, to include new arguments
!_
!_ Revision 1.6  1999/04/14 12:55:52  orr
!_ Changed units for input arguments for tracers:
!_ formerly in mol/metric ton (T); now in mol/m^3.
!_ Used 1024.5 kg/m^3 as a constant conversion factor.
!_ Modelers can now pass tracers on a per volume basis, as carried in
!_ models.
!_
!_ Revision 1.5  1999/04/06 16:57:37  orr
!_ Changed calc of dpCO2 to account for diff atm pressure
!_
!_ Revision 1.4  1999/04/06 13:17:58  orr
!_ Added 2 output arguments: pCO2surf and dpCO2
!_ (see section 4 of Biotic HOWTO)
!_
!_ Revision 1.3  1999/04/05 15:59:11  orr
!_ Cleaned up comments regarding units
!_
!_ Revision 1.2  1999/04/04 02:35:00  orr
!_ Changed units for input and output tracers:
!_ previously in umol/kg; now in mol/T
!_ Can also pass input and output in mol/m^3 with little error.
!_
!_ Revision 1.1  1999/04/03  21:59:56  orr
!_ Initial revision
!_
!_ ---------------------------------------------------------------------
!_
      subroutine co2calc_SWS(t,s,dic_in,ta_in,pt_in,sit_in
     &                  ,phlo,phhi,ph,xco2_in,atmpres
     &                  ,co2star,dco2star,pCO2surf,dpco2)

!-----------------------------------------------------------------------
! subroutine CO2CALC

! PURPOSE
!        Calculate delta co2* from total alkalinity and total CO2 at
! temperature (t), salinity (s) and "atmpres" atmosphere total pressure.

! USAGE
!       call co2calc(t,s,dic_in,ta_in,pt_in,sit_in
!    &                  ,phlo,phhi,ph,xco2_in,atmpres
!    &                  ,co2star,dco2star,pCO2surf,dpco2)

! INPUT
!        dic_in  = total inorganic carbon (mol/m^3)
!                where 1 T = 1 metric ton = 1000 kg
!        ta_in   = total alkalinity (eq/m^3)
!        pt_in   = inorganic phosphate (mol/m^3)
!        sit_in  = inorganic silicate (mol/m^3)
!        t       = temperature (degrees C)
!        s       = salinity (PSU)
!        phlo    = lower limit of pH range
!        phhi    = upper limit of pH range
!        xco2_in =atmospheric mole fraction CO2 in dry air (ppmv)
!        atmpres = atmospheric pressure in atmospheres
!                 (1 atm==1013.25mbar)

!       Note: arguments dic_in, ta_in, pt_in, sit_in, and xco2_in are
!             used to initialize variables dic, ta, pt, sit, and xco2.
!             * Variables dic, ta, pt, and sit are in the common block
!               "species".
!             * Variable xco2 is a local variable.
!             * Variables with "_in" suffix have different units
!               than those without.

! OUTPUT
!        co2star  = CO2*water (mol/m^3)
!        dco2star = delta CO2 (mol/m^3)
!       pco2surf = oceanic pCO2 (ppmv)
!       dpco2    = Delta pCO2, i.e, pCO2ocn - pCO2atm (ppmv)

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

! FILES and PROGRAMS NEEDED
!        drtsafe
!        ta_iter_1

!-----------------------------------------------------------------------

        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

!       ----------------------------------------------------------------
!       Change units from the input of mol/m^3 -> mol/kg:
!       (1 mol/m^3)  x (1 m^3/1024.5 kg)
!       where the ocean's mean surface density is 1024.5 kg/m^3
!       Note: mol/kg are actually what the body of this routine uses
!       for calculations.
!       ----------------------------------------------------------------
        permil = 1.0 / 1024.5
!       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

!       ----------------------------------------------------------------
!       Change units from uatm to atm. That is, atm is what the body of
!       this routine uses for calculations.
!       ----------------------------------------------------------------
        permeg=1.e-6
!       To convert input in uatm -> atm
        xco2=xco2_in*permeg
!       ----------------------------------------------------------------

!***********************************************************************
! Calculate all constants needed to convert between various measured
! carbon species. References for each equation are noted in the code.
! Once calculated, the constants are
! stored and passed in the common block "const". The original version of
! this code was based on the code by Dickson in Version 2 of "Handbook
! of Methods for the Analysis of the Various Parameters of the Carbon
! Dioxide System in Seawater", DOE, 1994 (SOP No. 3, p25-26).

! Derive simple terms used more than once

        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

! f = k0(1-pH2O)*correction term for non-ideality

! Weiss & Price (1980,Mar. Chem.,8,347-359; Eq 13 with table 6 values)

        ff = exp(-162.8301 + 218.2968/tk100  +
     &                90.9241*log(tk100) - 1.47696*tk1002 +
     &                s * (.025695 - .025225*tk100 +
     &                0.0049867*tk1002))

! K0 from Weiss 1974

        k0 = exp(93.4517/tk100 - 60.2409 + 23.3585 * log(tk100) +
     &                s * (.023517 - 0.023656 * tk100 + 0.0047036
     &                * tk1002))

! k1 = [H][HCO3]/[H2CO3]
! k2 = [H][CO3]/[HCO3]

! Millero p.664 (1995) using Mehrbach et al. data on seawater scale

        k1=10**(-1*(3670.7*invtk - 62.008 + 9.7944*dlogtk -
     &                0.0118 * s + 0.000116*s2))

        k2=10**(-1*(1394.7*invtk + 4.777 -
     &                0.0184*s + 0.000118*s2))

! kb = [H][BO2]/[HBO2]

! Millero p.669 (1995) using data from Dickson (1990)

        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)

! k1p = [H][H2PO4]/[H3PO4]

! DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)

        k1p = exp(-4576.752*invtk + 115.525 - 18.453 * dlogtk +
     &                (-106.736*invtk + 0.69171) * sqrts +
     &                (-0.65643*invtk - 0.01844) * s)

! k2p = [H][HPO4]/[H2PO4]

! DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))

        k2p = exp(-8814.715*invtk + 172.0883 - 27.927 * dlogtk +
     &                (-160.340*invtk + 1.3566) * sqrts +
     &                (0.37335*invtk - 0.05778) * s)

!-----------------------------------------------------------------------
! k3p = [H][PO4]/[HPO4]

! DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)

        k3p = exp(-3070.75*invtk - 18.141 +
     &                (17.27039*invtk + 2.81197) *
     &                sqrts + (-44.99486*invtk - 0.09984) * s)

!-----------------------------------------------------------------------
! ksi = [H][SiO(OH)3]/[Si(OH)4]

! Millero p.671 (1995) using data from Yao and Millero (1995)

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

!-----------------------------------------------------------------------
! kw = [H][OH]

! Millero p.670 (1995) using composite data

        kw = exp(-13847.26*invtk + 148.9652 - 23.6521 * dlogtk +
     &                (118.67*invtk - 5.977 + 1.0495 * dlogtk) *
     &                sqrts - 0.01615 * s)

!-----------------------------------------------------------------------
! ks = [H][SO4]/[HSO4]

! Dickson (1990, J. chem. Thermodynamics 22, 113)

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

!-----------------------------------------------------------------------
! kf = [H][F]/[HF]

! Dickson and Riley (1979) -- change pH scale to total

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

!-----------------------------------------------------------------------
! Calculate concentrations for borate, sulfate, and fluoride

! Uppstrom (1974)
        bt = 0.000232 * scl/10.811
! Morris & Riley (1966)
        st = 0.14 * scl/96.062
! Riley (1965)
        ft = 0.000067 * scl/18.9984
!***********************************************************************

! Calculate [H+] total when DIC and TA are known at T, S and 1 atm.
! The solution converges to err of xacc. The solution must be within
! the range x1 to x2.

! If DIC and TA are known then either a root finding or iterative method
! must be used to calculate htotal. In this case we use the
! Newton-Raphson "safe" method taken from "Numerical Recipes" (function
! "rtsafe.f" with error trapping removed).

! As currently set, this procedure iterates about 12 times. The x1 and
! x2 values set below will accomodate ANY oceanographic values. If an
! initial guess of the pH is known, then the number of iterations can be
! reduced to about 5 by narrowing the gap between x1 and x2. It is
! recommended that the first few time steps be run with x1 and x2 set as
! below. After that, set x1 and x2 to the previous value of the pH +/-
! ~0.5. The current setting of xacc will result in co2star accurate to 3
! significant figures (xx.y). Making xacc bigger will result in faster
! convergence also, but this is not recommended (xacc of 10**-9 drops
! precision to 2 significant figures).

! Parentheses added around negative exponents (Keith Lindsay)

        x1 = 10.0**(-phhi)
        x2 = 10.0**(-phlo)
        xacc = 1.e-10
        htotal = drtsafe(ta_iter_1,x1,x2,xacc)

! Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2,
! ORNL/CDIAC-74, Dickson and Goyet, eds. (Ch 2 p 10, Eq A.49)

        htotal2=htotal*htotal
        co2star=dic*htotal2/(htotal2 + k1*htotal + k1*k2)
        co2starair=xco2*ff*atmpres
        dco2star=co2starair-co2star
        ph=-log10(htotal)

!       ---------------------------------------------------------------
!c      Add two output arguments for storing pCO2surf
!c      Should we be using K0 or ff for the solubility here?
!       ---------------------------------------------------------------
        pCO2surf = co2star / ff
        dpCO2    = pCO2surf - xco2*atmpres

!  Convert units of output arguments
!      Note: co2star and dco2star are calculated in mol/kg within this
!      routine. Thus Convert now from mol/kg -> mol/m^3
       co2star  = co2star / permil
       dco2star = dco2star / permil

!      Note: pCO2surf and dpCO2 are calculated in atm above.
!      Thus convert now to uatm
       pCO2surf = pCO2surf / permeg
       dpCO2    = dpCO2 / permeg

        return
        end

!_ ---------------------------------------------------------------------
!_ RCS lines preceded by "c_ "
!_ ---------------------------------------------------------------------
!_
!_ $Source: /www/html/ipsl/OCMIP/phase2/simulations/Abiotic/Cchem/RCS/
!_          drtsafe.f,v $
!_ $Revision: 1.1 $
!_ $Date: 1999/04/03 22:00:42 $   ;  $State: Exp $
!_ $Author: orr $ ;  $Locker:  $
!_
!_ ---------------------------------------------------------------------
!_ $Log: drtsafe.f,v $
!_ Revision 1.1  1999/04/03 22:00:42  orr
!_ Initial revision
!_
!_ ---------------------------------------------------------------------
!_
      real FUNCTION DRTSAFE(FUNCD,X1,X2,XACC)

!     File taken from Numerical Recipes. Modified  R.M.Key 4/94

      MAXIT=100
      call FUNCD(X1,FL,DF)
      call FUNCD(X2,FH,DF)
      if (FL .lt. 0.0) then
        XL=X1
        XH=X2
      else
        XH=X1
        XL=X2
        SWAP=FL
        FL=FH
        FH=SWAP
      endif
      DRTSAFE=.5*(X1+X2)
      DXOLD=ABS(X2-X1)
      DX=DXOLD
      call FUNCD(DRTSAFE,F,DF)
      do 100, J=1,MAXIT
        if (((DRTSAFE-XH)*DF-F)*((DRTSAFE-XL)*DF-F) .ge. 0.0 .or.
     &              ABS(2.0*F) .gt. ABS(DXOLD*DF)) then
          DXOLD=DX
          DX=0.5*(XH-XL)
          DRTSAFE=XL+DX
          if (XL .eq. DRTSAFE) return
        else
          DXOLD=DX
          DX=F/DF
          TEMP=DRTSAFE
          DRTSAFE=DRTSAFE-DX
          if (TEMP .eq. DRTSAFE) return
        endif
        if (ABS(DX) .lt. XACC) return
        call FUNCD(DRTSAFE,F,DF)
        if (F .lt. 0.0) then
          XL=DRTSAFE
          FL=F
        else
          XH=DRTSAFE
          FH=F
        endif
  100 continue
      return
      end

!_ ---------------------------------------------------------------------
!_ RCS lines preceded by "c_ "
!_ ---------------------------------------------------------------------
!_
!_ $Source: /www/html/ipsl/OCMIP/phase2/simulations/Abiotic/Cchem/RCS/
!_          ta_iter_1.f,v $
!_ $Revision: 1.2 $
!_ $Date: 1999/09/01 17:55:41 $   ;  $State: Exp $
!_ $Author: orr $ ;  $Locker:  $
!_
!_ ---------------------------------------------------------------------
!_ $Log: ta_iter_1.f,v $
!_ Revision 1.2  1999/09/01 17:55:41  orr
!_ Fixed sign error in dfn/dx following remarks of C. Voelker
!_ (10/Aug/1999)
!_
!_ Revision 1.1  1999/04/03 22:00:42  orr
!_ Initial revision
!_
!_ ---------------------------------------------------------------------
!_
        subroutine ta_iter_1(x,fn,df)
        real k12,k12p,k123p
        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

! This routine expresses TA as a function of DIC, htotal and constants.
! It also calculates the derivative of this function with respect to
! htotal. It is used in the iterative solution for htotal. In the call
! "x" is the input value for htotal, "fn" is the calculated value for TA
! and "df" is the value for dTA/dhtotal

        x2=x*x
        x3=x2*x
        k12 = k1*k2
        k12p = k1p*k2p
        k123p = k12p*k3p
        c = 1.0 + st/ks
        a = x3 + k1p*x2 + k12p*x + k123p
        a2=a*a
        da = 3.0*x2 + 2.0*k1p*x + k12p
        b = x2 + k1*x + k12
        b2=b*b
        db = 2.0*x + k1

!     fn = hco3+co3+borate+oh+hpo4+2*po4+silicate+hfree+hso4+hf+h3po4-ta

        fn = k1*x*dic/b +
     &             2.0*dic*k12/b +
     &             bt/(1.0 + x/kb) +
     &             kw/x +
     &             pt*k12p*x/a +
     &             2.0*pt*k123p/a +
     &             sit/(1.0 + x/ksi) -
     &             x/c -
     &             st/(1.0 + ks/x/c) -
     &             ft/(1.0 + kf/x) -
     &             pt*x3/a -
     &             ta

!        df = dfn/dx

        df = ((k1*dic*b) - k1*x*dic*db)/b2 -
     &             2.0*dic*k12*db/b2 -
     &             bt/kb/(1.0+x/kb)**2. -
     &             kw/x2 +
     &             (pt*k12p*(a - x*da))/a2 -
     &             2.0*pt*k123p*da/a2 -
     &             sit/ksi/(1.0+x/ksi)**2. -
     &             1.0/c +
     &             st*(1.0 + ks/x/c)**(-2.)*(ks/c/x2) +
     &             ft*(1.0 + kf/x)**(-2.)*kf/x2 -
     &             pt*x2*(3.0*a-x*da)/a2

        return
        end