C**********************************************************************C
C                                                                      C
C   soot.f (02/15/02)      Ken Revzan and Michael Frenklach            C
C                                                                      C
C   Soot model includes:                                               C
C      - 6 concentration moments (5 size moments)                      C
C      - nucleation: A4 + A4 -> dimer (i.e., no PAH growth past A4)    C
C      - coagulation for all pressure regimes (harmonic mean)          C
C      - no particle aggregation                                       C
C      - PAH (i.e., A4) condensation on soot particle surface          C
C      - surface growth: basic HACA (C2H2)                             C
C      - surface oxidation by O2 and OH                                C
C                                                                      C
C   initSoot must be called before first call to any other routine     C
C                                                                      C
C**********************************************************************C

	subroutine initSoot(n_moments, iPAH, iC2H2, iCO, iH, iH2,
     &                      iH2O, iO2, iOH, iM0, freeMolOnly, M0_cutoff,
     &	                    PAH_carbon, PAH_diam)

C	Initialization of variables
C	Must be called before any other routine in this module

C	Input:
C	  n_moments:	number of moments
C	  iPAH:		index of inception PAH species in species array,
C			as input to function soot
C	  iC2H2:	index of C2H2
C	  iCO:		index of CO
C	  iH:		index of H
C	  iH2:		index of H2
C	  iH2O:		index of H2O
C	  iO2:		index of O2
C	  iOH:		index of OH
C	  iM0:		index of M0 (zeroth moment)
C	  freeMolOnly:	true if coagulation rates are to be calculated
C			for free molecular regime, false if for
C			transition regime
C	  M0_cutoff:	Only nucleation rates are calculated for values
C			of M0 below M0_cutoff
C	  PAH_carbon:	Number of carbon atoms in inception species
C	  PAH_diam:	Diameter of inception species (cm)

	implicit none

	include 'soot.h'

	integer n_moments, iPAH, iC2H2, iCO, iH, iH2, iH2O,
     &	                                     iO2, iOH, iM0
	logical freeMolOnly
	doubleprecision M0_cutoff, PAH_diam
	integer PAH_carbon

	integer n, m, r, l
	double precision CD1_PAH, CH2, CHS, CS2

C	Routine set for 6 moments only, for now

	if(n_moments .ne. 6) then
		write(6,*) 'Number of moments must be 6'
		stop
	endif

C	Common block variables

	nMoments = n_moments
	jPAH = iPAH
	jC2H2 = iC2H2
	jCO = iCO
	jH = iH
	jH2 = iH2
	jH2O = iH2O
	jO2 = iO2
	jOH = iOH
	jM0 = iM0
	fmOnly = freeMolOnly
	cutoff = M0_cutoff
	CPerPAH = PAH_carbon
	diamPAH = PAH_diam

C	Definition of basic constants - not included in header file
C	because of compiler complaints

C	particle diameter = CD1 * (number of C atoms)**(1/3)

	CD1 = (6 * C_MASS / PI / RHO_SOOT)**ONE_THIRD

C	cbNucl is PAH+PAH constant (for nucleation)

	cbNucl = 2.2 * SQRT(4 * PI * BOLTZMANN / C_MASS / CPerPAH)
     &           * diamPAH**2 * AVOGADRO**2              

	CPerDimer = 2 * CPerPAH

C	CB0 is beta for free-molecular coagulation

	CB0 = 2.2 * (3 * C_MASS / 4 / PI / RHO_SOOT)**ONE_SIXTH
     &          * SQRT(6 * BOLTZMANN / RHO_SOOT)    

C	CA is constant 'a' for continuum coagulation
C	(needs to be multiplied by mean free path)

	CA  = 2.514 / CD1

C	CBOH is OH+SOOT oxidation constant  

	CBOH = CD1**2 * SQRT(PI * BOLTZMANN / 2 / OH_MASS) * AVOGADRO

C	A constant needed for surface growth

	SURFC = CHI * PI * CD1**2

C	PAH diameter = CD1_PAH * (number of C atoms)**(1/2)  

	CD1_PAH = D1_PAH * SQRT(TWO_THIRDS)

C	A4+SOOT condensation constant

	CBCOND = 2.2 * SQRT(PI * BOLTZMANN / 2 / C_MASS) * AVOGADRO

