c******************************************************************************* c ... evaluate rate of progress of a given reaction kr function rop(kr,tc,tempref,cspl) IMPLICIT NONE include 'parameter.par' include 'chem.com' real*8 cgspl(nsmx),wspl(nsmx),cspl(nsmx),ropl(idim), 1 qfac(idim),rr_f(idim),hspl(nsmx) double precision tc,tempref,one,otc,vlntemp,ctot, 1 rr_r,xik,crprod,cvp,cpprod,vt,vot,cthb, 2 rklow,pr,prlog,fc,xp,fcent,fclog,xn, 3 cprlog,flog,clog,pcor,rop integer k,kr,ksp,isp,ip3,nyp,nnn,kthb,ispthb,kfl,kfc,ifp c------------------------------------------------------------------------------- c ... prelude one = 1.0 k = kr otc = 1./tc vlntemp = log(tc) call cgsp(tc,cgspl) ! gibbs free energies ctot = 0.0 do ksp=1,ns ctot = ctot + cspl(ksp) end do c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... initialize qfac array, product of concentrations to be used in c ... third body and fall off reactions qfac(k) = 1.0 c------------------------------------------------------------------------------- c ... evaluate forward reaction rate, rr_f if(itfdep(k).ne.0)then rr_f(k) = capb(k)*exp(par(2,k)*vlntemp - cape(k)*otc) else rr_f(k) = capb(k) endif c------------------------------------------------------------------------------- c ... nspec(k) > 0 iff reaction is reversible c ... xik = SUM_i=1,ns nu(i,k)*G_ik c ... evaluate (1/equilibrium_reaction_rate) = orr_e c ... evaluate reverse reaction rate, rr_r c ... kflagrev(k) = 1 -> reverse Arrhenius rate parameters are defined c ... use these rather than orr_e to find rr_r if(nspec(k).gt.0)then if(kflagrev(k).eq.1)then if(itrdep(k).ne.0)then rr_r = rcapb(k)*exp(rexpo(k)*vlntemp-rcape(k)*otc) else rr_r = rcapb(k) endif else xik = 0.0 do isp=1,nreac(k) xik = xik + anu(isp,k)*cgspl(nunk(isp,k)) end do do ip3=4,nprod(k)+3 xik = xik + anu(ip3,k)*cgspl(nunk(ip3,k)) end do rr_r = rr_f(k)*exp(xik*otc+anusum(k)*vlntemp) endif else rr_r = 0.0 ! reaction irreversible-> rr_r=0 endif c------------------------------------------------------------------------------- c ... find reactants concentrations product PROD(c_j**nu_jk^prime)_j=1,ns c ... it is assumed below that nyp is a positive integer (non-zero) c ... which should be true in general crprod = 1.0 do isp = 1,nreac(k) nyp = - nu(isp,k) cvp = cspl(nunk(isp,k)) do nnn=1,nyp crprod = crprod * cvp end do end do c------------------------------------------------------------------------------- c ... for reverse reaction, c ... find products concentrations product PROD(c_j**nu_jk^double-prime)_j=1,ns c ... nspec(k) > 0 if reaction is reversible c ... it is assumed below that nu(ip3,k) is a positive integer (non-zero) c ... which should be true in general if(nspec(k).gt.0)then cpprod = 1.0 do ip3 = 4,nprod(k)+3 cvp = cspl(nunk(ip3,k)) do nnn=1,nu(ip3,k) cpprod = cpprod * cvp end do end do else cpprod = 0.0 endif c------------------------------------------------------------------------------- c ... evaluate reaction rate-of-progress ropl(k) = rr_f(k)*crprod - rr_r*cpprod c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... evaluate third-body-scaled reaction rates-of-progress c ... nb. qfac is the product of concentrations, with possible 3rd c ... body efficiencies. c ... if the reaction is also a fall-off reaction with kfal=0, then c ... qfac is just used to find pr below, and then reset to 1.0 c ... otherwise, it is actually used to scale the rate of progress c ... later on below if(nthb.gt.0)then kthb = jthb(kr) if(kthb.gt.0)then qfac(kr) = ctot do ispthb = 1,ntbs(kthb) ksp = nktb(ispthb,kthb) qfac(kr) = qfac(kr)+(aik(ispthb,kthb)-1.d0)*cspl(ksp) end do endif endif c------------------------------------------------------------------------------- c ... evaluate fall-off pressure corrections c ... vt is dimensional temperature (K), vot = 1/vt c ... these are used in the computation of the dimensionless F if(nfal.gt.0)then kfl = jfal(kr) if(kfl.gt.0)then vt = tc*tempref vot = 1.0/vt kfc = kfal(kfl) ifp = ifop(kfl) if(kfc.eq.0)then cthb = qfac(kr) qfac(kr)= 1.0 else cthb = cspl(kfc) endif rklow = capbl(kr)*exp(fpar(2,kfl)*vlntemp-capel(kr)*otc) pr = rklow * cthb / rr_f(kr) prlog = log10(pr) if(ifp.eq.1)then ! Lindemann form fc = 1.0 else ! ifop = 2, 3, or 4 if(ifp.eq.2)then ! SRI form xp = 1.0/(1.0+prlog*prlog) fc = ( + (fpar(4,kfl) * exp(-fpar(5,kfl)*vot) + + exp(-vt/fpar(6,kfl)) + )**xp + ) * fpar(7,kfl) * vt**fpar(8,kfl) else ! 6 or 7-param TROE form fcent = (1.0-fpar(4,kfl))*exp(-vt/fpar(5,kfl)) + + fpar(4,kfl)*exp(-vt/fpar(6,kfl)) if(ifp.eq.4)then ! 7-parameter TROE form fcent = fcent + exp(-fpar(7,kfl)*vot) endif fclog = log10(fcent) xn = 0.75 - 1.27*fclog cprlog = prlog - (0.4+0.67*fclog) clog = cprlog/(xn-0.14*cprlog) flog = fclog/(1.0+clog*clog) fc = 10.0**flog endif endif pcor = fc*pr/(1.0+pr) ropl(kr) = ropl(kr) * pcor endif endif c------------------------------------------------------------------------------- c ... modify ropl due to 3rd bodies, given possible mod to qfac bec falloff if(nthb.gt.0)then kthb = jthb(kr) if(kthb.gt.0)then ropl(kr) = qfac(kr) * ropl(kr) endif endif c------------------------------------------------------------------------------- rop = ropl(kr) return end c*******************************************************************************