#include "options.inc"

#define use_fast_drtsafe_routine

#ifdef drive_test_co2
c return as well equilibrium DIC
c#define calculate_C_eq
c return also H+ concentration
c#define calc_htotal_also
#endif

      module pco2_module
c-----------------------------------------------------------------------
c   Module for surface pCO2 calculation (see below)
c-----------------------------------------------------------------------
      implicit none
      integer, parameter :: imax = 2000 ! max. vector length
      real, dimension(imax) :: k0,k1,k2,kw,kb,ks,kf,k1p,k2p,k3p,ksi,ff  ! Mehrbach constants
      real, dimension(imax) :: bt,st,ft   ! borate sulfate and fluor ??
      real, dimension(imax) :: pt,dic,ta  ! phosphate, DIC and alkalinity
      real, parameter :: rho0 = 1024.0    ! unit conversion factor
      real, parameter :: rrho02 = 1.0/rho0**2
      real, parameter :: sit=0.0077/rho0   ! silicate in mol/m**3, taken as constant
      real, parameter :: permeg=1.e-6, permil = 1.0/rho0 ! unit conversion
c      real, parameter :: atmpres = 1.0  ! atmospheric pressure in atm = 1013.25 mbar
      integer :: i
#ifdef  calculate_C_eq
      real, dimension(imax) :: co2star_eq
#endif
      end module pco2_module


      subroutine co2calc(t,s_in,dic_in,nit,slp,imt,
     &                   pco2_air,dco2,pco2_water
#ifdef  calculate_C_eq
     &                      , dic_eq
#endif
#ifdef calc_htotal_also
     &                      , htotal_out
#endif
     &)
c-------------------------------------------------------------------------
c   Calculate pco2* from total alkalinity and total CO2 at
c   temperature (t), salinity (s) and "atmpres" atmosphere total pressure. 
c--------------------------------------------------------------------------
       use timing_module
       use pco2_module
       implicit none
c      INPUT
       integer :: imt         ! vector length
       real   :: t(imt)        ! temperature deg C
       real   :: s_in(imt)     ! salinity in (PSU-35)/1000
       real   :: dic_in(imt)   ! diss. inorg. carbon  in mmol/m^3
       real   :: nit(imt)      ! nitrate in mmol/m^3
       real   :: slp(imt)      ! Sea level pressure in mbar
       real   :: pco2_air      ! atmospheric mole fraction CO2 in dry air (ppmv) 
c      OUTPUT
       real   :: dco2(imt)     ! difference in CO2 at sea surface and saturation CO2 in mol/m^3
       real   :: pco2_water(imt) ! partial pressure of CO2 at sea surface  (ppmv)
#ifdef  calculate_C_eq
       real   :: dic_eq(imt)   ! DIC in equilibrium mol/m^3
#endif
#ifdef calc_htotal_also
       real   :: htotal_out(imt)   ! H+ concentration
#endif

c      local work variables
       real :: tk,tk100,tk1002,invtk,dlogtk
       real :: is,is2,sqrtis,s2,sqrts,s15,scl
       real :: htotal(imt),s(imt),htotal2,co2star,pco2

c      stuff for the iteration
       real, external :: dic_iter
#ifdef  calculate_C_eq
       real, external :: dic_iter_C_eq
#endif

#ifdef use_fast_drtsafe_routine
       real :: drtsafe, uuu, x1, x2, fl, fh, df, swap, dx ,dxold, xl, xh
       real :: f, xacc, temp
       integer :: maxit, j
#else
       real, external :: drtsafe
#endif
       real :: lower_bound = 1.0e-6
       real :: upper_bound = 1.0e-9
       real :: accuracy = 1.e-10
       real :: atmpres

c       ---------------------------------------------------------------------

c        call tic('co2_flux')
        if (imt > imax) then
          print*,' imt larger than imax = ',imax
          print*,' increase imax in pco2_module'
          call halt_stop('in co2calc')
        endif

        do i=1,imt
         s(i) = s_in(i)*1000.0+35.0
         is2=s(i)*s(i)
         dic(i) = dic_in(i)          *rrho02
