************************************************************************ * MASTER MECHANISM - ROUTINE NAME : ho_rad2 * * * * * * PURPOSE : * * perform reaction of OH with VOC, i.e. : * * --> abstraction of H from CH * * --> abstraction of H from (OH) hydroxy * * --> abstraction of H from (OOH) peroxy * * --> addition of HO to Cd=Cd * * * * Note : Activation temperatures are estimated from the rate * * constant function: k = A T**2 exp(-B/T) exp(Ex/T) * * where Ex is from substituent groups and A, B is from * * leaving "H" group * * * * In HO_RADical reaction module each C-H bond containing also al- * * cohol, and hydroperoxide group and double-bonded carbon in the * * chain is treated independently so that all possible abstraction * * of H and addition of HO may occur. * * For the reaction rate constants for H abstraction and for HO- * * addition see subroutines RABSOH and OHADD_Cx * * * * INPUT: * * - rdct : name[a6]+formula[a120] of the species for which * * reaction with OH is considered * * - bond(i,j) : carbon-carbon bond matrix of chem * * - group(i) : groups at position (carbon) i * * - nring : min number of distinct rings in molecule * * - temp : reference temperature (more work on this needed) * * - cut_off : ratio below which a pathway is not considered * * - dbrch : NOT USED - MORE WORK ON THIS NEEDED * * - level : number of level (stable + radicals) that were * * necessary to produce the parent of rdct * * - stabl : number of stable level (no radical) that were * * necessary to produce the parent of rdct * * - nfn : total nb. of species having a fixed name * * - namfn(i) : table of the fixed name (6 character) * * - chemfn(i) : formula corresponding the ith species having a * * fixed name * * * * Data for the reactions with OH (see oh_rate.dat) * * arrhenius parameter (see oh_rate.dat) * * - nohdat : total number of known OH rate constant * * - ohchem(i) : chemical formula of ith species having a known OH * * rate constant * * - ohrate(i,3) : arrhenius parameter of the known OH rate constant * * corresponding to the ith species * * * * known reaction (see oh_prod.dat) * * - nkwoh : number of VOC having known products for the * * reaction with ox * * - nkwohpd(i) : number of products for the * * reaction of the ith VOC with ox * * - kwohrct(i) : chemical formula of the VOC having a known * * chemistry for the reaction with ox * * - kwohpd(i,j) : chemical formula of the jth product for the * * reaction of the ith VOC with ox * * - kwohcopd(i,j,k): name (6 character) of the kth coproduct of the jth* * product for the reaction of the ith VOC with ox * * - kwohyld(i,j) : yield of the jth product for the reaction of the * * ith VOC with ox * * * * INPUT/OUTPUT * * - dict(j) : dictionary line (name + formula + functional * * group info) of species number j * * - namlst(j) : name (lco=6 characters) of the species already * * used at position number j * * - nhldvoc : number of (stable) VOC in the stack * * - holdvoc(i) : list of the VOC in the stack * * - nhldrad : number of radical in the stack * * - holdrad(i) : list of the radicals in the stack * * * LITERATURE REFERENCES * * * Atkinson, R, DL Baulch, RA Cox, RF Hampson, JA Kerr, MJ Rossi, J Troe, * Evaluated kinetic and photochemical data for atmospheric chemistry, * organic species: Supplement VII, J. Phys. Chem. Ref. Dat., 28, * 2,191-393, https://doi.org/10.1063/1.556048, 1999. * (POSSIBLY) * * Kwok, E.S.C., & R. Atkinson, Estimation of hydroxyl radical * reaction-rate constants for gas-phase organic-compounds using a * structure-reactivity relationship - an update,, Atm. Env., * 29,14,1685-1695, https://doi.org/ 10.1016/1352-2310(95)00069-B, 1995. * * Jenkin 2015, MAGNIFY ************************************************************************ SUBROUTINE ho_rad2(rdct,bond,group,nring,brch,temp, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nohdat,ohchem,ohrate, & nfn,namfn,chemfn, & nkwoh,nkwohpd,kwohrct,kwohpd,kwohcopd,kwohyld) IMPLICIT NONE INCLUDE 'general.h' INCLUDE 'organic.h' INCLUDE 'common.h' * input: CHARACTER(LEN=lcf),INTENT(in) :: rdct CHARACTER(LEN=lgr),INTENT(in) :: group(mca) INTEGER,INTENT(in) :: bond(mca,mca) INTEGER,INTENT(in) :: nring REAL,INTENT(in) :: brch REAL,INTENT(in) :: cut_off REAL,INTENT(in) :: temp INTEGER,INTENT(in) :: nohdat CHARACTER(LEN=lfo),INTENT(in) :: ohchem(mrd) REAL,INTENT(in) :: ohrate(mrd,3) INTEGER,INTENT(in) :: nkwoh,nkwohpd(mkr) CHARACTER(LEN=lfo),INTENT(in) :: kwohrct(mkr), kwohpd(mkr,mnr) CHARACTER(LEN=lco),INTENT(in) :: kwohcopd(mkr,mnr,mcp) REAL,INTENT(in) :: kwohyld(mkr,mnr) * input/output CHARACTER(LEN=ldi),INTENT(inout) :: dict(mni) REAL,INTENT(inout) :: dbrch(mni) CHARACTER(LEN=lco),INTENT(inout) :: namlst(mni) INTEGER,INTENT(inout) :: level INTEGER,INTENT(inout) :: stabl CHARACTER(LEN=lst),INTENT(inout) :: holdvoc(mlv) INTEGER,INTENT(inout) :: nhldvoc CHARACTER(LEN=lst),INTENT(inout) :: holdrad(mra) INTEGER,INTENT(inout) :: nhldrad INTEGER,INTENT(inout) :: nfn CHARACTER(LEN=lco),INTENT(inout) :: namfn(mfn) CHARACTER(LEN=lfo),INTENT(inout) :: chemfn(mfn) * internal: allow for up to 20 reaction channels CHARACTER(LEN=lfo) :: pchem(mnr), tempkc, pchema(mnr) CHARACTER(LEN=lfo) :: chem CHARACTER(LEN=lgr) :: tgroup(mca), pold, pnew, tempkg CHARACTER(LEN=lco) :: coprod(mnr,mca), tprod(mca) INTEGER :: tbond(mca,mca), flag(mnr),carb INTEGER :: nr,nch,ich,i,j,k,l,j1,np,nc,ncd REAL :: sktemp, ktemp(mnr), smk298, sum REAL :: k298(mnr), tact(mnr) REAL :: brtio REAL :: rk298 REAL :: ratio(mnr) REAL :: garrhc(3),tarrhc(mnr,3),rarrhc(3) INTEGER :: known_species,nohprod INTEGER :: idreac, nlabel REAL :: xlabel,folow(3),fotroe(4) CHARACTER*1 :: a1,a2,a3,a4 CHARACTER(LEN=lco) :: r(3), p(mnp) REAL :: s(mnp),ar1,ar2,ar3,f298,fratio,A,nt,E CHARACTER*12 :: resu CHARACTER(LEN=lco) :: copchem REAL :: rdtcopchem INTEGER :: cdtable(4),tcdtable(4),cdsub(4),conjug,ncdcase INTEGER :: cdcarbo(4,2),cdeth(4,2) INTEGER :: cpcdtable(4) INTEGER :: rjg(mri,2) ! ring-join group pairs INTEGER :: cnum,nca CHARACTER(LEN=lfo) :: rdckprod(mca),pchem_del(mnr) CHARACTER(LEN=lco) :: rdcktprod(mca,mca) INTEGER :: nip,flag_del(mnr) CHARACTER(LEN=lco) :: coprod_del(mnr,mca) REAL :: sc(mca),sc_del(mnr,mca) INTEGER :: rxnflg CHARACTER(LEN=lsb) :: progname='*ho_rad2*' CHARACTER(LEN=ler) :: mesg * write info for finding bugs IF (wtflag.NE.0) WRITE(*,*) progname,': input : ',rdct(lco+1:) * ----------- * initialize * ----------- copchem=' ' rdtcopchem=0. conjug=0 nc = 0 nca= 0 nch = 0 nr = 0 rk298=0. known_species=0 nohprod = 0 ktemp = 0 chem=' ' pchem = ' ' pchema = ' ' flag = 0 k298 = 0. tact = 0. ratio = 0. tgroup = group tbond = bond coprod = ' ' tarrhc=0. rarrhc=0. tcdtable=0 cdtable=0 coprod_del = ' ' pchem_del = ' ' flag_del = 0 sc_del = 0. rxnflg = 1 * IF RINGS EXIST remove ring-join characters from groups IF (nring.gt.0) CALL rjgrm(nring,group,rjg) tgroup = group tbond = bond i = lco + INDEX(rdct(lco+1:lcf),' ') chem=rdct(lco+1:i) nc = INDEX(chem,' ') - 1 IF (nc.LT.1) RETURN DO i=1,mca IF (tgroup(i).NE.' ') nca=nca+1 ENDDO * ------------------------------------------------------------ * check if the species is known in the rate constant data base * ------------------------------------------------------------ DO i=1,nohdat IF (chem.EQ.ohchem(i)) THEN rk298=ohrate(i,1)*298**ohrate(i,2)*exp(-ohrate(i,3)/298.) rarrhc(:)=ohrate(i,:) ENDIF ENDDO * ----------------------------------------------------------- * Check file ho_prod.dat for any known products & branching ratios * * WARNING: for those reactions, yields < 5% are not automatically * * eliminated * * ----------------------------------------------------------- * check radical and find co-products: c CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) c pchem(j) = rdckprod(1) c sc_del(j,1) = sc(1) c coprod(j,:) = rdcktprod(1,:) IF (rk298.NE.0.) THEN DO i=1,nkwoh IF (chem.EQ.kwohrct(i)) THEN known_species = 1 nohprod = nkwohpd(i) * loop through known products DO j=1,nohprod CALL addreac(nr,progname,rdct(lco+1:lcf),flag) ! nr = nr + 1 ! IF (nr.GT.mnr) THEN ! mesg = 'too many reactions created for species :' ! CALL errexit(progname,mesg,rdct(lco+1:lcf)) ! ENDIF ! flag(nr) = 1 pchema(j) = kwohpd(i,j) ratio(j) = kwohyld(i,j) * check the products IF (index(pchema(j),'.').ne.0) THEN IF(wtflag.NE.0)PRINT*,"radchk1" CALL radchk(pchema(j),rdckprod,rdcktprod,nip,sc) pchem(j)=rdckprod(1) IF (nip.EQ.2) THEN flag_del(nr) = 1 pchem_del(nr) = rdckprod(2) sc_del(nr,1) = sc(1) sc_del(nr,2) = sc(2) coprod_del(nr,:) = rdcktprod(2,:) CALL stdchm(pchem_del(nr)) ENDIF ***MC c coprod(j,1:mca) = tprod(1:mca) coprod(j,1:mca) = rdcktprod(1,1:mca) ***MC ELSE pchem(j)= pchema(j) ENDIF CALL stdchm(pchem(j)) * add known coproduct(s) if necessary ***MC c DO k=1,mcp ! loop possible coproducts c IF (kwohcopd(i,j,k).NE.' ') THEN c coprod(j,k) = kwohcopd(i,j,k) c ELSE c GOTO 450 c ENDIF ! coproduct is non-zero c ENDDO ! coproducts loop c c450 CONTINUE c ENDDO c ENDIF ! reactant has known products c ENDDO c IF (known_species.ne.0) GOTO 800 c ENDIF DO k=1,mcp ! loop possible coproducts IF (kwohcopd(i,j,k).NE.' ') THEN DO 455 l=1,mca IF (coprod(j,l).ne.' ') GOTO 455 coprod(j,l) = kwohcopd(i,j,k) GOTO 456 455 CONTINUE ELSE GOTO 450 ENDIF ! coproduct is non-zero 456 CONTINUE ENDDO ! coproducts loop 450 CONTINUE ENDDO ENDIF ! reactant has known products ENDDO IF (known_species.ne.0) GOTO 800 ENDIF ***MC ************************************************************************ **** H ABSTRACTION *** ************************************************************************ * --------------------------------- * H-ABSTRACTION FROM CH * --------------------------------- * find all C-H bonds - loop over each C DO i=1,nca IF (INDEX(group(i),'CH').EQ.0) CYCLE CALL addreac(nr,progname,rdct(lco+1:lcf),flag) * compute associated rate constant: k=A*T**2*exp(-E/T) garrhc(:)=0. CALL rabsoh(tbond,tgroup,i,garrhc,nring) tarrhc(nr,:)=garrhc(:) * remove one H atom IF (group(i)(1:3).EQ.methyl) THEN pold = methyl pnew = primary ELSE IF (group(i)(1:3).EQ.primary) THEN pold = primary pnew = secondary ELSE IF (group(i)(1:3).EQ.aldehyde) THEN pold = aldehyde pnew = carbonyl ELSE IF (group(i)(1:2).EQ.secondary) THEN pold = secondary pnew = 'C' ENDIF CALL swap(group(i),pold,tgroup(i),pnew) * add radical dot at end of group: ! CALL adddot(tgroup(i)) ! (alternative: shorter) nc = INDEX(tgroup(i),' ') tgroup(i)(nc:nc) = '.' * re-build molecule, including ring if present: CALL rebond(tbond,tgroup,tempkc,nring) * check radical and find co-products: IF(wtflag.NE.0)PRINT*,"radchk2" CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) pchem(nr) = rdckprod(1) sc_del(nr,1) = sc(1) IF (nip.EQ.2) THEN flag_del(nr) = 1 pchem_del(nr) = rdckprod(2) sc_del(nr,2) = sc(2) coprod_del(nr,:) = rdcktprod(2,:) CALL stdchm(pchem_del(nr)) ENDIF coprod(nr,:) = rdcktprod(1,:) * find standard name: CALL stdchm(pchem(nr)) * reset group: tgroup(i) = group(i) ENDDO *------------------------------------------------ * H ABSTRACTION FROM (OH) - BUT NOT CO(OH): *------------------------------------------------ * find all hydroxy group but not carboxylic acid DO 120 i=1,mca IF (INDEX(group(i),hydroxy).EQ.0) GO TO 120 IF (INDEX(group(i),carboxylic_acid).NE.0) GO TO 120 CALL addreac(nr,progname,rdct(lco+1:lcf),flag) * assign rate constant by Kwok and Atkinson, Atm. Env., 1695, 1995 ! Jenkin 2015, MAGNIFY tarrhc(nr,1)=1.28E-12 tarrhc(nr,2)=0. tarrhc(nr,3)=660. * replace hydroxy group by alkoxy group: pold = hydroxy pnew = alkoxy CALL swap(group(i),pold,tgroup(i),pnew) * rebuild, check, and find coproducts: CALL rebond(tbond,tgroup,tempkc,nring) IF(wtflag.NE.0)PRINT*,"radchk3" CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) pchem(nr) = rdckprod(1) IF (nip.NE.1) WRITE(99,*) '2 produits ho_rad2.f 2' IF (nip.NE.1) STOP coprod(nr,:) = rdcktprod(1,:) * rename CALL stdchm(pchem(nr)) * reset group: tgroup(i) = group(i) 120 CONTINUE *------------------------------------------------ * H-ATOM ABSTRACTION FROM (OOH) *------------------------------------------------ * find hydroperoxide groups: (Peroxy acid is not treated separately!) DO 130 i=1,mca IF (INDEX(group(i),hydro_peroxide).EQ.0) GO TO 130 CALL addreac(nr,progname,rdct(lco+1:lcf),flag) * assign rate constants for alkyl hydroperoxide using the * rate constant for the CH3OOH+OH -> CH3O2 provided by * Atkinson, J. Phys. Chem. Ref. Data, 521, 1999 tarrhc(nr,1)=1.9E-12 tarrhc(nr,2)=0. tarrhc(nr,3)=-190. * change (OOH) to (OO.) pold = hydro_peroxide pnew = alkyl_peroxy CALL swap(group(i),pold,tgroup(i),pnew) * rebuild, check, and assign coproducts: CALL rebond(tbond,tgroup,tempkc,nring) IF(wtflag.NE.0)PRINT*,"radchk4" CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) pchem(nr) = rdckprod(1) IF (nip.NE.1) WRITE(6,*) '2 produits ho_rad2.f 3' coprod(nr,:) = rdcktprod(1,:) * rename: CALL stdchm(pchem(nr)) * reset group: tgroup(i) = group(i) 130 CONTINUE *------------------------------------------------ * OH ADDITION TO THE -ONO2 FUNCTIONAL GROUP *------------------------------------------------ * REACTION REMOVED - NOT FOUND IN RECENT REVIEW ! * See archive of the generator (before april 2004) if this reaction * need to be introduced again in the generator *------------------------------------------------ * H ABSTRACTION (OH ADDITION ?) TO THE -COOH GROUP *------------------------------------------------ * IF RINGS EXIST remove ring-join characters from groups IF (nring.gt.0) CALL rjgrm(nring,group,rjg) * find carboxylic groups: DO 430 i=1,mca IF (INDEX(group(i),carboxylic_acid).EQ.0) GO TO 430 CALL addreac(nr,progname,rdct(lco+1:lcf),flag) * assign rate constants for OH addition using the rate constant of Kwok * and Atkinson, Atm. Env., 1695, 1995. reaction lead to H abstraction, * i.e. : RCOOH + OH -> RCOO + H2O tarrhc(nr,1)=0.52E-12 tarrhc(nr,2)=0. tarrhc(nr,3)=0. * change CO(OH) to CO(O.) pold = carboxylic_acid pnew = acyl_oxy CALL swap(group(i),pold,tgroup(i),pnew) *rebuild, check, and assign coproducts: CALL rebond(tbond,tgroup,tempkc,nring) IF(wtflag.NE.0)PRINT*,"radchk5" CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) pchem(nr) = rdckprod(1) c IF (nip.NE.1) WRITE(6,*) '2 produits ho_rad2.f 4' sc_del(nr,1) = sc(1) IF (nip.EQ.2) THEN flag_del(nr) = 1 pchem_del(nr) = rdckprod(2) sc_del(nr,2) = sc(2) coprod_del(nr,:) = rdcktprod(2,:) ENDIF coprod(nr,:) = rdcktprod(1,:) * add H2O to the list of coproducts * DO j=1,mca * IF (coprod(nr,j)(1:2).EQ.' ') THEN * coprod(nr,j)='H2O' * GOTO 440 * ENDIF * ENDDO *440 CONTINUE * rename: CALL stdchm(pchem(nr)) * reset group: tgroup(i) = group(i) 430 CONTINUE ************************************************************************ **** OH ADDITION TO >C=C< BOND *** ************************************************************************ * IF RINGS EXIST remove ring-join characters from groups IF (nring.gt.0) THEN CALL rjgrm(nring,group,rjg) ENDIF * find double-bonded carbons in the chain (label 600 => end C=C) IF (INDEX(rdct(lco+1:lcf),'Cd').EQ.0) GOTO 600 IF (wtflag.NE.0) WRITE(*,*) & "*oh addition to double bond* ------------------------- " * The parameterization used to assign a rate constant for OH addition * to C=C bonds depends on whether the C=C bond is conjugated or not * with a C=O bond (i.e. structure of type -C=C-C=O). The subroutine * return the "case" the species belong. Five cases are considered : * * CASE 1 : regular Cd molecule : only >C=C< and >C=C-C=C< bonds in * the molecule (i.e, without conjugated C=C-C=O) * CASE 2 : are for structure containing the >C=C-C=O structure * but no C=C-C=C * CASE 3 : are for the -CO-C=C-C=C-C=O structure only (i.e. containing * carbonyl at both side of the conjugated C=C-C=C) * CASE 4 : Two double bonds non conjugated (i.e. C=C-C-C=C) but one * containing at least one C=C-C=O * CASE 5 : are for the -CO-C=C-C=C< structure (i.e. containing * carbonyl at only one side of the conjugated C=C-C=C * Ludo: * CASE 7 : addition of OH to C=C bond of vinyl ether (-O-C=C) * and dihydrofurans * get the case for chem (cdcase) CALL cdcase(chem,bond,group,rxnflg,nring, & ncd,ncdcase,conjug, & cdtable,tcdtable,cdsub,cdcarbo,cdeth) * case 1 : C=C bonds conjugated or not, without any C=C-C=O structure * ------- IF (ncdcase.eq.1) THEN CALL hoadd_c1(chem,bond,group,ncd,cdtable,conjug, & nr,flag,tarrhc,pchem,coprod,flag_del, & pchem_del,coprod_del,sc_del) * case 2 : C=C-C=O bonds, only 1 C=C (i.e no conjugated Cds) * ------- ELSE IF (ncdcase.eq.2) THEN CALL hoadd_c2(chem,bond,group, & ncd,conjug,cdtable,cdsub,cdcarbo, & nr,flag,tarrhc,pchem,coprod) * case 3 : -C=O-C=C-C=C-C=O- structure * ------- ELSE IF (ncdcase.eq.3) THEN CALL hoadd_c3(chem,bond,group, & ncd,conjug,cdtable,tcdtable,cdsub,cdcarbo, & nr,flag,tarrhc,pchem,coprod,flag_del, & pchem_del,coprod_del,sc_del) * case 4 : -C=C-C-C=C-C=O- * ------- ELSE IF (ncdcase.eq.4) THEN * perform reaction for each C=C, ignoring the other C=C bond DO i=1,3,2 * check everything is ok IF (cdtable(i).eq.0) THEN mesg = 'case 4 encountered but Cd not found for' CALL errexit(progname,mesg,chem) WRITE(99,*) 'ho_rad2',chem !STOP ENDIF * make a copy of cdtable and "kill" the Cd bond not considered cpcdtable = 0 cpcdtable(i)=cdtable(i) cpcdtable(i+1)=cdtable(i+1) * send to reaction subroutine IF (cdcarbo(i,1).ne.0) THEN CALL hoadd_c2(chem,bond,group, & ncd,conjug,cpcdtable,cdsub,cdcarbo, & nr,flag,tarrhc,pchem,coprod) ELSE IF (cdcarbo(i+1,1).ne.0) THEN CALL hoadd_c2(chem,bond,group, & ncd,conjug,cpcdtable,cdsub,cdcarbo, & nr,flag,tarrhc,pchem,coprod) ELSE CALL hoadd_c1(chem,bond,group,ncd,cpcdtable,conjug, & nr,flag,tarrhc,pchem,coprod,flag_del, & pchem_del,coprod_del,sc_del) ENDIF ENDDO c CALL hoadd_c4(chem,bond,group,ncd,tcdtable, c & nr,flag,tarrhc,pchem,coprod) * case 5 : -CO-C=C-C=C< * ------- ELSE IF (ncdcase.eq.5) THEN CALL hoadd_c5(chem,bond,group, & ncd,conjug,cdtable,cdcarbo, & nr,flag,tarrhc,pchem,coprod,flag_del, & pchem_del,coprod_del,sc_del) * case 6 : -C=C=O * ------- ELSE IF (ncdcase.eq.6) THEN CALL hoadd_c6(chem,bond,group, ncd, & nr,flag,tarrhc,pchem,coprod) * Case 7 : -C=C-O- * ---------------- ELSE IF (ncdcase.EQ.7) THEN CALL hoadd_c7(chem,bond,group,ncd,cdtable,cdeth,conjug, & nr,flag,tarrhc,pchem,coprod) ENDIF * end of OH addition to C=C 600 CONTINUE !------------------------------------------------ ! OH addition on aromatic ring !----------------------------------------------- IF (INDEX(chem,'c').NE.0) THEN CALL hoadd_arom(chem,bond,group,nca, & nr,flag,tarrhc,pchem,coprod) ENDIF 650 CONTINUE *------------------------------------------------ * END REACTIONS -SELECT THE VARIOUS PATHWAYS *------------------------------------------------ * compute the rate at 298 K (do not exceed collision rate) sktemp=0. smk298=0. DO i=1,nr k298(i) = tarrhc(i,1) * 298.**tarrhc(i,2) * & EXP(-tarrhc(i,3)/298.) k298(i) = AMIN1(k298(i),2.0E-10) ktemp(i) = tarrhc(i,1) * temp**tarrhc(i,2) * & EXP(-tarrhc(i,3)/temp) sktemp = sktemp + ktemp(i) smk298 = smk298 + k298(i) ENDDO * preliminary version : remove the temperature dependency for C-H * abstraction reaction (except -OH (Ea=85) and -OOH(Ea=-190) ) DO i=1,nr IF (tarrhc(i,3).ne.-190.) THEN IF (tarrhc(i,3).ne.85.) THEN tarrhc(i,:)=0. tarrhc(i,1)=k298(i) ENDIF ENDIF ENDDO * if no reactions, warning error (this happened for CHO-O-CHO for * example, for which the rate constant is expected to be low but not * available - abstraction forms -O-CHO = 0) IF (sktemp.EQ.0) THEN WRITE(99,'(a)') '--warning-- (STOP)' WRITE(99,'(a)') 'OH rate constant estimated as 0 in ho_rad' WRITE(99,'(a)') 'reaction with OH is not considered for:' WRITE(99,'(a)') chem RETURN ENDIF * if total rate constant less than 1.0E-13, then warning IF (sktemp.LT.1.0E-13) THEN WRITE(99,'(a)') '--warning--' WRITE(99,'(a)') 'OH rate constant estimated less than 1E-13' WRITE(99,'(a)') 'reaction with OH is not considered for:' WRITE(99,'(a)') chem ENDIF * keep only 5 numbers after comma before compare k298 DO i=1,nr resu=' ' WRITE(resu,'(e12.5)') k298(i) READ (resu,'(e12.5)') k298(i) ENDDO * lump what is identical in the product, coproduct, and rate constant IF (nr.ge.2) THEN DO 200 i=1,nr-1 DO 210 j=i+1,nr IF (pchem(i).EQ.pchem(j)) THEN IF (k298(i).ne.k298(j)) GOTO 210 DO k=1,mca IF (coprod(i,k).ne.coprod(j,k)) GO TO 210 ENDDO IF (flag(i).eq.0) GOTO 210 flag(j) = 0 flag(i)=flag(i)+1 ENDIF 210 CONTINUE 200 CONTINUE ENDIF WHERE(flag.gt.1) k298 = k298 * REAL(flag) ktemp = ktemp * REAL(flag) tarrhc(:,1) = tarrhc(:,1) * REAL(flag) flag=1 ENDWHERE * write estimated versus prescribed rate constants IF ((wtflag.NE.0).AND.(rk298.NE.0.)) & WRITE(73,*) rk298,smk298,' ',chem(1:index(chem,' ')) * remove pathways below cutoff WHERE (k298.LT.cut_off*smk298) flag = 0 sum=0. DO i=1,nr IF (flag(i).EQ.1) THEN nch = nch + 1 sum = sum + k298(i) ENDIF ENDDO * set branching ratio for each channel DO i=1,nr IF (flag(i).EQ.1) THEN ratio(i)=k298(i)/sum ENDIF ENDDO * reentry point if both rate constant and reaction products are known 800 CONTINUE * set the reaction rate of each channel * prelimary version : if the rate constant is assigned then multiply * the first arrhenius coef. by the branching ratio IF (rk298.NE.0.) THEN IF(rarrhc(1).GT.0.) WRITE(72,*) rdct(1:6),' ',rarrhc(:) & ,' ',chem(1:index(chem,' ')) DO i=1,nr IF (flag(i).EQ.1) THEN tarrhc(i,1)=rarrhc(1)*ratio(i) tarrhc(i,2)=rarrhc(2) tarrhc(i,3)=rarrhc(3) ENDIF ENDDO * if the rate constant is unknown then rescale tarrhc(i,1) ELSE IF(smk298.GT.0.) WRITE(72,*) rdct(1:6),' ',smk298,rarrhc(2:3) & ,' ',chem(1:index(chem,' ')) DO i=1,nr IF (flag(i).EQ.1) THEN tarrhc(i,1)=tarrhc(i,1)*smk298/sum ENDIF ENDDO ENDIF * remove pathways below cutoff * ------------------------ * OUTPUT - WRITE REACTION * ------------------------ nt =0. E = 0. ich = 10 A = 0. * write information required for operator IF (wtopeflag.EQ.1) THEN DO 11 i=1,mnr IF (flag(i).NE.0) THEN A = A + tarrhc(i,1) IF ((tarrhc(i,2).NE.nt).OR.(tarrhc(i,3).NE.E)) THEN A = sum nt = 0. E = 0. GOTO 11 ELSE nt=tarrhc(i,2) E=tarrhc(i,3) ENDIF ENDIF 11 CONTINUE WRITE(10,'(A15,A1,A6,A4,1X,ES10.3,1X,f4.1,1X,f7.0)') & '****INIT GHO +','G',rdct(1:lco),'****',A,nt,E ENDIF DO 300 i=1,nr IF (flag(i).EQ.0) GO TO 300 * initialize reaction CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ich = ich + 1 IF (nch.GT.1) a4 = alfa(ich:ich) a1 = rdct(1:1) ar1=tarrhc(i,1) ar2=tarrhc(i,2) ar3=tarrhc(i,3) IF (rk298.NE.0.) THEN f298=rk298 ELSE f298=smk298 ENDIF fratio=ratio(i) r(1) = rdct(1:lco) r(2) = 'HO ' s(1) = 1. IF (known_species.eq.0) THEN brtio = brch * ktemp(i)/sktemp ELSE brtio = brch * ratio(i) ENDIF * CALL bratio(pchem(i),brtio,p(1),copchem,rdtcopchem, * print*,'i, ratio = ',i,k298(i)/smk298 CALL bratio(pchem(i),brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) IF (wtopeflag.EQ.1) !write information required for operator & WRITE(10,'(f5.3,2X,A1,A6)') ratio(i), 'G',p(1) * second product (if delocalisation allowed) np = 1 IF (flag_del(i).NE.0) THEN np = np + 1 s(1) = sc_del(i,1) s(np) = sc_del(i,2) CALL bratio(pchem_del(i),brtio,p(np), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ENDIF * write out coproducts c np = 1 DO j=1,mca IF (coprod(i,j)(1:1).NE.' ') THEN np = np + 1 IF (np.GT.mnp) THEN mesg = 'np > mnp (too many products in the reaction : 1)' CALL errexit(progname,mesg,chem) ENDIF s(np) = 1.0 p(np) = coprod(i,j) IF (wtopeflag.EQ.1) & WRITE(10,'(f5.3,2X,A1,A6)')s(np)*ratio(i),'G',p(np) ENDIF ENDDO DO j=1,mca IF ((coprod_del(i,j)(1:1).NE.' ').AND. & (flag_del(i).NE.0)) THEN np = np + 1 IF (np.GT.mnp) then mesg = 'np > mnp (too many products in the reaction : 2)' CALL errexit(progname,mesg,chem) ENDIF s(np) = sc_del(i,j) p(np) = coprod_del(i,j) IF (wtopeflag.EQ.1) write(10,'(20X,A1,A6)') 'G',p(np) ENDIF ENDDO IF (rdtcopchem.GT.0.) THEN np = np + 1 IF (np.GT.mnp) THEN mesg = 'np > mnp (too many products in the reaction : 3)' CALL errexit(progname,mesg,chem) ENDIF s(np) = rdtcopchem p(np) = copchem IF (wtopeflag.EQ.1) & WRITE(10,'(f5.3,2X,A1,A6)')s(np)*ratio(i),'G',p(np) ENDIF * write out: * reactions with OH are thermal reactions (idreac=0) and do not require labels CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) 300 CONTINUE IF (wtopeflag.EQ.1) WRITE(10,*)'end',progname * end of ho_rad 999 CONTINUE !IF (wtflag.EQ.1) WRITE(*,*)progname,': end :' RETURN END