C	Constants used in continuum coagulation
C	Results from (d_pah + d_soot)^2

	CH2 = CBCOND * CD1_PAH**2               ! const * C_h^2
	CHS = CBCOND * 2 * CD1 * CD1_PAH        ! const * 2 * C_h * C_s
	CS2 = CBCOND * CD1**2                   ! const * C_s^2

C	This for PAH surface condensation

	do n=1,nMoments-1
	  powPAH(1,n) = CH2 * CPerPAH**(dfloat(2*n+1)/2.0D0)
	  powPAH(2,n) = CHS * CPerPAH**n
	  powPAH(3,n) = CS2 * CPerPAH**(dfloat(2*n-1)/2.0D0)
	enddo

C	Vector for use in Lagrangian

	do n=1,nMoments
	  x(n) = dfloat(n-1)
	enddo
	do n=1,nMoments
	  do m=1,n
	    prime(n,m) = 1.0
	    do l=1,n
	      if(l .ne. m) prime(n,m) = prime(n,m) * (x(m) - x(l))
	    enddo
	  enddo
	enddo

C	Binomial coefficients

	do r=0,nMoments
	  binom(0,r) = 1
	  binom(r,r) = 1
	  do l=1,r
	    binom(l,r) = binom(l-1,r-1) + binom(l,r-1)
	  enddo
	enddo

	first = .true.

	return

	end

C**********************************************************************C
C
	integer function soot(Y, T, P, rate, nucRate, coagRate,
     &	                      surfRate, surfPAH, surfC2H2, surfO2,
     &	                      surfOH, prodPAH, prodC2H2, prodCO, prodH,
     &	                      prodH2, prodH2O, prodO2, prodOH)
C                                                                      
C	Function SOOT calculates the rates of the first 6
C	concentration moments of the soot particle size distribution
C	function and surface production rates for gaseous species
C	participating in soot nucleation, growth, and oxidation.
C                                                                      
C	Input:
C	  Y:		Array of chemical species and soot moments;
C			chemical species concentrations are in mol/cm^3;
C			soot moments are in number/cm^3
C	  T:		Temperature (K)
C	  P:		Pressure (atm)
C	Output:
C	  rate:		Array of moment rates
C	  nucRate:	Array of nucleation rates
C	  coagRate:	Array of coagulation rates
C	  surfRate:	Array of total surface growth/oxidation/
C			condensation rates
C	  surfPAH:	Array of condensation rates
C	  surfC2H2:	Array of C2H2 surface growth rates
C	  surfO2:	Array of O2 oxidation rates
C	  surfOH:	Array of OH oxidation rates
C			All rates in cm^-3
C	  prodPAH:	Surface production rate for inception species
C	  prodC2H2:	Surface production rate for C2H2
C	  prodCO:	Surface production rate for CO
C	  prodH:	Surface production rate for H
C	  prodH2:	Surface production rate for H2
C	  prodH2O:	Surface production rate for H2O
C	  prodO2:	Surface production rate for O2
C	  prodOH:	Surface production rate for OH
C			All production rates in mol cm^-3
C			Negative production rates indicate consumption
C	Return values:
C	   0:		Normal return
C	  -1:		Dispersion <= 1
C	  -2:		Moment M(r) <= M(r-1)
C----------------------------------------------------------------------
C 
	implicit none

	include 'soot.h'

	double precision Y(*), T, P, rate(0:*), nucRate(0:*),
     &	                 coagRate(0:*), surfRate(0:*), surfPAH(0:*),
     &	                 surfC2H2(0:*), surfO2(0:*), surfOH(0:*),
     &	                 prodPAH, prodC2H2, prodCO, prodH, prodH2,
     &	                 prodH2O, prodO2, prodOH

	integer n, highMoment
	doubleprecision srT, f(MAX_MOMENTS), prodRate, fitpos, fitneg,
     &	                dispersion, a, b, rsq

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	highMoment = nMoments - 1

C	Initialization of output vars

	do n=0,highMoment
	  nucRate(n) = 0.0
	  coagRate(n) = 0.0
	  surfRate(n) = 0.0
	  surfPAH(n) = 0.0
	  surfC2H2(n) = 0.0
	  surfO2(n) = 0.0
	  surfOH(n) = 0.0
	enddo
	prodPAH = 0
	prodC2H2 = 0
	prodCO = 0
	prodH = 0
	prodH2 = 0
	prodH2O = 0
	prodO2 = 0
	prodOH = 0

	srT = sqrt(T)