c         ta(i)  = (2310.0*s(i)/35.0) *rrho02  ! old OMIP fit
         ta(i) = ( -0.5487*is2*s(i) + 59.919*is2   ! new fit from K. Friis for North Atl.
     &            - 2122.3*S(i) + 26722 )*rrho02
        
	 pt(i)  = nit(i)/16.0        *rrho02
        enddo

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

        do i=1,imt
         tk = 273.15 + t(i)
         tk100 = tk/100.0
         tk1002=tk100*tk100
         invtk=1.0/tk
         dlogtk=log(tk)
         is=19.924*s(i)/(1000.-1.005*s(i))
         is2=is*is
         sqrtis=sqrt(is)
         s2=s(i)*s(i)
         sqrts=sqrt(s(i))
         s15=s(i)**1.5
         scl=s(i)/1.80655
C        f = k0(1-pH2O)*correction term for non-ideality
C        Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)
         ff(i) = exp(-162.8301 + 218.2968/tk100  +
     &		90.9241*log(tk100) - 1.47696*tk1002 +
     &		s(i)*(.025695 - .025225*tk100 + 
     &		0.0049867*tk1002))
C        K0 from Weiss 1974
         k0(i) = exp(93.4517/tk100 - 60.2409 + 23.3585 * log(tk100) +
     &		s(i)*(.023517 - 0.023656 * tk100 + 0.0047036 * tk1002))
C        k1 = [H][HCO3]/[H2CO3]
C        k2 = [H][CO3]/[HCO3]
C        Millero p.664 (1995) using Mehrbach et al. data on seawater scale 
         k1(i)=10**(-1*(3670.7*invtk - 62.008 + 9.7944*dlogtk -
     &		0.0118 * s(i) + 0.000116*s2))
         k2(i)=10**(-1*(1394.7*invtk + 4.777 - 
     &		0.0184*s(i) + 0.000118*s2))
C        kb = [H][BO2]/[HBO2]
C        Millero p.669 (1995) using data from Dickson (1990)
         kb(i)=exp((-8966.90 - 2890.53*sqrts - 77.942*s(i) +
     &		1.728*s15 - 0.0996*s2)*invtk +
     &		(148.0248 + 137.1942*sqrts + 1.62142*s(i)) +
     &		(-24.4344 - 25.085*sqrts - 0.2474*s(i)) *
     &		dlogtk + 0.053105*sqrts*tk)
C        k1p = [H][H2PO4]/[H3PO4]
C        DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
         k1p(i) = exp(-4576.752*invtk + 115.525 - 18.453 * dlogtk +
     &		(-106.736*invtk + 0.69171) * sqrts +
     &		(-0.65643*invtk - 0.01844) * s(i))
C        k2p = [H][HPO4]/[H2PO4]
C        DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))
         k2p(i) = exp(-8814.715*invtk + 172.0883 - 27.927 * dlogtk +
     &		(-160.340*invtk + 1.3566) * sqrts +
     &		(0.37335*invtk - 0.05778) * s(i) )
C        k3p = [H][PO4]/[HPO4]
C        DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)
         k3p(i) = exp(-3070.75*invtk - 18.141 + 
     &		(17.27039*invtk + 2.81197) *
     &		sqrts + (-44.99486*invtk - 0.09984) * s(i))
C        ksi = [H][SiO(OH)3]/[Si(OH)4]
C        Millero p.671 (1995) using data from Yao and Millero (1995)
         ksi(i) = 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(i)))
C        kw = [H][OH]
C        Millero p.670 (1995) using composite data
         kw(i) = exp(-13847.26*invtk + 148.9652 - 23.6521 * dlogtk +
     &		(118.67*invtk - 5.977 + 1.0495 * dlogtk) *
     &		sqrts - 0.01615 * s(i))
C        ks = [H][SO4]/[HSO4]
C        Dickson (1990, J. chem. Thermodynamics 22, 113)
         ks(i)=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(i)))
C        kf = [H][F]/[HF]
C        Dickson and Riley (1979) -- change pH scale to total
         kf(i)=exp(1590.2*invtk - 12.641 + 1.525*sqrtis +
     &		log(1.0 - 0.001005*s(i)) + 
     &		log(1.0 + (0.1400/96.062)*(scl)/ks(i)))