C	Nucleation; don't allow negative values

	if(Y(jPAH) .gt. 0.0) then
	  nucRate(0) = cbNucl * srT * Y(jPAH) * Y(jPAH)
	else
	  nucRate(0) = 0.0
	endif

C	Consumption rate of PAH due to nucleation

	prodPAH = -2 * nucRate(0) / AVOGADRO

	do n=1,highMoment
	  nucRate(n) = CPerDimer * nucRate(n-1)
	enddo

C	f is vector of size moments for Lagrange interp
C	Index of highest fractional moment is nMoments+1/6

	f(1) = 1.0
	do n=1,highMoment
	  f(n+1) = Y(jM0+n) / Y(jM0)
	enddo
	do n=LO_FRAC,-1
	  mu(n) = fitneg(dfloat(n)/6.0, 3, f)
	enddo
	do n=0,6*highMoment+1
	  if(mod(n,6) .ne. 0) then
	    mu(n) = fitpos(dfloat(n)/6.0, nMoments, f)
	  else
	    mu(n) = f(n/6+1)
	  endif
	enddo
	do n=6*highMoment+2,HI_FRAC
	  mu(n) = 0.0
	enddo

C	User-chosen cutoff to allow buildup without coagulation or
C	surface growth/oxidation. Place after calculation of mu to
C	avoid problems with Fv, D, etc.

	if(Y(jM0) .lt. cutoff) then
	  do n=0,highMoment
	    rate(n) = nucRate(n)
	  enddo
	  soot = 0
	  return
	endif

	call coagulation(P, T, Y(jM0), coagRate)

	call condensation(highMoment, srT, Y(jM0), Y(jPAH),
     &	                  surfPAH, prodRate)
	prodPAH = prodPAH + prodRate

	call surface(highMoment, T, Y(jC2H2), Y(jH), Y(jH2), Y(jH2O),
     &	             Y(jO2), Y(jOH), Y(jM0), surfC2H2, surfO2, surfOH,
     &	             surfRate, prodC2H2, prodCO, prodH, prodH2,
     &	             prodH2O, prodO2, prodOH)
	surfRate(0) = 0.0	! No surface growth contrib to M0

C       Oxidation regime - base surface rates on M1, ...

	if(surfRate(1) .lt. 0.0) then
	  do n=1,highMoment
	    f(n) = log(mu(6*n))
	  enddo
	  call linear(highMoment, f, a, b, rsq)
	  do n=4,6*highMoment-2,6
	    mu(n) = exp(a + b * n / 6.0)
	  enddo
          call surface(highMoment, T, Y(jC2H2), Y(jH), Y(jH2),
     &	               Y(jH2O), Y(jO2), Y(jOH), Y(jM0), surfC2H2,
     &	               surfO2, surfOH, surfRate, prodC2H2, prodCO,
     &	               prodH, prodH2, prodH2O, prodO2, prodOH)
	  do n=1,highMoment
	    surfRate(n) =
     &                  surfRate(n) * Y(jM0+n) / exp(a + b * n) / Y(jM0)
	  enddo
	  surfRate(0) = 0.0
	endif

	do n=0,highMoment
	  surfRate(n) = surfRate(n) + surfPAH(n)
	enddo

C	Rates of the moments

	do n=0,highMoment
	  rate(n) = nucRate(n) + coagRate(n) + surfRate(n)
	enddo

C	Error on non-physical dispersion

	if(dispersion() .lt. 1.0) then
	  soot = -1
	  return
	endif

C	Error on non-physical moments

	do n=1,highMoment
	  if(Y(jM0+n) .le. Y(jM0)) then
	    soot = -2
	    return
	  endif
	enddo

	soot = 0

	return

	end

C**********************************************************************C

	doubleprecision function avgSootDiam()

	implicit none

	include 'soot.h'

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	avgSootDiam = CD1 * mu(2)

	return

	end

C**********************************************************************C

	doubleprecision function betaC1(i, j)

C	First part of collision rate used in continuum coagulation

	implicit none

	include 'soot.h'

	integer i, j

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	integer i6, j6

	i6 = 6 * i
	j6 = 6 * j

	betaC1 =     mu(-2+i6) * mu(2+j6) +
     &	         2 * mu(i6)    * mu(j6) +
     &	             mu(2+i6)  * mu(-2+j6)

	return

	end

C**********************************************************************C

	doubleprecision function betaC2(i, j)

C	Second (viscosity) part of collision rate used in continuum
C	coagulation

	implicit none

	include 'soot.h'

	integer i, j

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	integer i6, j6

	i6 = 6 * i
	j6 = 6 * j

	betaC2 = mu(-4+i6) * mu(2+j6) +
     &	         mu(-2+i6) * mu(j6) +
     &	         mu(i6)    * mu(-2+j6) +
     &	         mu(2+i6)  * mu(-4+j6)

	return

	end

C**********************************************************************C

	subroutine coagulation(P, T, M0, rate)

C	Coagulation rates calculated using Method II of M. Frenklach
C	and S.Harris, J. Colloid Interface Sci., v.118, n.1, 1987,
C	pp. 252-261. See also M. Frenklach and H. Wang, "Soot Formation
C	in Combustion: Mechanisms and Models" (H. Bockhorn, Ed.),
C	Springer-Verlag, Heidelberg, 1994, p. 165.

C	Input:
C	  P:		pressure (atm)
C	  T:		temperature (K)
C	  fmOnly:	if true, free molecular regime, else transition
C	Output:
C	  rate:		array of coagulation rates

	implicit none

	include 'soot.h'

	doubleprecision P, T, M0, rate(0:*)

	doubleprecision CKT, CAT, crk, crk2, crk3, crk4A, crk4B, crk5A,
     &	                crk5B, crkc, crk2c, crk3c, crk4Ac, crk4Bc,
     &                  crk5Ac, crk5Bc, betaC1, betaC2, M02
	doubleprecision viscosity, freePath
	doubleprecision CB0_sT, f(4), gridFun, fitpos
	integer n

C	Free molecular regime

	CB0_sT = CB0 * sqrt(T)

	do n = 1,4
	  f(n) = gridFun(n-1, 0, 0)
	enddo
	crk = CB0_sT * fitpos(0.5D0, 4, f)
	do n = 1,4
	  f(n) = gridFun(n-1, 1, 1)
	enddo
	crk2 = CB0_sT * fitpos(0.5D0, 4, f)
	do n = 1,4
	  f(n) = gridFun(n-1, 1, 2)
	enddo
	crk3 = CB0_sT * fitpos(0.5D0, 4, f)
	do n = 1,3
	  f(n) = gridFun(n-1, 1, 3)
	enddo
	crk4a = CB0_sT * fitpos(0.5D0, 3, f)
	do n = 1,4
	  f(n) = gridFun(n-1, 2, 2)
	enddo
	crk4b = CB0_sT * fitpos(0.5D0, 4, f)
	do n = 1,2
	  f(n) = gridFun(n-1, 1, 4)
	enddo
	crk5a = CB0_sT * fitpos(0.5D0, 2, f)
	do n = 1,3
	  f(n) = gridFun(n-1, 2, 3)
	enddo
	crk5b = CB0_sT * fitpos(0.5D0, 3, f)

	if(.not.fmOnly) then

C	 Continuum regime

	  CKT = CK * T / viscosity(T)
	  CAT = CA * freePath(P,T)

	  crkc   = CKT * (betaC1(0, 0) + CAT * betaC2(0, 0))
	  crk2c  = CKT * (betaC1(1, 1) + CAT * betaC2(1, 1))
	  crk3c  = CKT * (betaC1(1, 2) + CAT * betaC2(1, 2))
	  crk4Ac = CKT * (betaC1(1, 3) + CAT * betaC2(1, 3))
	  crk4Bc = CKT * (betaC1(2, 2) + CAT * betaC2(2, 2))
	  crk5Ac = CKT * (betaC1(1, 4) + CAT * betaC2(1, 4))
	  crk5Bc = CKT * (betaC1(2, 3) + CAT * betaC2(2, 3))