C        Calculate concentrations for borate, sulfate, and fluoride
C        Uppstrom (1974)
	 bt(i) = 0.000232 * scl/10.811
C        Morris & Riley (1966)
	 st(i) = 0.14 * scl/96.062
C        Riley (1965)
	 ft(i) = 0.000067 * scl/18.9984
       enddo

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
c       call toc('co2_flux')

c       call tic('co2_flux_iter')

#ifdef  calculate_C_eq
c    calculate DIC equilibrium concentration
       do i=1,imt
        atmpres=slp(i)/1013.25 ! mbar to atm
        co2star_eq(i)= pco2_air*permeg*ff(i)*atmpres 
       enddo
       do i=1,imt
	htotal(i) = drtsafe(dic_iter_C_eq,lower_bound,upper_bound,accuracy)
       enddo
       do i=1,imt
        dic_eq(i)=co2star_eq(i)*
     &            (1.0+k1(i)/htotal(i)+k1(i)*k2(i)/htotal(i)**2)
     &             /rrho02/1000.0  ! in mol/m^3
       enddo
#endif

#ifdef use_fast_drtsafe_routine
         MAXIT=100
         x1 = lower_bound
         x2 = upper_bound
         xacc = accuracy
       do i=1,imt
c        htotal(i) = drtsafe(dic_iter,lower_bound,upper_bound,accuracy)
c        REAL FUNCTION DRTSAFE(FUNCD,X1,X2,XACC)
C        File taken from Numerical Recipes. Modified  R.M.Key 4/94
         uuu= dic_iter(X1,FL,DF)
         uuu= dic_iter(X2,FH,DF)
         IF(FL .LT. 0.0) THEN
           XL=X1
           XH=X2
         ELSE
           XH=X1
           XL=X2
           SWAP=FL
           FL=FH
           FH=SWAP
         END IF
         DRTSAFE=.5*(X1+X2)
         DXOLD=ABS(X2-X1)
         DX=DXOLD
         uuu= dic_iter(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) goto 101
           ELSE
             DXOLD=DX
             DX=F/DF
             TEMP=DRTSAFE
             DRTSAFE=DRTSAFE-DX
             IF(TEMP .EQ. DRTSAFE) goto 101
	   END IF
           IF(ABS(DX) .LT. XACC) goto 101
           uuu= dic_iter(DRTSAFE,F,DF)
           IF(F .LT. 0.0) THEN
             XL=DRTSAFE
             FL=F
           ELSE
             XH=DRTSAFE
             FH=F
           END IF
  100    CONTINUE
  101    continue
c        RETURN
c        END function drtsafe
         htotal(i) = drtsafe
       enddo
#else
       do i=1,imt
	htotal(i) = drtsafe(dic_iter,lower_bound,upper_bound,accuracy)
       enddo
#endif
#ifdef calc_htotal_also
        htotal_out(1:imt)=htotal(1:imt)
#endif
c       call toc('co2_flux_iter')

c       call tic('co2_flux')

       do i=1,imt
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)
	htotal2=htotal(i)*htotal(i)
	co2star=dic(i)*htotal2/(htotal2 + k1(i)*htotal(i) + k1(i)*k2(i)) ! mol/kg
        atmpres=slp(i)/1013.25 ! mbar to atm
        dco2(i)= ( pco2_air*permeg*ff(i)*atmpres - co2star )/permil   ! mol/m^3
        pco2_water(i) = co2star/ff(i)/permeg/atmpres                  ! ppm or uatm

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
c              co2star  = co2star / permil
c             Note: pCO2surf and dpCO2 are calculated in atm above. 
c             Thus convert now to uatm
c              pCO2(i) = pCO2(i) / permeg
C
      enddo
c      call toc('co2_flux')

      end subroutine co2calc





      subroutine co2calc_fast(t,s_in,dic_in,nit,slp,imt,
     &                   pco2_air,dco2,pco2_water
     &)
c-------------------------------------------------------------------------
c   Calculate pco2* from total alkalinity and total CO2 at
c   temperature (t), salinity (s) and "atmpres" atmosphere total pressure. 
c   using a fit of pH as a function of DIC, SALT and TEMP
c   valid range is 1940<DIC<2200 mmol/m^3
c                  -5<T<25 degC and 31.5<S<37.5 PSU
c   accuracy within this range is as good as for other routine
c--------------------------------------------------------------------------
       use timing_module
       implicit none