C	 Harmonic mean of free-molecular and continuum betas

	  crk   = crk   * crkc   / (crk   + crkc)
	  crk2  = crk2  * crk2c  / (crk2  + crk2c)
	  crk3  = crk3  * crk3c  / (crk3  + crk3c)
	  crk4A = crk4A * crk4Ac / (crk4A + crk4Ac)
	  crk4B = crk4B * crk4Bc / (crk4B + crk4Bc)
	  crk5A = crk5A * crk5Ac / (crk5A + crk5Ac)
	  crk5B = crk5B * crk5Bc / (crk5B + crk5Bc)

        endif

	M02 = M0 * M0
	rate(0) = -0.5 * crk * M02
	rate(1) = 0.
	rate(2) = crk2 * M02
	rate(3) = 3 * crk3 * M02
	rate(4) = (4 * crk4A +  3 * crk4B) * M02
	rate(5) = (5 * crk5A + 10 * crk5B) * M02

	return

	end

C**********************************************************************C

	subroutine condensation(highMoment, srT, M0, PAH,
     &	                        rate, prodRate)

C	Calculates the moment rates due to condensation of the PAH
C	inception species on the surface of soot particles and the
C	production rate of PAH

C	Input:
C	  highMoment:		highest moment to be calculated
C	  srT:			square root of temperature
C	  M0:			0th moment
C	  PAH:			PAH concentration (mol/cm^3)
C	Output:
C	  rate:			array of moment rates
C	  prodRate:		PAH surface production rate

	implicit none

	include 'soot.h'

	integer highMoment
	doubleprecision srT, M0, PAH, rate(0:highMoment), prodRate

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	doubleprecision coeff
	integer r, l, l6

	coeff = srT * PAH * M0

	rate(0) = 0.0
	do r=1,highMoment
          rate(r) = 0.0
	  do l=0,r-1
	    l6 = 6 * l
	    rate(r) = rate(r) + binom(l,r) * (
     &                          powPAH(1,r-l) * mu(l6) +
     &                          powPAH(2,r-l) * mu(l6+2) +
     &                          powPAH(3,r-l) * mu(l6+4))
	  enddo
          rate(r) = coeff * rate(r)
	enddo

C	Consumption rate of PAH due to condensation

	prodRate = - rate(1) / CPerPAH / AVOGADRO

	return

	end

C**********************************************************************C

	doubleprecision function dispersion()

	implicit none

	include 'soot.h'

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	dispersion = mu(12) / mu(6) / mu(6)

	return

	end

C**********************************************************************C

	doubleprecision function fitneg(point, numTerms, y)

C	Perform Lagrange interpolation at `point', using values in y
C	x is in common

C	This routine kept separate from fitpos to accommodate
C	future development

	implicit none

	include 'soot.h'

	doubleprecision point
	integer numTerms
	doubleprecision y(numTerms)

	doubleprecision a, prod
	integer n

	prod = 1.0
	do n=1,numTerms
	  prod = prod * (point - x(n))
	enddo

	fitneg = 1.0
	do n=1,numTerms
	  a = prod / ((point - x(n)) * prime(numTerms,n))
	  fitneg = fitneg * (y(n) ** a)
	enddo

	return

	end

C**********************************************************************C

	doubleprecision function fitpos(point, numTerms, y)

C	Perform Lagrange interpolation at `point', using values in y
C	x is in common

	implicit none

	include 'soot.h'

	doubleprecision point
	integer numTerms
	doubleprecision y(numTerms)

	doubleprecision a, prod
	integer n

	prod = 1.0
	do n=1,numTerms
	  prod = prod * (point - x(n))
	enddo

	fitpos = 1.0
	do n=1,numTerms
	  a = prod / ((point - x(n)) * prime(numTerms,n))
	  fitpos = fitpos * (y(n) ** a)
	enddo

	return

	end

C**********************************************************************C

	doubleprecision function freePath(P, T)

C	Mean free path for air in cm (from AK, myf 1997)
C	P in atm;   T in K

	implicit none

	double precision P,T

	freePath = 2.3701D-03 * T / P / 1.103D+05

	return

	end

C**********************************************************************C

	doubleprecision function gridFun(k, n, m)

C	This function used in calculation of free-molecular collision
C	rates, which, for any (n,m), are determined by Lagrange
C	interpolation using two or more gridFun(k,n,m)

	implicit none

	include 'soot.h'

	integer k, n, m

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	integer i, j, l 

	gridFun = 0.0
	do l=0,k
	  i = 6 * (k - l + n)
	  j = 6 * (l + m)
	  gridFun = gridFun + binom(l,k) * (
     &                      mu(i+1) * mu(j-3) +
     &                  2 * mu(i-1) * mu(j-1) +
     &	                    mu(i-3) * mu(j+1) )
	enddo

	return

	end