c      INPUT
       integer :: imt         ! vector length
       real   :: t(imt)        ! temperature deg C
       real   :: s_in(imt)     ! salinity in (PSU-35)/1000
       real   :: dic_in(imt)   ! diss. inorg. carbon  in mmol/m^3
       real   :: nit(imt)      ! nitrate in mmol/m^3 not used anyway
       real   :: slp(imt)      ! Sea level pressure in mbar
       real   :: pco2_air      ! atmospheric mole fraction CO2 in dry air (ppmv) 
c      OUTPUT
       real   :: dco2(imt)     ! difference in CO2 at sea surface and saturation CO2 in mol/m^3
       real   :: pco2_water(imt) ! partial pressure of CO2 at sea surface  (ppmv)

      real, parameter :: rho0 = 1024.0    ! unit conversion factor
      real, parameter :: rrho02 = 1.0/rho0**2
      real, parameter :: permeg=1.e-6, permil = 1.0/rho0 ! unit conversion
c      local work variables
       real :: tk,tk100,tk1002,invtk,dlogtk,s2
       real :: htotal,s,htotal2,co2star,pco2
       real :: dic,ff,k1,k2,atmpres,pH_fit,xx,yy,zz
       real :: xx2,yy2,zz2
       integer :: i

       do i=1,imt
        s = s_in(i)*1000.0+35.0
        dic = dic_in(i)*rrho02
        tk = 273.15 + t(i)
        tk100 = tk/100.0
        tk1002=tk100*tk100
        invtk=1.0/tk
        dlogtk=log(tk)
        s2=s*s
C        f = k0(1-pH2O)*correction term for non-ideality
C        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))
C        k1 = [H][HCO3]/[H2CO3]
C        k2 = [H][CO3]/[HCO3]
C        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))
c    calculate pH from estimated fit
      xx=dic_in(i)-2000.0; yy=s-35.0; zz=t(i);
      xx2=xx*xx; yy2=yy*yy; zz2=zz*zz;
      pH_fit=8.456655e+00+(-1.688880e-03)*xx 
     & +(7.390822e-02)*yy+(-1.739867e-02)*zz 
     & +(-1.864020e-06)*xx2+(-2.034953e-03)*yy2 
     & +(6.480540e-05)*zz2 
     & +(8.483447e-07)*xx*zz+(1.945036e-04)*xx*yy 
     & +(-4.480909e-05)*yy*zz 
     & +(-2.946383e-09)*xx**3+(-1.420486e-03)*yy**3 
     & +(-4.430561e-07)*zz**3 
     & +(-1.515261e-05)*xx*yy2+(5.880059e-07)*yy*xx2  
     & +(4.673103e-08)*xx*zz2+(1.772955e-08)*zz*xx2 
     & +(-1.897002e-06)*yy*zz2+(2.194509e-05)*zz*yy2 
     & +(-1.493190e-06)*xx*yy*zz   
	htotal = 10**(-pH_fit)
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)
	htotal2=htotal*htotal
	co2star=dic*htotal2/(htotal2 + k1*htotal + k1*k2) ! mol/kg
        atmpres=slp(i)/1013.25 ! mbar to atm
        dco2(i)= ( pco2_air*permeg*ff*atmpres - co2star )/permil   ! mol/m^3
        pco2_water(i) = co2star/ff/permeg/atmpres                  ! ppm or uatm
      enddo
      end subroutine co2calc_fast



      real function dic_iter(x,fn,df)
c-------------------------------------------------------------------------
c   This routine expresses TA as a function of DIC, htotal and constants.
c   It also calculates the derivative of this function with respect to 
c   htotal. It is used in the iterative solution for htotal. In the call
c   "x" is the input value for htotal, "fn" is the calculated value for TA
c   and "df" is the value for dTA/dhtotal
c-------------------------------------------------------------------------
       use pco2_module
       implicit none
       real :: x
       real :: fn,df
       real :: k12,k12p,k123p,x2,x3,a,c,a2,da,b,b2,db

	x2=x*x
	x3=x2*x
	k12 = k1(i)*k2(i)
	k12p = k1p(i)*k2p(i)
	k123p = k12p*k3p(i)
	c = 1.0 + st(i)/ks(i)
	a = x3 + k1p(i)*x2 + k12p*x + k123p
	a2=a*a
	da = 3.0*x2 + 2.0*k1p(i)*x + k12p
	b = x2 + k1(i)*x + k12
	b2=b*b
	db = 2.0*x + k1(i)