C**********************************************************************C

        subroutine linear(npts, y, a, b, rsq)

        implicit none

        integer npts
        doubleprecision y(npts), a, b, rsq

        integer n
        doubleprecision sum, x1, x2, y1, y2, D,
     &                  ymean, ypred, ess, yss

        sum = 1.0 * npts
        x1 = 0.0
        x2 = 0.0
        y1 = 0.0
        y2 = 0.0
        do n=1,npts
          x1 = x1 + float(n)
          x2 = x2 + float(n) * float(n)
          y1 = y1 + y(n)
          y2 = y2 + y(n) * float(n)
        enddo
        D = sum * x2 - x1 * x1
        a = (y1 * x2 - x1 * y2) / D
        b = (sum * y2 - y1 * x1) / D
        ymean = y1 / npts
        ess = 0.0
        yss = 0.0
        do n=1,npts
          ypred = a + b * n
          ess = ess + (ymean - ypred) * (ymean - ypred)
          yss = yss + (y(n) - ymean) * (y(n) - ymean)
        enddo
        rsq = ess / yss

        return

        end

C**********************************************************************C

	doubleprecision function sootSurfaceArea(M0)

	implicit none

	include 'soot.h'

	doubleprecision M0

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	sootSurfaceArea = PI * CD1 * CD1 * mu(4) * M0

	return

	end

C**********************************************************************C

	doubleprecision function sootVolumeFraction(M0)

	implicit none

	include 'soot.h'

	doubleprecision M0

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	sootVolumeFraction = PI * CD1 * CD1 * CD1 * mu(6) * M0 / 6

	return

	end

C**********************************************************************C

	subroutine sums(highMoment, massAdded, coef, sum)
C
C       Calculates sums for the surface growth terms
C     
	implicit none

	include 'soot.h'

	integer highMoment, massAdded
	doubleprecision coef, sum(*)

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	integer r, l

	do r=1,highMoment
	  sum(r) = 0.d0
	  do l=0,r-1
	    sum(r) = sum(r) + binom(l,r) * massAdded**(r-l) * mu(6*l+4)
	  enddo
	  sum(r) = coef * sum(r)
	enddo

	return

	end

C**********************************************************************C

	subroutine surface(highMoment, T, C2H2, H, H2, H2O, O2, OH, M0,
     &	                   rateC2H2, rateO2, rateOH, rate, prodRateC2H2,
     &	                   prodRateCO, prodRateH, prodRateH2,
     &	                   prodRateH2O, prodRateO2, prodRateOH)

C	Calculate the surface growth / oxidation rates for moments and
C	the production rates for gas-phase species

C	Surface growth is currently described by the following scheme:
C	  1.  Csoot-H + H = Csoot* + H2         (fR1, rR1)
C	  2.  Csoot-H + OH = Csoot* + H2O	(fR2, rR2)
C	  3.  Csoot*  + H -> Csoot              (fR3)
C	  4.  Csoot*  + C2H2 -> Csoot-H + H     (fR4)
C	  5.  Csoot*  + O2 -> products          (fR5)
C	  6.  Csoot   + OH -> products          (fR6)

C	Input:
C	  highMoment:		the highest rate to be calculated
C	  T:			temperature (K)
C	  C2H2:			C2H2 concentration (mol cm^-3)
C	  H:			H concentration
C	  H2:			H2 concentration
C	  H2O:			H2O concentration
C	  O2:			O2 concentration
C	  OH:			OH concentration
C	  M0:			Zeroth moment
C	Output:
C	  rateC2H2:		moment rates due to C2H2 surface growth
C	  rateO2:		moment rates due to O2
C	  rateOH:		moment rates due to OH
C	  rate:			total moment rates due to surface growth
C	  prodRateC2H2:		production rate for C2H2
C	  prodRateCO:		production rate for CO
C	  prodRateH:		production rate for H
C	  prodRateH2:		production rate for H2
C	  prodRateH2O:		production rate for H2O
C	  prodRateO2:		production rate for O2
C	  prodRateOH:		production rate for OH
C				Negative rates indicate consumption

	implicit none

	include 'soot.h'

	integer highMoment
	doubleprecision T, C2H2, H, H2, H2O, O2, OH, M0,
     &	                rateC2H2(0:*), rateO2(0:*), rateOH(0:*),
     &	                rate(0:*), prodRateC2H2, prodRateCO, prodRateH,
     &	                prodRateH2, prodRateH2O, prodRateO2, prodRateOH

	doubleprecision mu(LO_FRAC:HI_FRAC)
	common /MOMENTS/ mu

	doubleprecision RT, fR1, rR1, fR2, rR2, fR3, fR4, fR5, fR6,
     &	                alpha, par_a, par_b, denom,
     &	                ssRatio, Carea, Crad, sC2H2, sO2, sOH,
     &	                D, dispersion
	integer n