C
C	fn = hco3+co3+borate+oh+hpo4+2*po4+silicate+hfree+hso4+hf+h3po4-ta
C
	fn = k1(i)*x*dic(i)/b +
     &	     2.0*dic(i)*k12/b +
     &	     bt(i)/(1.0 + x/kb(i)) +
     &	     kw(i)/x +
     &	     pt(i)*k12p*x/a +
     &	     2.0*pt(i)*k123p/a +
     &	     sit/(1.0 + x/ksi(i)) -
     &	     x/c -
     &	     st(i)/(1.0 + ks(i)/x/c) -
     &	     ft(i)/(1.0 + kf(i)/x) -
     &	     pt(i)*x3/a -
     &	     ta(i)
C
C	df = dfn/dx
C
	df = ((k1(i)*dic(i)*b) - k1(i)*x*dic(i)*db)/b2 -
     &	     2.0*dic(i)*k12*db/b2 -
     &	     bt(i)/kb(i)/(1.0+x/kb(i))**2 -
     &	     kw(i)/x2 +
     &	     (pt(i)*k12p*(a - x*da))/a2 
     &	     - 2.0*pt(i)*k123p*da/a2 
     &	     - sit/ksi(i)/(1.0+x/ksi(i))**2 - 1.0/c 
     &	     + st(i)*(1.0 + ks(i)/x/c)**(-2)*(ks(i)/c/x2) 
     &	     + ft(i)*(1.0 + kf(i)/x)**(-2)*kf(i)/x2 
     &	     - pt(i)*x2*(3.0*a-x*da)/a2

      dic_iter = 0.0

      end function dic_iter


#ifdef  calculate_C_eq

      real function dic_iter_C_eq(x,fn,df)
c-------------------------------------------------------------------------
c   This routine expresses TA as a function of CO2star, htotal and constants.
c   It also calculates the derivative of this function with respect to 
c   htotal. It is used in the iterative solution for htotal. In the call
c   "x" is the input value for htotal, "fn" is the calculated value for TA
c   and "df" is the value for dTA/dhtotal
c-------------------------------------------------------------------------
       use pco2_module
       implicit none
       real :: x
       real :: fn,df
       real :: k12,k12p,k123p,x2,x3,a,c,a2,da,b,b2,db

	x2=x*x
	x3=x2*x
	k12 = k1(i)*k2(i)
	k12p = k1p(i)*k2p(i)
	k123p = k12p*k3p(i)
	c = 1.0 + st(i)/ks(i)
	a = x3 + k1p(i)*x2 + k12p*x + k123p
	a2=a*a
	da = 3.0*x2 + 2.0*k1p(i)*x + k12p
	b = x2 + k1(i)*x + k12
	b2=b*b
	db = 2.0*x + k1(i)
C
C	fn = hco3+co3+borate+oh+hpo4+2*po4+silicate+hfree+hso4+hf+h3po4-ta
C
	fn = 
     &       k1(i)/x*co2star_eq(i) +
     &       2*k12/x2*co2star_eq(i) +
     &	     bt(i)/(1.0 + x/kb(i)) +
     &	     kw(i)/x +
     &	     pt(i)*k12p*x/a +
     &	     2.0*pt(i)*k123p/a +
     &	     sit/(1.0 + x/ksi(i)) -
     &	     x/c -
     &	     st(i)/(1.0 + ks(i)/x/c) -
     &	     ft(i)/(1.0 + kf(i)/x) -
     &	     pt(i)*x3/a -
     &	     ta(i)
C
C	df = dfn/dx
C
	df = 
     &       -k1(i)/x2*co2star_eq(i) 
     &       -4*k12/x3*co2star_eq(i) +
     &	     bt(i)/kb(i)/(1.0+x/kb(i))**2 -
     &	     kw(i)/x2 +
     &	     (pt(i)*k12p*(a - x*da))/a2 
     &	     - 2.0*pt(i)*k123p*da/a2 
     &	     - sit/ksi(i)/(1.0+x/ksi(i))**2 - 1.0/c 
     &	     + st(i)*(1.0 + ks(i)/x/c)**(-2)*(ks(i)/c/x2) 
     &	     + ft(i)*(1.0 + kf(i)/x)**(-2)*kf(i)/x2 
     &	     - pt(i)*x2*(3.0*a-x*da)/a2

      dic_iter_C_eq = 0.0

      end function dic_iter_C_eq

#endif


      REAL FUNCTION DRTSAFE(FUNCD,X1,X2,XACC)
C     File taken from Numerical Recipes. Modified  R.M.Key 4/94
      MAXIT=100
      uuu= FUNCD(X1,FL,DF)
      uuu= 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
      END IF
      DRTSAFE=.5*(X1+X2)
      DXOLD=ABS(X2-X1)
      DX=DXOLD
      uuu= 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
	END IF
        IF(ABS(DX) .LT. XACC)RETURN
        uuu= FUNCD(DRTSAFE,F,DF)
        IF(F .LT. 0.0) THEN
          XL=DRTSAFE
          FL=F
        ELSE
          XH=DRTSAFE
          FH=F
        END IF
  100  CONTINUE
      RETURN
      END function drtsafe




#ifdef drive_test_co2

      program test
       use timing_module
      implicit none
      real, allocatable :: t(:),s(:),dic(:),pco2_air(:),nit(:)
      integer :: imt = 1,n1,n2,n3,n4,i,i1,i2,i3,i4,ii,n5,i5

      real :: trange(3) = (/-5.,25.,1.0/)
      real :: srange(3) = (/32.,38.,1.0/) ! start, end and spacing in TS space
      real :: nitrange(3) = (/980.,1030.,10.0/) ! start, end and spacing in TS space
      real :: dicrange(3) = (/1950.,2200.0,25.0/) 
      real :: pco2range(3) = (/278.0,388.0,10.0/) 
      real  :: pco2_water,dco2,dic_eq,htotal

      character(len=80) :: name,unit

#include "netcdf.inc"
      integer :: ncid,iret,tdim,sdim,tid,sid,pco2id,dicid
      integer :: dicdim,dco2id,pco2dim,pco2_waterid
      integer :: dic_eqid,htotid,pco2_water2id,dco22id

      n1=(trange(2)-trange(1))/trange(3)
      n2=(srange(2)-srange(1))/srange(3)
      n3=(dicrange(2)-dicrange(1))/dicrange(3)
      n4=(pco2range(2)-pco2range(1))/pco2range(3)
      n5=(nitrange(2)-nitrange(1))/nitrange(3)
      print*,' n1=',n1,' n2=',n2,' n3=',n3,' n4=',n4,' n5=',n5
      allocate( t(n1), s(n2), dic(n3), pco2_air(n4),nit(n5) )

      t(1)=trange(1)
      do i=2,n1; t(i)=t(i-1)+trange(3); enddo
      s(1)=srange(1)
      do i=2,n2; s(i)=s(i-1)+srange(3); enddo
      dic(1)=dicrange(1)
      do i=2,n3; dic(i)=dic(i-1)+dicrange(3); enddo
      pco2_air(1)=pco2range(1)
      do i=2,n4; pco2_air(i)=pco2_air(i-1)+pco2range(3); enddo
      nit(1)=nitrange(1)
      do i=2,n5; nit(i)=nit(i-1)+nitrange(3); enddo

      ncid   = nccre ('co2_test.cdf', NCCLOB, iret)
      iret   = nf_set_fill(ncid, NF_NOFILL, iret)
      tdim   = ncddef(ncid, 'x', n1, iret)
      sdim   = ncddef(ncid, 'y', n2, iret)
      dicdim = ncddef(ncid, 'dic', n3, iret)

      tid  = ncvdef (ncid,'x',NCFLOAT,1,tdim,iret)
      name = 'Potential temparature'
      call ncaptc(ncid,tid,'long_name',NCCHAR,len_trim(name),name,iret) 
      name = 'deg C'
      call ncaptc(ncid,tid,'units',NCCHAR,len_trim(name),name,iret) 

      sid  = ncvdef (ncid,'y',NCFLOAT,1,sdim,iret)
      name = 'Salinity'
      call ncaptc(ncid,sid,'long_name',NCCHAR,len_trim(name),name,iret) 
      name = 'psu'
      call ncaptc(ncid,sid,'units',NCCHAR,len_trim(name),name,iret) 

      dicid  = ncvdef (ncid,'dic',NCFLOAT,1,dicdim,iret)
      name = 'Dissolved inorganic carbon'
      call ncaptc(ncid,dicid,'long_name',
     &            NCCHAR,len_trim(name),name,iret) 
      name = 'mmol/m**3'
      call ncaptc(ncid,dicid,'units',NCCHAR,len_trim(name),name,iret) 

      pco2dim = ncddef(ncid, 'pco2', n4, iret)