C	The following rate coefficients are from Appel, Bockhorn, Frenklach 1999

	RT = 1.987D-3 * T        ! R in kcal/mole
	fR1 = 4.17D+13 * DEXP(-13.0 / RT) * H
	rR1 = 3.9D+12 * DEXP(-11.0 / RT) * H2
	fR2 = 1.0D+10 * (T**0.734D0) * DEXP(-1.43 / RT) * OH
	rR2 = 3.68E+08 * (T**1.139D0) * DEXP(-17.1 / RT) * H2O
	fR3 = 2.D+13 * H
	fR4 = 8.0D+7 * (T**1.56D0) * DEXP(-3.8 / RT) * C2H2
	fR5 = 2.2D+12 * DEXP(-7.5 / RT) * O2
	fR6 = 0.13D0 * OH	     ! gamma = 0.13 from Neoh et al.

	par_a = 12.65D0 - 5.63D-03 * T
	par_b = -1.38D0 + 6.80D-04 * T
	alpha = DTANH(par_a / DLOG10(mu(6)) + par_b)

C	Fraction of radical sites via steady-state:

	denom = rR1 + rR2 + fR3 + fR4 + fR5
	if(denom .ne. 0.0) then
	  ssRatio = (fR1 + fR2) / denom
	  Carea = alpha * surfC * M0
	  Crad  = Carea * ssRatio

	  rateC2H2(0) = 0.0
	  sC2H2 = fR4 * Crad
	  call sums(highMoment, 2, sC2H2, rateC2H2(1))

	  D = dispersion()

	  rateO2(0) = 0.0
	  sO2   = fR5 * Crad
	  call sums(highMoment, -2, sO2, rateO2(1))

	  rateOH(0) = 0.0
	  sOH   = fR6 * CBOH * sqrt(T) * M0
	  call sums(highMoment, -1, sOH, rateOH(1))

C	  Production rate for gas-phase species

	  prodRateC2H2 = -rateC2H2(1) / 2 / AVOGADRO
	  prodRateO2 =  rateO2(1) / 2 / AVOGADRO
	  prodRateOH =  ((rR2 * Crad - fR2 * Carea) * mu(4) +
     &	                                        rateOH(1)) / AVOGADRO
	  prodRateH =  ((rR1 + fR4) * Crad - fR1 * Carea - fR3 * Crad)
     &                    * mu(4) / AVOGADRO
	  prodRateH2 = (fR1 * Carea - rR1 * Crad) * mu(4) / AVOGADRO
	  prodRateH2O = (fR2 * Carea - rR2 * Crad) * mu(4) / AVOGADRO
	  prodRateCO = -(rateO2(1) + rateOH(1)) / AVOGADRO
	else
	  do n=0,highMoment
	    rateC2H2(n) = 0.0
	    rateO2(n) = 0.0
	    rateOH(n) = 0.0
	  enddo
	  prodRateC2H2 = 0.0
	  prodRateCO = 0.0
	  prodRateO2 = 0.0
	  prodRateOH = 0.0
	  prodRateH = 0.0
	  prodRateH2 = 0.0
	endif

	do n=0,highMoment
          rate(n) = rateC2H2(n) + rateO2(n) + rateOH(n)
	enddo

	return

	end

C**********************************************************************C

	double precision function viscosity(T)

C	Viscosity of air as a function of temperature (K)
C	(from AK, myf 1997)

	implicit none

	double precision T
	viscosity = 14.58D-06 * (T**1.5D0) / (T + 100.4D0)

	return

	end