c      pco2dim = ncddef(ncid, 'nit', n5, iret)
      pco2id  = ncvdef (ncid,'pco2',NCFLOAT,1,pco2dim,iret)
c      pco2id  = ncvdef (ncid,'nit',NCFLOAT,1,pco2dim,iret)
      name = 'CO2 partial pressure'
c      name = 'SLP'
      call ncaptc(ncid,pco2id,'long_name',
     &            NCCHAR,len_trim(name),name,iret) 
      name = 'ppmv'
c      name = 'mbar'
      call ncaptc(ncid,pco2id,'units',NCCHAR,len_trim(name),name,iret) 

      pco2_waterid  = ncvdef (ncid,'pCO2_water_fit',NCFLOAT,4,
     &                   (/tdim,sdim,dicdim,pco2dim/),iret)
      name = 'Oceanic pCO2 at sea surface'
      call ncaptc(ncid,pco2_waterid,'long_name',NCCHAR,
     &            len_trim(name),name,iret) 
      name = 'ppmv'
      call ncaptc(ncid,pco2_waterid,'units',
     &  NCCHAR,len_trim(name),name,iret) 


      pco2_water2id  = ncvdef (ncid,'pCO2_water',NCFLOAT,4,
     &                   (/tdim,sdim,dicdim,pco2dim/),iret)
      name = 'Oceanic pCO2 at sea surface from fit'
      call ncaptc(ncid,pco2_water2id,'long_name',NCCHAR,
     &            len_trim(name),name,iret) 
      name = 'ppmv'
      call ncaptc(ncid,pco2_water2id,'units',
     &  NCCHAR,len_trim(name),name,iret) 


      dco2id  = ncvdef (ncid,'dCO2',NCFLOAT,4,
     &                   (/tdim,sdim,dicdim,pco2dim/),iret)
      name = 'Air-sea diff of CO2 concentration'
      call ncaptc(ncid,dco2id,'long_name',NCCHAR,
     &            len_trim(name),name,iret) 
      name = 'mol/m^3'
      call ncaptc(ncid,dco2id,'units',NCCHAR,len_trim(name),name,iret) 

      dco22id  = ncvdef (ncid,'dCO2_fit',NCFLOAT,4,
     &                   (/tdim,sdim,dicdim,pco2dim/),iret)
      name = 'Air-sea diff of CO2 concentration from fit'
      call ncaptc(ncid,dco22id,'long_name',NCCHAR,
     &            len_trim(name),name,iret) 
      name = 'mol/m^3'
      call ncaptc(ncid,dco22id,'units',NCCHAR,len_trim(name),name,iret) 


#ifdef  calculate_C_eq
      dic_eqid  = ncvdef (ncid,'DIC_eq',NCFLOAT,4,
     &                   (/tdim,sdim,dicdim,pco2dim/),iret)
      name = 'DIC equilibrium concentration'
      call ncaptc(ncid,dic_eqid,'long_name',NCCHAR,
     &            len_trim(name),name,iret) 
      name = 'mol/m^3'
      call ncaptc(ncid,dic_eqid,'units',NCCHAR,len_trim(name),name,iret) 
#endif
#ifdef calc_htotal_also
      htotid  = ncvdef (ncid,'htotal',NCFLOAT,4,
     &                   (/tdim,sdim,dicdim,pco2dim/),iret)
      name = 'H+ concentration'
      call ncaptc(ncid,htotid,'long_name',NCCHAR,
     &            len_trim(name),name,iret) 
      name = 'mol/m^3'
      call ncaptc(ncid,htotid,'units',NCCHAR,len_trim(name),name,iret) 
#endif

      call ncendf(ncid, iret)

      call ncvpt(ncid, tid, 1, n1,t,iret)
      call ncvpt(ncid, sid, 1, n2,s,iret)
      call ncvpt(ncid, dicid, 1, n3,dic,iret)
      call ncvpt(ncid, pco2id, 1, n4,pco2_air,iret)
c      call ncvpt(ncid, pco2id, 1, n5,nit,iret)


      do i1=1,n1
      do i2=1,n2
      do i3=1,n3
      do i4=1,n4
c      do i4=1,n5
       call tic('old')
        CALL co2calc
     &  (t(i1),(s(i2)-35)/1000.0,dic(i3),2.0,1013.25,1,
     &                 pco2_air(i4),dco2,pco2_water
c     &  (t(i1),(s(i2)-35)/1000.0,dic(i3),2.0,nit(i4),1,
c     &                 pco2_air(1),dco2,pco2_water
#ifdef  calculate_C_eq
     &              , dic_eq
#endif
#ifdef calc_htotal_also
     &              , htotal
#endif
     &)
       call toc('old')


        call ncvpt(ncid, pco2_waterid, (/i1,i2,i3,i4/), (/1,1,1,1/),
     &              pCO2_water,iret)
        call ncvpt(ncid, dco2id, (/i1,i2,i3,i4/), (/1,1,1,1/),dCO2,iret)
#ifdef  calculate_C_eq
        call ncvpt(ncid,dic_eqid, (/i1,i2,i3,i4/), 
     &                         (/1,1,1,1/),dic_eq,iret)
#endif
#ifdef calc_htotal_also
        call ncvpt(ncid,htotid, (/i1,i2,i3,i4/), 
     &                         (/1,1,1,1/),htotal,iret)
#endif


       call tic('new')
        CALL co2calc_fast
     &  (t(i1),(s(i2)-35)/1000.0,dic(i3),2.0,1013.25,1,
     &                 pco2_air(i4),dco2,pco2_water
c     &  (t(i1),(s(i2)-35)/1000.0,dic(i3),2.0,nit(i4),1,
c     &                 pco2_air(1),dco2,pco2_water
#ifdef  calculate_C_eq
     &              , dic_eq
#endif
#ifdef calc_htotal_also
     &              , htotal
#endif
     &)
       call toc('new')

        call ncvpt(ncid, pco2_water2id, (/i1,i2,i3,i4/), (/1,1,1,1/),
     &              pCO2_water,iret)
        call ncvpt(ncid, dco22id,(/i1,i2,i3,i4/),(/1,1,1,1/),dCO2,iret)

      enddo
      enddo
      enddo
      enddo

      call ncclos (ncid, iret)


      print*,' old time summary    = ',
     &   timing_secs('old') ,' s'
      print*,' new time summary    = ',
     &   timing_secs('new') ,' s'

#ifdef notdef
      t        = 10.0
      s        = (35.0-35.0)/1000.0
      Nit      = 5.0
      dic      = 1980.0
      pco2_air = 270.0

      print*,' Input : '
      print*,' Temperature = ',t,' deg C'
      print*,' Salinity    = ',s*1000.0+35,' psu'
      print*,' DIC         = ',dic,' mmol/m^3'
      print*,' Nitrate     = ',nit,' mmol/m^3'
      print*,' Pco2 air    = ',pco2_air,' ppmv'
      CALL co2calc(t,s,dic,nit,1013.25,1,pco2_air,dco2,pco2_water)
      print*,' Output : '
      print*,' Pco2 water  = ',pco2_water,' ppmv'
      print*,' Air-sea diff. of CO2 conc. = ',dco2,' mol/m^3'
#endif 

      end program test


      subroutine halt_stop(string)
      character*(*) :: string
      print*,string
      stop
      end 
      subroutine barrier()
      end 


#endif

