! MASTER MECHANISM V.3.0 ROUTINE NAME - RO * ! * ! ! -------------------------------------------------------- ! ! SAR used for the rate constants estimates : ! ! - Isomerization : Choice between Atkinson, 2007 (Atmospheric ! ! Environment) and Vereecken and Peeters (2010) ! ! - Decomposition : Choice between Atkinson, 2007 (Atmospheric ! ! Environment) and Vereecken and Peeters (2009) ! ! - RO+O2: Values from Atkinson, 2007 (Atmospheric Environment) ! ! ------------------------------------------------------- ! ! ! INCLUDE general.h includes all information about glo- * ! variables: CUT_OFF, ALFA, and func- * ! tional groups * ! ! ! COMMON BLOCK CHMLST - Information about the chemicals in * ! the dictionary * ! CHMDAT - Information about the groups and the * ! bond-matrix of the chemical in RDCT * ! * ! LOCAL CONSTANTS... * ! * ! N maximum number of products allowed in HVDISS * ! * ! LOCAL VARIABLES... * ! * ! INTERNAL: * ! * ! TGROUP contents the groups of RDCT * ! TBOND contents the bond-matrix of RDCT * ! PCHEM temporary variable names for products * ! FLAG,NR,NP flag indicates how many products (NR) exit * ! IND index of chemical in CHMLST - is 0, if not ex. * ! NDB number of double bonds in chain * ! NCH, ICH number and index of non-identical products * ! I,J,K DO-LOOPs indices * ! RATE reaction rates of product channels * ! group: 2.5e-15 molecule per sec * ! TOTR total sum of rates of reactions with NO3 & RDCT * ! RTOT total rate for reaction of double-bonded carbon * ! FRACT fraction deduced by position of d-bonded carbon * ! PP,SS products and stoichiometric coefficients after * ! fragmentation of excited Criegee radicals * ! * ! OUTPUT: * ! * ! A1 & A4 information about type and channel of reaction * ! R(1-2) reagent: input chemical and NO3 * ! S(N) stoichiometric coefficients of products * ! P(N) array of products: acyl radicals, HNO3 ,and co- * ! products * ! TA activation energy is 0. * ! * ! LITERATURE REFERENCES * ! * ! Atkinson, R., Rate constants for the atmospheric reactions of ! alkoxy radicals: An updated estimation method, Atmos. Environ., ! 41,38,8468-8485, https://doi.org/10.1016/j.atmosenv.2007.07.002, 2007 ! SAPRC99 ! Tuazon (98) (qualitative data) ***************************************************************** SUBROUTINE ro(rdct,bond,group,nring,brch,temp,tbm, & nbson,bsongrp,bsonval, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn, & nrokr,rokrct,rokprd,roarrh,rocost) IMPLICIT NONE INCLUDE 'general.h' INCLUDE 'organic.h' INCLUDE 'common.h' ! photolysis reactions: ! input: CHARACTER(LEN=lcf) rdct CHARACTER(LEN=lgr) group(mca) INTEGER bond(mca,mca) INTEGER nring REAL brch REAL cut_off REAL temp REAL tbm INTEGER nbson CHARACTER(LEN=lgb) bsongrp(mbg) REAL bsonval(mbg) CHARACTER(LEN=lfo) rokrct(mkr,2),rokprd(mkr,3) INTEGER nrokr REAL roarrh(mkr,3),rocost(mkr,3) ! input/output CHARACTER(LEN=ldi) dict(mni) REAL dbrch(mni) CHARACTER(LEN=lco) namlst(mni) INTEGER level INTEGER stabl CHARACTER(LEN=lst) holdvoc(mlv) INTEGER nhldvoc CHARACTER(LEN=lst) holdrad(mra) INTEGER nhldrad INTEGER nfn CHARACTER(LEN=lco) namfn(mfn) CHARACTER(LEN=lfo) chemfn(mfn) ! internal REAL heat CHARACTER(LEN=lfo) pest(2),po2, piso(mca), pdec(mca,2), tempkc CHARACTER(LEN=lfo) tempo2,tempdec CHARACTER(LEN=lgr) tgroup(mca), tempkg, pold, pnew CHARACTER(LEN=lco) copiso(mca,mca), copdec(mca,mca), tprod(mca) CHARACTER(LEN=lco) copest(mca),copo2(mca) INTEGER tbond(mca,mca), flo2, fliso(mca), fldec(mca) INTEGER np,nc,ich,nich,ia,ia0,ii,i,j,k,l,m,n, flest INTEGER imol,irad,itype,niso,ndec,icount INTEGER oring, oringa,oringb, oringc INTEGER specialcase INTEGER icase, ilen REAL brtio,sum_o2,ea,Eaa,Eb REAL dhfro, dhfp1, dhfp2, dhnet,dhfho2 REAL rsum,rsum298 REAL ar1est,ar2est,ar3est REAL ar1o2,ar2o2,ar3o2 REAL ar1iso(mca),ar2iso(mca),ar3iso(mca) REAL ar1dec(mca),ar2dec(mca),ar3dec(mca) REAL kiso(mca), kdec(mca), kest CHARACTER*1 a1, a2, a3, a4 CHARACTER(LEN=lco) p(mnp), r(3) REAL s(mnp), ar1,ar2,ar3,f298,fratio CHARACTER*12 resu INTEGER idreac, nlabel REAL xlabel,folow(3),fotroe(4) CHARACTER(LEN=lco) copchem,copchem1,copchem2 REAL rdtcopchem,rdtcopchem1,rdtcopchem2 INTEGER nca,ncr,rjg(mri,2),ring(mca) INTEGER begrg,endrg,rngflg,dbflg,opflg CHARACTER(LEN=lfo) rdckprod(mca),pdec2(mca),piso2(mca) CHARACTER(LEN=lco) rdcktprod(mca,mca),copdec2(mca,mca) CHARACTER(LEN=lco) copiso2(mca,mca) INTEGER nip INTEGER fldec2(mca),fliso2(mca),fldecether REAL sc(mca),siso(mca,mca),sdec(mca,mca) INTEGER nrg,dbf,bddec(mca,mca),ring_old CHARACTER(LEN=lgr) grdec(mca) CHARACTER(LEN=lfo) chem,tchem(2) INTEGER known_species REAL :: FSD_tmp(5),FSD(mca,5) CHARACTER(LEN=lsb) :: progname='*ro! ' CHARACTER(LEN=ler) :: mesg ! heat of formation for HO2 (in kcal.mol-1) DATA dhfho2/3.49/ IF (wtflag.NE.0) write(6,*) '*ro* : ',rdct(lco+1:lcf) ! check if species is allowed in this routine IF (INDEX(rdct(lco+1:lcf),alkoxy).EQ.0) THEN WRITE(6,'(a)') '--error--, in subroutine: ro' WRITE(6,'(a)') 'this routine was called with' WRITE(6,'(a)') 'a species with no R(O.) function:' WRITE(6,'(a)') rdct(lco+1:lcf) WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF ! initialize reaction line and transcribe groups and bonds nca=0 DO i=1,mca IF (group(i)(1:1).NE.' ') nca=nca+1 ENDDO copchem= ' ' copchem1=' ' copchem2=' ' rdtcopchem=0. rdtcopchem1=0. rdtcopchem2=0. specialcase=0 po2=' ' pnew=' ' piso(:)=' ' pdec(:,:)=' ' tprod(:)=' ' pdec2(:)=' ' piso2(:)=' ' rdcktprod(:,:)=' ' copdec(:,:)=' ' copdec2(:,:)=' ' copiso(:,:)=' ' copest(:)= ' ' copo2(:)=' ' siso(:,:) = 0 sdec(:,:) = 0 flest = 0 fldecether = 0 flo2 = 0 niso = 0 ndec = 0 tempo2 =' ' fliso(:) = 0 fldec(:) = 0 fldec2(:)=0 fliso2(:)=0 ar1iso(:) = 0. ar2iso(:) = 0. ar3iso(:) = 0. ar1dec(:) = 0. ar2dec(:) = 0. ar3dec(:) = 0. kiso(:)=0. kdec(:)=0. Eaa=0. ea=0. dhfp1=0. dhfp2=0. dhnet=0. ring(:)=0 ar1est = 0 ar2est = 0 ar3est = 0 pest(:) = ' ' ar1o2 = 0. ar2o2 = 0. ar3o2 = 0. dbflg=0 opflg=0 rngflg=0 tempkc=' ' nrg=0 dbf=0 ring_old=0 grdec(:)=' ' bddec(:,:)=0 known_species = 0 FSD_tmp(:)=0 FSD(:,:)=0 ********************************************************************** ! check if the reaction is known in the database i = lco + INDEX(rdct(lco+1:lcf),' ') chem=rdct(lco+1:i) nc = INDEX(chem,' ') - 1 IF (nc.LT.1) RETURN DO i=1,nrokr IF (chem.EQ.rokrct(i,1)) THEN known_species = 1 WRITE(6,*) 'REACTION DE RO CONNU' ! WRITE REACTION CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ! reactant a1 = '1' r(1) = rdct(1:lco) ar1 = roarrh(i,1) ar2 = roarrh(i,2) ar3 = roarrh(i,3) f298 = ar1*(298.**ar2)*exp(-ar3/298.) * second reactant IF (rokrct(i,2).NE.' ') THEN r(2) = rokrct(i,2)(1:6) ENDIF ! brtio = brch ! force to highest yield ! first product brtio = brch * s(1) s(1) = rocost(i,1) CALL bratio(rokprd(i,1),brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) np=1 *second product IF (rokprd(i,2).NE.' ') THEN np=np+1 s(np) = rocost(i,2) brtio = brch * s(np) CALL bratio(rokprd(i,2),brtio,p(np), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ENDIF *third product IF (rokprd(i,3).NE.' ') THEN np=np+1 s(np) = rocost(i,3) brtio = brch * s(np) CALL bratio(rokprd(i,3),brtio,p(np), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ENDIF ! write reaction CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ENDIF ENDDO IF (known_species.ne.0) RETURN ********************************************************************** ! IF RINGS EXIST remove ring-join characters from groups, IF (nring.gt.0) THEN CALL rjgrm(nring,group,rjg) ENDIF tgroup = group tbond = bond ! locate alkoxy group DO i=1,mca IF (INDEX(group(i),alkoxy).NE.0) ia = i ENDDO ia0=ia IF (wtflag.NE.0) write(26,*) rdct ! identify if primary, secondary, tertiary: IF (group(ia)(1:3).EQ.primary) THEN itype = 1 ELSE IF (group(ia)(1:2).EQ.secondary) THEN itype = 2 ELSE IF (group(ia)(1:1).EQ.'C') THEN itype = 3 ELSE IF (group(ia)(1:1).EQ.'c') THEN itype = 4 ELSE WRITE(6,'(a)') '--warning, (stop) in ro' WRITE(6,'(a)') 'the following species has no type' WRITE(6,'(a)') '(e.g. primary, secondary, ...)' WRITE(6,'(a)') rdct(lco+1:lcf) WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF ********************************************************************** ! JMLT, Feb 2020: updated to include moieties other than NO2 ! If the alkoxy is borne by a carbon also bearing one of the following groups, ! we assume that there is only decomposition, which releases that group. IF (INDEX(tgroup(ia),'(NO2)').NE.0 .OR. & INDEX(tgroup(ia),'(ONO2)').NE.0 .OR. & INDEX(tgroup(ia),'(OOH)').NE.0 ) GOTO 111 ! If the alkoxy is borne by a carbon also bearing an (OOH) group, ! we assume that there is only decomposition, releasing HO2 ********************************************************************** ! Special case 1 : >C=CR-CR(polar)-C(O.)< ! assumed to only decompose to >C=CR-CR(polar)(.) + >C=O ! idem for >C=CR-CO-C(O.) structures DO i=1,mca IF (bond(i,ia).ne.0) THEN DO j=1,mca IF ((bond(i,j).eq.1).AND.(group(j)(1:2).EQ.'Cd')) THEN IF ((INDEX(tgroup(i),hydroxy).NE.0).OR. & (INDEX(tgroup(i),nitrate).NE.0).OR. & (INDEX(tgroup(i)(1:3),'CO ').NE.0).OR. & (INDEX(tgroup(i),hydro_peroxide).NE.0)) THEN specialcase = 1 GOTO 400 ENDIF ENDIF ENDDO ENDIF ENDDO ********************************************************************** ! Special case 2 : >C=CR-CO-O-CH(O.)- ! assumed to only be rearranged by the alpha-ester way DO i=1,mca IF (bond(i,ia).eq.3) THEN DO j=1,mca IF ((bond(i,j).eq.3).AND.(group(j)(1:2).EQ.'CO')) THEN DO k=1,mca IF ((bond(k,j).eq.1).AND.(group(k)(1:2).EQ.'Cd')) THEN specialcase = 2 ENDIF ENDDO ENDIF ENDDO ENDIF ENDDO ********************************************************************** ! PHENOXY RADICALS !! !!!!!!!!!!!!!!!!!!!!! ! For phenoxy radicals, three pathways are considered : ! 1 - reaction with O3 to form a peroxy radical ! 2 - reaction with NO2 to form a nitro phenol (if there is a H ! available close to the O. radical ! 3 - reaction with HO2 to form a phenolic product + O2 tgroup(:)=group(:) tbond(:,:)=bond(:,:) IF (itype.EQ.4) THEN IF (wtflag.NE.0) WRITE(6,*) 'phenoxy radical found: ', chem ! first pathway : ozone reaction IF (tgroup(ia).EQ.'c(O.)') THEN tgroup(ia)='c(OO.)' ELSE WRITE(6,*) 'error in ro.f - phenoxy + O3 ' STOP ENDIF CALL rebond(tbond,tgroup,tempkc,nring) CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) tempkc = rdckprod(1) IF (wtflag.NE.0)PRINT*,'radchk1' IF (nip.NE.1) STOP '2 produits phenoxy ro.f' CALL stdchm (tempkc) tgroup(:)=group(:) tbond(:,:)=bond(:,:) ! write reaction CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) a1 = rdct(1:1) r(1) = rdct(1:lco) r(2) = 'O3' s(1) = 1. ar1 = 2.86E-13 ar2 = 0 ar3 = 0 f298 = ar1*(298.**ar2)*exp(-ar3/298.) fratio=1. brtio=1 CALL bratio(tempkc,brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) WRITE(70,*) rdct(1:6),' ',ar1,ar2,ar3 & ,' ',chem(1:index(chem,' ')) ! second pathaway : nitro formation ! GOTO 963 tchem(:)=' ' j=0 DO i=1,nca IF ((tbond(i,ia).EQ.1).AND.(tgroup(i)(1:3).EQ.'cH ')) THEN j=j+1 tgroup(ia)='c(OH)' tgroup(i)='c(NO2)' CALL rebond(tbond,tgroup,tempkc,nring) CALL stdchm (tempkc) tchem(j)=tempkc tgroup(:)=group(:) tbond(:,:)=bond(:,:) ENDIF ENDDO IF (tchem(1).EQ.tchem(2)) tchem(2)=' ' ! write reactions: DO i=1,2 IF (tchem(i).NE.' ') THEN CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) a1 = rdct(1:1) r(1) = rdct(1:lco) r(2) = 'NO2' s(1) = 1. ar1 = 2.08E-12 ar2 = 0 ar3 = 0 f298 = ar1*(298.**ar2)*exp(-ar3/298.) fratio=1. brtio=1 CALL bratio(tchem(i),brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) WRITE(71,*) rdct(1:6),' ',ar1,ar2,ar3 & ,' ', chem(1:index(chem,' ')) ENDIF ENDDO !963 CONTINUE ! reaction with HO2 : phenoxy + HO2 -> phenolic products IF (tgroup(ia).EQ.'c(O.)') THEN tgroup(ia)='c(OH)' ELSE WRITE(6,*) 'error in ro.f - phenoxy + HO2 ' STOP ENDIF CALL rebond(tbond,tgroup,tempkc,nring) CALL stdchm (tempkc) tgroup(:)=group(:) tbond(:,:)=bond(:,:) ! write reaction CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) a1 = rdct(1:1) r(1) = rdct(1:lco) r(2) = 'HO2' s(1) = 1. ar1 = 1E-11 !personal communication, Mike Jenkin, MAGNIFY ar2 = 0 ar3 = 0 f298 = ar1*(298.**ar2)*exp(-ar3/298.) fratio=1. brtio=1 CALL bratio(tempkc,brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) RETURN ENDIF ********************************************************************** *---------------------------* * alpha-ester rearrangement * *---------------------------* * reaction rates are based on the SAPRC99 chemical scheme. * R-CO-O-CH(O.)R' => R-C(OH)=O + .CO-R' * Rate constants are estimated using k=Aexp(-Ea/RT), where * A=8.0E+10 and Ea/R=3723 in accordance with qualitative data of * Tuazon (98). * Only 1-4 shift is considered. 1-5 assumed to be not significant * NB : Because of lack of information, no such reaction for formate. ring_old=nring DO i=1,mca IF ((tbond(i,ia).EQ.3).AND.(tgroup(ia)(1:2).EQ.'CH')) THEN IF(nring.GT.0)THEN CALL findring(i,ia,nca,tbond,rngflg,ring) ENDIF DO j=1,mca IF ((tbond(j,i).NE.0).AND.(tgroup(j).EQ.'CO')) THEN flest = 1 pold = 'CO' pnew = 'CO(OH)' CALL swap (group(j),pold,tgroup(j),pnew) pold = '-O-' pnew = ' ' CALL swap (group(i),pold,tgroup(i),pnew) IF (itype.EQ.1) THEN pold = 'CH2(O.)' pnew = 'CHO.' ELSE IF (itype.EQ.2) THEN pold = 'CH(O.)' pnew = 'CO.' ENDIF CALL swap (group(ia),pold,tgroup(ia),pnew) tbond(i,j) = 0 tbond(j,i) = 0 tbond(ia,i) = 0 tbond(i,ia) = 0 IF (ring(ia).NE.0) THEN CALL rebond(tbond,tgroup,tempkc,nring) CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) pest(1) = rdckprod(1) IF (wtflag.NE.0)PRINT*,'radchk2' IF (nip.NE.1) WRITE(6,*) '2 produits a-ester ro.f' tprod(1:mca) = rdcktprod(1,1:mca) copest(1:mca) = tprod(1:mca) CALL stdchm(pest(1)) pest(2)=' ' ELSE * end modif IF(wtflag.NE.0)PRINT*,'fragm1' CALL fragm(tbond,tgroup,pest(1),pest(2)) DO k=1,2 IF (index(pest(k),'.').NE.0) THEN tempkc = pest(k) CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) pest(k) = rdckprod(1) IF (wtflag.NE.0)PRINT*,'radchk3' IF (nip.NE.1) WRITE(6,*) '2 produits a-ester ro.f' tprod(1:mca) = rdcktprod(1,1:mca) copest(1:mca) = tprod(1:mca) ENDIF CALL stdchm (pest(k)) ENDDO ENDIF ar1est = 8.E+10 ar2est = 0. ar3est = 3.723E+03 ! reset tgroup = group tbond = bond nring=ring_old ENDIF ENDDO ENDIF ENDDO IF (specialcase.EQ.2) GOTO 500 *************** Heat of formation of alkoxy radical ******** ! calculate heat of formation of alkoxy radical: dhfro = heat(rdct(lco+1:lcf),nbson,bsongrp,bsonval) ************************************************************** ! ------------------- ! reaction with O2 ! ------------------- * Reaction rates are based on the SAPRC99 chemical scheme. ! Rate constants are estimated using k=Aexp(-Ea/RT), where ! A=6.0E-14 and 1.5E-14 respectively for all primary & secondary RO ! and Ea=6.96+0.183DHr (in kcal.mol-1). Note that Ea is taken as ! the maximum : max(0.4;6.96+0.183DHr) ! primary RO : R-CH2(O.) ! UPDATE from Atkinson 2007 kO2=2.5E-14 exp(-300/T) for both ! primary and secondary alkoxy radicals (ric) IF (itype.EQ.1) THEN flo2 = 1 ar1o2 = 2.5E-14 tgroup(ia) = aldehyde ! secondary RO : R-CH(O.)-R ! disallow reaction for rings c-C5 or less (only decomposition is known) ELSE IF (itype.EQ.2) THEN IF(nring.GT.0)THEN DO i=1,mca IF(tbond(i,ia).NE.0) & CALL findring(i,ia,nca,tbond,rngflg,ring) ENDDO ncr = 0 DO i=1,mca IF(ring(i).NE.0)ncr=ncr+1 ENDDO IF(ncr.LE.5) THEN !flo2 = 0 !PRINT*,'O2 reaction disallowed for c-C5 and smaller' GOTO 100 ENDIF ENDIF !PRINT*,'O2 reaction allowed for c-C6 and larger' flo2 = 1 ar1o2 = 2.5E-14 pold = secondary pnew = carbonyl CALL swap(group(ia),pold,tempkg,pnew) pold = alkoxy pnew = ' ' CALL swap(tempkg,pold,tgroup(ia),pnew) ENDIF *rebuild and rename: IF (flo2.EQ.1) THEN CALL rebond(tbond,tgroup,po2,nring) CALL stdchm(po2) ENDIF ! note: carbonyl group may shift during re-writing of cyclic, but ! the groups are reset before ring decomposition so it doesnt matter ! compute the heat of reaction and activation energy ! dhfp1=heat(po2,nbson,bsongrp,bsonval) ! dhnet=(dhfp1+dhfho2)-(dhfro) ! ea=max(0.4,6.96 + 0.183*dhnet) ! ar3o2=(ea*1000.)/1.9872 ar3o2 = 300 ! If po2 is R-C(=O)ONO2 , assume decomposition => R. + CO2 + NO2 ! nb: doesnt involve breaking a ring bond IF(index(po2,'CO(ONO2)').ne.0) THEN DO i= 1,mca IF (index(tgroup(i),'CO(ONO2)').NE.0) THEN DO j=1,mca IF (bond(i,j).EQ.1) THEN tbond(i,j) = 0 tbond(j,i) = 0 nc = index(tgroup(j),' ') tgroup(j)(nc:nc) = '.' ! UPDATE MAY 2016 : following line commented, it shouldnt rebond before ! fragmenting ! CALL rebond(tbond,tgroup,tempkc,nring) IF(wtflag.NE.0)PRINT*,'fragm2' CALL fragm (tbond,tgroup,tempo2,tempkc) ! check radical from R. (not CO(ONO2) side) IF(INDEX(tempo2,'CO(ONO2)').NE.0) THEN CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) po2 = rdckprod(1) IF (wtflag.NE.0)PRINT*,'radchk4' IF (nip.NE.1) WRITE(6,*) '2 produits O2 ro.f' copo2(1:mca) = rdcktprod(1,1:mca) ELSE CALL radchk(tempo2,rdckprod,rdcktprod,nip,sc) po2 = rdckprod(1) IF (wtflag.NE.0)PRINT*,'radchk5' IF (nip.NE.1) WRITE(6,*) '2 produits O2 ro.f' copo2(1:mca) = rdcktprod(1,1:mca) ENDIF CALL stdchm(po2) DO k=1,mca IF (copo2(k)(1:1).EQ.' ')THEN copo2(k) = 'CO2' copo2(k+1) = 'NO2' GOTO 100 ENDIF IF (k.GE.mca) THEN mesg = 'k > mca (too many coproducts in the O2 reaction)' CALL errexit(progname,mesg,rdct(lco+1:lcf)) ENDIF ENDDO ENDIF ENDDO ENDIF ENDDO ENDIF 100 CONTINUE ! reset tgroup = group tbond = bond ! find nodes that are already part of any ring IF (nring.GT.0) THEN DO i=1,nring begrg=rjg(i,1) endrg=rjg(i,2) CALL findring(begrg,endrg,nca,tbond,rngflg,ring) ENDDO ENDIF ! ------------------- ! isomerization (for non-ring nodes only): ! ------------------- ! look down-chain from group JJ [if it is not already on a ring] ! first, five-membered rings : isomerisation rate for 5 membered ring ! are expected to be too slow to have a significant contribution and ! is not considered anymore below (fortran line must be comment lines DO 110 i=1,mca IF (bond(ia,i).EQ.0.OR.ring(i).GT.0) GO TO 110 oringa=0 IF (group(i)(1:3).EQ.ether) oringa=oringa+1 IF (group(i)(1:2).EQ.carbonyl) THEN oringa=oringa+1 DO l=1,mca IF (bond(i,l).EQ.3) oringa=oringa-1 ENDDO ENDIF IF (group(i)(1:3).EQ.aldehyde) oringa=oringa+1 DO 115 j=1,mca ! avoid step-back IF ( (bond(i,j).EQ.0) .OR. (j.EQ.ia) ) GO TO 115 oringb=0 IF (group(j)(1:3).EQ.ether) oringb=oringb+1 IF (group(j)(1:2).EQ.carbonyl) THEN oringb=oringb+1 DO l=1,mca IF (bond(j,l).EQ.3) oringb=oringb-1 ENDDO ENDIF IF (group(j)(1:3).EQ.aldehyde) oringb=oringb+1 ! look deeper for six-membered rings: DO 125 k=1,mca oringc=0 IF ( (bond(k,j).EQ.0) .OR. (k.eq.i) ) GO TO 125 IF (group(k)(1:3).EQ.ether) oringc=oringc+1 IF (group(k)(1:2).EQ.carbonyl) THEN oringc=oringc+1 DO l=1,mca IF (bond(k,l).EQ.3) oringc=oringc-1 ENDDO ENDIF IF (group(k)(1:3).EQ.aldehyde) oringc=oringc+1 IF (INDEX(group(k),'CH').NE.0) THEN niso = niso + 1 IF (niso.gt.mca) then WRITE(6,'(a)') '--error--, in subroutine ro' WRITE(6,'(a)') 'niso is greater than mca' WRITE(6,'(a)') 'need to modify size of the table' WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF fliso(niso) = 1 IF (kisomfg.EQ.2) THEN CALL isom_rate(tbond,tgroup,k,5,ar1,ar2,ar3,FSD_tmp) FSD(niso,:)=FSD_tmp(:) ELSE CALL rabsisom(tbond,tgroup,k,ar1,ar2,ar3) ENDIF ar1iso(niso)=ar1 ar2iso(niso)=ar2 ar3iso(niso)=ar3 ! add additional activation energy linked to ring strain oring=oringa+oringb+oringc ! IF (oring.gt.0) THEN ! ar3iso(niso)=ar3iso(niso)+3500/1.9872 ! ENDIF IF (group(k)(1:3).EQ.methyl) THEN pold = methyl pnew = primary ELSE IF (group(k)(1:3).EQ.primary) THEN pold = primary pnew = secondary ELSE IF (group(k)(1:3).EQ.aldehyde) THEN pold = aldehyde pnew = carbonyl ELSE IF (group(k)(1:2).EQ.secondary) THEN pold = secondary pnew = 'C' ELSE WRITE(6,'(a)') '--error--' WRITE(6,'(a)') 'from ROUTINE : rabsisom' WRITE(6,'(a)') 'group not allowed to react' WRITE(6,'(a)') 'by isomerization' WRITE(6,'(a)') tgroup(k) WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF CALL swap(group(k),pold,tgroup(k),pnew) nc = INDEX(tgroup(k),' ') tgroup(k)(nc:nc) = '.' pold = alkoxy pnew = hydroxy CALL swap(group(ia),pold,tgroup(ia),pnew) CALL rebond(tbond,tgroup,tempkc,nring) CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) piso(niso) = rdckprod(1) siso(niso,1) = sc(1) IF (wtflag.NE.0)PRINT*,'radchk6' IF (nip.EQ.2) THEN fliso2(niso) = 1 piso2(niso) = rdckprod(2) siso(niso,2) = sc(2) copiso2(niso,1:mca) = rdcktprod(2,1:mca) ENDIF copiso(niso,1:mca) = rdcktprod(1,1:mca) CALL stdchm(piso(niso)) ! reset tgroup(k) = group(k) tgroup(ia) = group(ia) rngflg=0 ENDIF 125 CONTINUE 115 CONTINUE 110 CONTINUE 111 CONTINUE ! ----------------------------- ! decomposition / ring opening ! ----------------------------- ! Special case for >C(O.)(NO2) : decomposes to release NO2 ! Special case for >C(O.)(ONO2): decomposes to release NO3 ! Special case for >C(O.)(OOH) : decomposes to release HO2 ! IF ((INDEX(tgroup(ia),'(O.)').NE.0).AND. ! & (INDEX(tgroup(ia),'(NO2)').NE.0)) THEN IF (INDEX(tgroup(ia),'(O.)').NE.0.AND. & (INDEX(tgroup(ia),'(NO2)' ).NE.0).OR. & (INDEX(tgroup(ia),'(ONO2)').NE.0).OR. & (INDEX(tgroup(ia),'(OOH)' ).NE.0))THEN IF (INDEX(tgroup(ia),'(NO2)' ).NE.0) icase = 1 IF (INDEX(tgroup(ia),'(ONO2)').NE.0) icase = 2 IF (INDEX(tgroup(ia),'(OOH)' ).NE.0) icase = 3 ndec = ndec + 1 IF (ndec.gt.mca) then WRITE(6,'(a)') '--error --, in subroutine ro' WRITE(6,'(a)') 'ndec is greater than mca' WRITE(6,'(a)') 'need to modify size of table' WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF fldec(ndec) = 1 pold='(O.)' pnew=' ' CALL swap(tgroup(ia),pold,tempkg,pnew) tgroup(ia)=tempkg IF (INDEX(tgroup(ia),'CH').NE.0) THEN pold='CH' pnew='CHO' ELSE pold='C' pnew='CO' ENDIF CALL swap(tgroup(ia),pold,tempkg,pnew) SELECT CASE (icase) CASE (1) ! (NO2) -> NO2 nc=INDEX(tempkg,'(NO2)') ilen = 5 pdec(ndec,2)='NO2' CASE (2) ! (ONO2) -> NO3 nc=INDEX(tempkg,'(ONO2)') ilen = 6 pdec(ndec,2)='NO3' CASE (3) ! (OOH) -> HO2 nc=INDEX(tempkg,'(OOH)') ilen = 5 pdec(ndec,2)='HO2' END SELECT tempkg(nc:)=tempkg(nc+ilen:) tgroup(ia)=tempkg CALL rebond(tbond,tgroup,pdec(ndec,1),nring) ar1dec(ndec)=1E12 ar2dec(ndec)=0 ar3dec(ndec)=0 GOTO 500 ENDIF ! find all C-C(O.) bonds. If the carbon is a Cd (C=C-C(O.)-), ! then assume no dissociation DO 150 ii=1,mca i=ii ia=ia0 IF (bond(i,ia).EQ.0) GO TO 150 IF (tgroup(i)(1:2).eq.'Cd') GOTO 150 ! post aromatic chemistry RIC 2016 magnify ! decomposition releasing phenyl radicals can be assumed not to occur IF (tgroup(i)(1:1).eq.'c') GOTO 150 ndec = ndec + 1 IF (ndec.gt.mca) then WRITE(6,'(a)') '--error --, in subroutine ro' WRITE(6,'(a)') 'ndec is greater than mca' WRITE(6,'(a)') 'need to modify size of table' WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF fldec(ndec) = 1 ! compute rate constant (Vereecken and Peeters, 2009) IF (wtflag.NE.0) PRINT*,"decomp_ro1" CALL decomp_ro(bond,tgroup,ia,i,Eb,nring) ! If the molecule is a C=C-C-CO. then place a message in the 22 file, ! to point out that thermochemical groups used to calculate the ! Ea are crude estimations. IF (wtflag.NE.0) THEN DO j=1,mca IF ((bond (i,j).eq.1).and.(group(j)(1:2).eq.'Cd')) THEN WRITE (22,*)'The thermochemical groups (bsongrpcarter.dat)' WRITE(22,*)'used to calculate Ea for beta decomposition' WRITE (22,*)'of RO. are crude estimations' WRITE (22,*) rdct(lco+1:lcf) ENDIF ENDDO ENDIF ! check if ring opening required IF(nring.GT.0)THEN IF (wtflag.NE.0) PRINT*,"nring = ",nring CALL findring(i,ia,nca,tbond,rngflg,ring) IF (wtflag.NE.0) PRINT*,"rngflg = ",rngflg ENDIF ! LINEAR DECOMPOSITION: IF (rngflg.EQ.0) THEN ! change (O.) to carbonyl: IF (itype.EQ.1) THEN tgroup(ia) = 'CH2O' ELSE IF (itype.EQ.2) THEN IF (wtflag.NE.0) print*,"group(ia) ",group(ia) pold = alkoxy pnew = ' ' CALL swap(group(ia),pold,tempkg,pnew) pold = secondary pnew = aldehyde CALL swap(tempkg,pold,tgroup(ia),pnew) ELSE IF (itype.EQ.3) THEN pold = alkoxy pnew = ' ' CALL swap(group(ia),pold,tempkg,pnew) pold = 'C' pnew = carbonyl CALL swap(tempkg,pold,tgroup(ia),pnew) IF (wtflag.NE.0) print*,"tgroup(ia) ",tgroup(ia) ENDIF ! break bond matrix, making two new fragments: tbond(i,ia) = 0 tbond(ia,i) = 0 ! linear decomposition: case of R-O-C(O.)< IF (tgroup(i)(1:3).EQ.'-O-') THEN DO j=1,mca IF ((tbond(i,j).EQ.3).AND.(j.NE.ia)) THEN nc = index(tgroup(j),' ') tgroup(j)(nc:nc+3) = '(O.)' tgroup(i)=' ' tbond(i,j) = 0 tbond(j,i) = 0 ! march 2015 if there were two ether groups, they must be converted to a peroxy group ! JMLT, Oct'15: case of R-O--O-C(O.)< : fragments to >C=O + R(OO.) IF (tgroup(j).EQ.'-O-(O.)') THEN DO k=1,nca IF ((tbond(j,k).NE.0).AND.(k.NE.i)) THEN tbond(k,j) = 0 tbond(j,k) = 0 tgroup(j) = ' ' nc = INDEX(tgroup(k),' ') tgroup(k)(nc:nc+4)='(OO.)' ENDIF ENDDO ENDIF ! CALL rebond(tbond,tgroup,tempkc,nring) IF(wtflag.NE.0)PRINT*,'fragm3' CALL fragm(tbond,tgroup,pdec(ndec,1),pdec(ndec,2)) ! JMLT, Oct'15: case of CHO-O--O-C(O.)< : fragments to >C=O + CO + HO2 IF(INDEX(pdec(ndec,1),'CHO(OO').NE.0) THEN pdec(ndec,1) = pdec(ndec,2) pdec(ndec,2) = ' ' copdec(ndec,1) = 'CO ' copdec(ndec,2) = 'HO2 ' sdec(ndec,1)=1. sdec(ndec,2)=1. ELSE IF(INDEX(pdec(ndec,2),'CHO(OO').NE.0) THEN pdec(ndec,2) = ' ' copdec(ndec,1) = 'CO ' copdec(ndec,2) = 'HO2 ' sdec(ndec,1)=1. sdec(ndec,2)=1. ENDIF ! JMLT, Nov'15: case of CO(OH)-O--O-C(O.)< : fragments to >C=O + CO2 + HO2 IF(INDEX(pdec(ndec,1),'CO(OH)(OO').NE.0) THEN pdec(ndec,1) = pdec(ndec,2) pdec(ndec,2) = ' ' copdec(ndec,1) = 'CO2 ' copdec(ndec,2) = 'HO2 ' sdec(ndec,1)=1. sdec(ndec,2)=1. ELSE IF(INDEX(pdec(ndec,2),'CO(OH)(OO').NE.0) THEN pdec(ndec,2) = ' ' copdec(ndec,1) = 'CO2 ' copdec(ndec,2) = 'HO2 ' sdec(ndec,1)=1. sdec(ndec,2)=1. ENDIF Eaa = 7.5 IF (wtflag.NE.0) WRITE(22,*)' ' IF (wtflag.NE.0) WRITE(22,*)'pdec(1)',pdec(ndec,1), & 'pdec(2)',pdec(ndec,2) dhfp1 = heat(pdec(ndec,1),nbson,bsongrp,bsonval) dhfp2 = heat(pdec(ndec,2),nbson,bsongrp,bsonval) dhnet = dhfp1 + dhfp2 - dhfro IF (wtflag.NE.0) WRITE(22,*)'dHnet=',dhnet GOTO 175 ENDIF ENDDO ELSE ! linear decomposition: other cases ! add dot to alkyl radical fragment nc = index(tgroup(i),' ') tgroup(i)(nc:nc) = '.' IF(wtflag.NE.0)PRINT*,'fragm4' CALL fragm(tbond,tgroup,pdec(ndec,1),pdec(ndec,2)) ! check if product is substituted alkene ! (and decompose it if necessary) DO k=1,2 IF (index(pdec(ndec,k),'.').EQ.0) THEN IF(INDEX(pdec(ndec,k),'Cd(O').NE.0 .OR. & INDEX(pdec(ndec,k),'CdH(O').NE.0) THEN CALL alkcheck(pdec(ndec,k),tprod(ii)) ENDIF ENDIF ENDDO ENDIF ELSE ! CYCLIC-TO-LINEAR RING OPENING: ! C1-C-C-C(O.)-C1 => C.-C-C-C-C=O ! [C(O.) was converted to carbonyl above] ! .. suppress opening if ONLY the OTHER alpha-carbon is substituted ! .. (ring-opening is preferred between alkoxy and substituted carbon) ! BA-RV : this test is not ok in the frame of terpene chemistry due ! to (C)(C)C(O.)-Cd= structure and lead to reaction without product ! and rate constant assigned. The selection should be made a the rate ! constant (must be lower for the decomposition on the unsubstituted ! carbon), which is the case for the acyclic parametrisation. This ! test is removed ! DO j=1,mca ! IF(j.NE.i.AND.tbond(j,ia).NE.0)THEN ! IF(INDEX(tgroup(j),methyl).EQ.0.AND. ! & INDEX(tgroup(j),primary).EQ.0.AND. ! & INDEX(tgroup(i),primary).NE.0) THEN ! !PRINT*,'suppress opening: ',ii,ia0 ! GOTO 112 ! ENDIF ! !PRINT*,'opening ring at: ',ii,ia0 ! ENDIF ! ENDDO ! CYCLIC ETHERS: SPECIAL CASE: ! John Orlando, pers. comm. 2020 ! IF group(i) is -O- AND group(ia) is >CH(O.) ! do not break ring, convert >CH(O.) to >C=O + HO2 (collision with O2) ! (Any other substituent on the >C(O.) DOES force ring-breaking: ! see later) IF (tgroup(i).EQ.'-O-' .AND. & INDEX(tgroup(ia),'H(O.)').GT.0) THEN ! .. change (O.) to carbonyl: pold = alkoxy pnew = ' ' CALL swap(tgroup(ia),pold,tempkg,pnew) pold = secondary pnew = carbonyl CALL swap(tempkg,pold,tgroup(ia),pnew) ! .. standardise formula of new molecule CALL rebond(tbond,tgroup,tempkc,nring) CALL stdchm(tempkc) nc = INDEX(tempkc,' ') - 1 pdec(ndec,1)=tempkc po2 = tempkc ! .. HO2 coproduct, fast reaction copdec(ndec,1)='HO2' flo2 = 1 ar1o2 = 2.5E-14 ar2o2 = 0. ar3o2 = 300. GO TO 176 ELSE ! .. break bond matrix, opening the ring: tbond(i,ia) = 0 tbond(ia,i) = 0 opflg=1 nring=nring-1 IF (wtflag.NE.0) PRINT*,"nring = ",nring ! CYCLIC ETHERS: OLD WAY: can still do this IF group(ia) is functionalised ! if group(i) is -O-, it must be converted to (O.) 2009 IF (tgroup(i).EQ.'-O-') THEN DO j=1,mca IF ((tbond(i,j).NE.0).AND.(j.NE.ia)) THEN tbond(i,j) = 0 tbond(j,i) = 0 tgroup(i) = ' ' nc = INDEX(tgroup(j),' ') tgroup(j)(nc:nc+3)='(O.)' ! if group(j) is -O-, we have -O--O- => convert to (OO.) !2015 IF (tgroup(j).EQ.'-O-(O.)') THEN DO k=1,nca IF ((tbond(j,k).NE.0).AND.(k.NE.i)) THEN tbond(k,j) = 0 tbond(j,k) = 0 tgroup(j) = ' ' nc = INDEX(tgroup(k),' ') tgroup(k)(nc:nc+4)='(OO.)' ENDIF ENDDO ENDIF fldecether = 1 ENDIF ENDDO ELSE nc = INDEX(tgroup(i),' ') tgroup(i)(nc:nc) = '.' ENDIF ! .. change (O.) to carbonyl: IF (itype.EQ.2) THEN pold = alkoxy pnew = ' ' CALL swap(tgroup(ia),pold,tempkg,pnew) pold = secondary pnew = aldehyde CALL swap(tempkg,pold,tgroup(ia),pnew) ELSE IF (itype.EQ.3) THEN pold = alkoxy pnew = ' ' CALL swap(tgroup(ia),pold,tempkg,pnew) pold = 'C' pnew = carbonyl CALL swap(tempkg,pold,tgroup(ia),pnew) ENDIF ! .. find artificial break point DO j=1,nca IF(ring(j).NE.0) endrg=j ENDDO DO j=nca,1,-1 IF(ring(j).NE.0) begrg=j ENDDO ! .. find if ring opens at artificial break point ! IF((i.EQ.begrg.AND.ia.EQ.endrg) ! $ .OR.(ia.EQ.begrg.AND.i.EQ.endrg))THEN ! CONTINUE ! ELSE ! .. if opens elsewhere, rearrange new linear molecule ! CALL rejoin(nca,i,ia,m,n,tbond,tgroup) ! i = m ! ia = n ! ENDIF ! .. standardise formula of new linear molecule CALL rebond(tbond,tgroup,tempkc,nring) CALL stdchm(tempkc) IF (wtflag.NE.0) print*,"tempkc ",tempkc nc = INDEX(tempkc,' ') - 1 pdec(ndec,1)=tempkc ! radchk gets called later (a couple of screens-worth down) ENDIF ENDIF ! end special case for cyclic ether with >CH(O.) ! compute rate constants: ! rate constants are taken from the SAPRC99 chemical scheme. ! Rate is : k=2E14exp(-Ea/RT), where Ea is given as : ! Ea=Eaa+0.44DHr=-8.73+2.35IP+0.44DHr with ! DHr : heat of the reaction (in kcal) ! IP : ionization potential (in eV) ! In what follows, values for Eaa (in kcal) are directly provided and ! values are taken from the SAPRC99 documentation ! UPDATE Eea taken from Atkinson 2007 (ric) ! compute heat of formation of each product IF (wtflag.NE.0) WRITE(22,*)' ' IF (wtflag.NE.0) WRITE(22,*)'pdec(1)',pdec(ndec,1),'pdec(2)' & ,pdec(ndec,2) dhfp1 = heat(pdec(ndec,1),nbson,bsongrp,bsonval) dhfp2 = heat(pdec(ndec,2),nbson,bsongrp,bsonval) dhnet = dhfp1 + dhfp2 - dhfro IF (wtflag.NE.0) WRITE(22,*)'dHnet=',dhnet ! set the Eaa factor IF (tgroup(i)(1:4).EQ.'CH3.') THEN Eaa=12.9 ELSE IF(tgroup(i)(1:3).EQ.'CH2') THEN IF (INDEX(tgroup(i),hydroxy).NE.0) THEN Eaa=6.8 ELSE Eaa=9.5 ENDIF ELSE IF(tgroup(i)(1:2).EQ.'CH') THEN IF (tgroup(i)(1:3).EQ.'CHO') THEN Eaa=9.99 ELSE IF (INDEX(tgroup(i),hydroxy).NE.0) THEN Eaa = 5.2 ELSE IF (INDEX(tgroup(i),nitrate).NE.0) THEN Eaa = 7.46 ELSE IF (INDEX(tgroup(i),hydro_peroxide).NE.0) THEN Eaa = 7.46 ELSE Eaa = 8.2 ENDIF ELSE IF(tgroup(i)(1:2).EQ.'CO') THEN Eaa = 5. ELSE IF(tgroup(i)(1:1).EQ.'C') THEN IF (INDEX(tgroup(i),hydroxy).NE.0) THEN Eaa = 4.8 ELSE IF (INDEX(tgroup(i),nitrate).NE.0) THEN Eaa = 5.58 ELSE IF (INDEX(tgroup(i),hydro_peroxide).NE.0) THEN Eaa = 5.58 ELSE Eaa = 7.1 ENDIF ELSE IF (fldecether.EQ.1) THEN Eaa = 7.5 ELSE IF((tgroup(i)(1:1).EQ.' ').AND. & (copdec(ndec,1).EQ.'CO2')) THEN Eaa = 5. ELSE WRITE(6,'(a)') 'tgroup(i)=',tgroup(i) WRITE(6,'(a)') 'copdec=',copdec(ndec,1:5) WRITE(6,'(a)') '--error--, in ro.f' WRITE(6,'(a)') 'Eaa factor can not be attributed' WRITE(6,'(a)') 'when treating the decomposition reaction' WRITE(6,'(a)') 'for the species :' WRITE(6,'(a)') rdct(lco+1:lcf) WRITE(99,*) 'ro',rdct(lco+1:lcf) STOP ENDIF 175 CONTINUE IF (wtflag.NE.0) WRITE(22,*) 'Eaa=',eaa IF (kdissfg.EQ.2) THEN ar1dec(ndec)=1.12E09 ar2dec(ndec)=1.7 ar3dec(ndec)=Eb ELSE ar1dec(ndec)=5.0E13 ar2dec(ndec)=0. ea=(Eaa+0.40*dhnet) ar3dec(ndec)=max(0.75,ea) IF (wtflag.NE.0) WRITE(22,*)'ar3dec(ndec)',ar3dec(ndec) ar3dec(ndec)=(ar3dec(ndec)*1000.)/1.9872 ENDIF 176 CONTINUE IF (wtflag.NE.0) WRITE(22,*) 'ar1dec=',ar1dec(ndec) IF (wtflag.NE.0) WRITE(22,*) 'ar3dec=',ar3dec(ndec) ! identify radical fragment (KK), run through radical checker routine IF (INDEX(pdec(ndec,1),'.').NE.0) THEN irad = 1 imol = 2 ELSE irad = 2 imol = 1 ENDIF IF(INDEX(pdec(ndec,irad)," ").GT.1)THEN tempkc = pdec(ndec,irad) IF (wtflag.NE.0)PRINT*,'radchk7' CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) pdec(ndec,irad) = rdckprod(1) sdec(ndec,1) = sc(1) IF (nip.EQ.2) THEN fldec2(ndec) = 1 pdec2(ndec) = rdckprod(2) sdec(ndec,2) = sc(2) DO j=1,mca copdec2(ndec,j) = rdcktprod(2,j) ENDDO ENDIF DO j=1,mca copdec(ndec,j) = rdcktprod(1,j) ENDDO CALL stdchm(pdec(ndec,irad)) ENDIF ! non-radical fragment (IMOL): IF (pdec(ndec,imol).NE.'CHO(ONO2)') & CALL stdchm(pdec(ndec,imol)) * If pdec is R-C(=O)ONO2 , we assume decomposition to R. + CO2 + NO2 ! Ludo: This part of the decomposition causes problems because the former ! value of tempkc is still considered and the final products are not ! the good ones. IF(index(pdec(ndec,imol),'CO(ONO2)').NE.0) THEN ! DO j=1,mca ! IF (index(tgroup(j),'.').NE.0) i=j ! ENDDO nc=index(pdec(ndec,imol),' ')-1 CALL grbond(pdec(ndec,imol),nc,grdec,bddec,dbf,nrg) DO k=1,mca ! IF (index(tgroup(k),'CO(ONO2)').NE.0) THEN IF (index(grdec(k),'CO(ONO2)').NE.0) THEN DO l=1,mca IF (bddec(l,k).EQ.1) THEN ! IF (tbond(l,k).EQ.1) THEN !IF (bond(l,k).EQ.1) THEN bddec(l,k) = 0 bddec(k,l) = 0 nc = index(grdec(l),' ') grdec(l)(nc:nc) = '.' ! tbond(l,k) = 0 ! tbond(k,l) = 0 ! nc = index(tgroup(l),' ') ! tgroup(l)(nc:nc) = '.' ! make sure non-radical fragment is listed first in disconnected bond matrix ! IF(k.GT.i) CALL rejoin(nca,i,k,m,n,tbond,tgroup) ! CALL rebond(tbond,tgroup,tempkc,nring) ! CALL fragm (tbond,tgroup,tempdec,tempkc) IF(wtflag.NE.0)PRINT*,'fragm5' CALL fragm(bddec,grdec,tempdec,tempkc) IF(INDEX(tempdec,'CO(ONO2)').NE.0) THEN ! CALL radchk(tempkc,pdec(ndec,k),tprod) CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) pdec(ndec,imol) = rdckprod(1) IF (wtflag.NE.0)PRINT*,'radchk8' IF (nip.NE.1) WRITE(6,*) '2 produits DEC 2 ro.f' DO m = 1,mca tprod(m) = rdcktprod(1,m) ENDDO ELSE ! CALL radchk(tempdec,pdec(ndec,k),tprod) CALL radchk(tempdec,rdckprod,rdcktprod,nip,sc) pdec(ndec,imol) = rdckprod(1) IF (wtflag.NE.0)PRINT*,'radchk9' IF (nip.NE.1) WRITE(6,*) '2 produits DEC 3 ro.f' DO m = 1,mca tprod(m) = rdcktprod(1,m) ENDDO ENDIF CALL stdchm(pdec(ndec,imol)) DO 113 j=1,mca IF (tprod(j).NE.' ') THEN DO n=1,mca IF (copdec(ndec,n).EQ.' ') THEN copdec(ndec,n) = tprod(j) GOTO 113 ENDIF ENDDO ENDIF 113 CONTINUE DO m=1,mca IF (copdec(ndec,m)(1:1).EQ.' ')THEN copdec(ndec,m) = 'CO2' copdec(ndec,m+1) = 'NO2' GOTO 112 ENDIF IF (m.GE.mca) THEN WRITE(6,*)'From ro.f, too many coprods ' WRITE(6,*) 'in decomposition reaction' WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF ENDDO ENDIF ENDDO ENDIF ENDDO ELSE IF (index(pdec(ndec,imol),'CHO(ONO2)').ne.0) THEN pdec(ndec,imol) = 'CO2' DO m=1,mca IF (copdec(k,m)(1:1).EQ.' ')THEN copdec(k,m) = 'NO2' copdec(k,m+1) = 'HO2' GOTO 112 ENDIF ENDDO ENDIF 112 CONTINUE ! reset: IF (opflg.EQ.1) THEN nring=nring+1 rngflg=0 opflg=0 ENDIF ia = ia0 tgroup = group tbond = bond nrg=0 dbf=0 grdec=' ' bddec=0 150 CONTINUE * ------------------------------ * check for duplicate reactions: * ------------------------------ ! isomerisation DO 200 i=1,niso-1 DO 205 j=i+1,niso IF (piso(j).EQ.piso(i)) THEN ! ric modif to avoid duplicate reaction in case of CH3C(CH3)(CH3)C(CH3)(CH3)CH3 IF (fliso(j).EQ.0) GOTO 205 ! end modif DO k=1,mca IF (copiso(i,k).NE.copiso(j,k)) GOTO 205 ENDDO IF (ar1iso(i).ne.ar1iso(j)) GOTO 205 IF (ar3iso(i).ne.ar3iso(j)) GOTO 205 fliso(j) = 0 fliso(i) = fliso(i)+1 ENDIF 205 CONTINUE 200 CONTINUE DO i=1,niso IF (fliso(i).gt.1) THEN ar1iso(i)=REAL(fliso(i))*ar1iso(i) fliso(i)=1 ENDIF ENDDO ! decomposition ! keep only 5 numbers after comma before compare ar3dec DO i=1,ndec resu=' ' write(resu,'(e12.5)') ar3dec(i) read(resu,'(e12.5)') ar3dec(i) ENDDO DO 210 i=1,ndec-1 DO 215 j=i+1,ndec IF ( (pdec(j,1).EQ.pdec(i,1) .AND. pdec(j,2).EQ.pdec(i,2)) | .OR. | (pdec(j,1).EQ.pdec(i,2) .AND. pdec(j,2).EQ.pdec(i,1)) | ) THEN DO k=1,mca IF (copdec(i,k).NE.copdec(j,k)) GO TO 215 ENDDO IF (ar1dec(i).ne.ar1dec(j)) GOTO 215 IF (ar3dec(i).ne.ar3dec(j)) GOTO 215 IF (fldec(i).eq.0) GOTO 215 fldec(j) = 0 fldec(i) = fldec(i)+1 ENDIF 215 CONTINUE 210 CONTINUE DO i=1,ndec IF (fldec(i).gt.1) THEN ar1dec(i)=REAL(fldec(i))*ar1dec(i) fldec(i)=1 ENDIF ENDDO ******************************************************** * Decomposition for : Cd-Ci(polar)-Cia-O. structures * ******************************************************** 400 CONTINUE IF (specialcase.eq.1) THEN ndec = 1 fldec(ndec) = 1 ! change (O.) to carbonyl: IF (itype.EQ.1) THEN tgroup(ia) = 'CH2O' ELSE IF (itype.EQ.2) THEN pold = alkoxy pnew = ' ' CALL swap(group(ia),pold,tempkg,pnew) pold = secondary pnew = aldehyde CALL swap(tempkg,pold,tgroup(ia),pnew) ELSE IF (itype.EQ.3) THEN pold = alkoxy pnew = ' ' CALL swap(group(ia),pold,tempkg,pnew) pold = 'C' pnew = carbonyl CALL swap(tempkg,pold,tgroup(ia),pnew) ENDIF ! check if bond that will be broken belongs to a ring. CALL findring(i,ia,nca,tbond,rngflg,ring) ! break bond matrix, make two new fragments: tbond(i,ia) = 0 tbond(ia,i) = 0 ! add dot to alkyl radical fragment: nc = index(group(i),' ') tgroup(i)(nc:nc) = '.' IF(rngflg.EQ.1)THEN CALL rebond(tbond,tgroup,tempkc,nring) ENDIF IF (rngflg.EQ.0) THEN IF(wtflag.NE.0)PRINT*,"fragm6" CALL fragm(tbond,tgroup,pdec(ndec,1),pdec(ndec,2)) ELSE pdec(ndec,1)=tempkc pdec(ndec,2)=' ' ENDIF ! rate constant is assumed to be k=2E14exp(-Ea/RT), ! where Ea = 0.75 ! RV Vereecken SAR added IF (kdissfg.EQ.2) THEN CALL decomp_ro(bond,tgroup,ia,i,Eb,nring) IF (wtflag.NE.0) PRINT*,"decomp_ro2" ar1dec(ndec)=1.12E09 ar2dec(ndec)=1.7 ar3dec(ndec)=Eb ELSE ar1dec(ndec)=2.0E14 ar2dec(ndec)=0. ar3dec(ndec)=0.75 ar3dec(ndec)=(ar3dec(ndec)*1000.)/1.9872 ENDIF ! identify radical fragment (KK), run through radical checker routine IF (INDEX(pdec(ndec,1),'.').NE.0) THEN irad = 1 imol = 2 ELSE irad = 2 imol = 1 ENDIF tempkc = pdec(ndec,irad) CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) IF(wtflag.NE.0)PRINT*,'radchk1' pdec(ndec,irad) = rdckprod(1) sdec(ndec,1) = sc(1) IF (nip.EQ.2) THEN fldec2(ndec) = 1 pdec2(ndec) = rdckprod(2) sdec(ndec,2) = sc(2) DO j=1,mca copdec2(ndec,j) = rdcktprod(2,j) ENDDO ENDIF DO l = 1,mca tprod(l) = rdcktprod(1,l) ENDDO CALL stdchm(pdec(ndec,irad)) DO j=1,mca copdec(ndec,j) = tprod(j) ENDDO ! non-radical fragment (IMOL): CALL stdchm(pdec(ndec,imol)) ! reset: tbond(i,ia) = bond(i,ia) tbond(ia,i) = bond(ia,i) tgroup(i) = group(i) ENDIF ! ------------------------------------------------ ! compare rate constants, keep only if 2% or more. ! ------------------------------------------------ 500 CONTINUE rsum=0. IF (flest.EQ.1) THEN kest = ar1est*(temp**ar2est)*exp(-ar3est/temp) rsum = rsum + kest ENDIF sum_o2 = ar1o2*(temp**ar2o2)*exp(-ar3o2/temp)*tbm*0.209 rsum = rsum + sum_o2 DO i=1,niso IF (fliso(i).EQ.1) THEN kiso(i) = ar1iso(i)*(temp**ar2iso(i))*exp(-ar3iso(i)/temp) IF (kisomfg.EQ.2) THEN kiso(i)=kiso(i)*(FSD(i,1)*temp**4 + FSD(i,2)*temp**3 + & FSD(i,3)*temp**2 + FSD(i,4)*temp + FSD(i,5)) ENDIF rsum = rsum + kiso(i) ENDIF ENDDO DO i=1,ndec IF (fldec(i).EQ.1) THEN kdec(i) = ar1dec(i)*(temp**ar2dec(i))*exp(-ar3dec(i)/temp) rsum = rsum + kdec(i) ENDIF ENDDO rsum = 0.02*rsum IF (kest.LT.rsum) flest = 0 IF (sum_o2.LT.rsum) flo2 = 0 DO i=1,niso IF (kiso(i).LT.rsum) fliso(i) = 0 ENDDO DO i=1,ndec IF (kdec(i).LT.rsum) fldec(i) = 0 ENDDO ! compute sum of the rate constant at temp and at 298 K rsum=0. rsum298=0. nich = 0 IF (flest.eq.1) THEN nich = nich+1 rsum = kest rsum298 = ar1est*(298.**ar2est)*exp(-ar3est/298.) ENDIF IF (flo2.eq.1) THEN nich=nich+1 rsum= rsum + sum_o2 rsum298 = rsum298 + ar1o2*(298.**ar2o2)*exp(-ar3o2/298.) & *2.45E19*0.209 ENDIF DO i=1,mca IF (fliso(i).EQ.1) THEN nich = nich + 1 rsum = rsum + kiso(i) rsum298 = rsum298 + & ar1iso(i)*(298.**ar2iso(i))*exp(-ar3iso(i)/298.) ENDIF IF (fldec(i).EQ.1) THEN nich = nich + 1 rsum = rsum + kdec(i) rsum298 = rsum298 + & ar1dec(i)*(298.**ar2dec(i))*exp(-ar3dec(i)/298.) ENDIF ENDDO IF (wtopeflag.EQ.1) & write(10,*)'***********! alkoxy ','G',rdct(1:6),'********' ! write out alpha ester rearrangement: ! ------------------------------------- ich = 10 IF (flest.EQ.1) THEN CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ich = ich + 1 IF (nich.GT.1) a4 = alfa(ich:ich) ! brtio = brch*kest/rsum brtio = brch ! first product s(1) = 1. CALL bratio(pest(1),brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) *second product s(2) = 1. CALL bratio(pest(2),brtio,p(2), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) IF (wtopeflag.EQ.1) & write(10,21)'ester rearrgt /',brch,'G',p(1),' + ','G',p(2) 21 FORMAT(A15,5X,f5.3,2X,A1,A6,A3,A1,A6) ! other products np = 2 DO j=1,mca IF (copest(j)(1:1).NE.' ') THEN np = np + 1 IF (np.GT.mnp) then WRITE(6,'(a)') '--error--1 in ro' WRITE(6,'(a)') 'np is greater than mnp' WRITE(6,'(a)') '(too many products in the reaction)' WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF s(np) = 1. p(np) = copest(j) IF (wtopeflag.EQ.1) write(10,'(20X,A1,A6)') 'G',p(np) ENDIF ENDDO IF (rdtcopchem1.GT.0.) THEN np = np + 1 IF (np.GT.mnp) then WRITE(6,'(a)') '--error-- in ro' WRITE(6,'(a)') 'np is greater than mnp' WRITE(6,'(a)') '(too many products in the reaction)' WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF s(np) = rdtcopchem1 p(np) = copchem1 IF (wtopeflag.EQ.1) & WRITE(10,'(20X,A1,f2.0,A1,A6)')'#',s(np),'G',p(np) ENDIF IF (rdtcopchem2.GT.0.) THEN np = np + 1 IF (np.GT.mnp) then WRITE(6,'(a)') '--error-- in ro' WRITE(6,'(a)') 'np is greater than mnp' WRITE(6,'(a)') '(too many products in the reaction)' WRITE(99,*) 'ro',rdct(lco+1:lcf) !STOP ENDIF s(np) = rdtcopchem2 p(np) = copchem2 IF (wtopeflag.EQ.1) & WRITE(10,'(20X,A1,f2.0,A1,A6)')'#',s(np),'G',p(np) ENDIF ! reactant a1 = '1' r(1) = rdct(1:lco) ar1 = ar1est ar2 = ar2est ar3 = ar3est f298 = ar1*(298.**ar2)*exp(-ar3/298.) fratio = f298/rsum298 CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ENDIF ! write out O2 reaction: ! ----------------------- IF (flo2.EQ.1) THEN CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ich = ich + 1 IF (nich.GT.1) a4 = alfa(ich:ich) s(1) = 1. ! brtio = brch*sum_o2/rsum brtio = brch CALL bratio(po2,brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) IF (wtopeflag.EQ.1) write(10,22)'O2 -> HO2 /',sum_o2/rsum,'G',p(1) 22 FORMAT(A14,4X,f5.3,2X,A1,A6) ! second product: HO2 s(2) = 1. p(2) = 'HO2 ' ! other products np = 2 DO j=1,mca IF (copo2(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,rdct(lco+1:lcf)) ENDIF s(np) = 1. p(np) = copo2(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: 2)' CALL errexit(progname,mesg,rdct(lco+1:lcf)) ENDIF s(np) = rdtcopchem p(np) = copchem IF (wtopeflag.EQ.1) & WRITE(10,'(23X,A1,f5.3,1X,A1,A6)')'#',sum_o2*s(np)/rsum, & 'G',p(np) ENDIF ! reactant - use label 300 for the reaction with O2 a1 = '1' r(1) = rdct(1:lco) ! r(2) = 'O2 ' ! r(2) = 'EXTRA ' r(2) = 'OXYGEN' ! nlabel=300 ar1 = ar1o2 ar2 = ar2o2 ar3 = ar3o2 f298 = ar1*(298.**ar2)*exp(-ar3/298.)*2.45E19*0.209 fratio = f298/rsum298 ! write out - extra reaction => idreac=2 (nlabel was set above) ! idreac=2 idreac=0 CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ENDIF ! write out isomerization ! ----------------------- DO 300 i=1,niso IF (fliso(i).EQ.0) GO TO 300 CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ich = ich + 1 IF (nich.GT.1) a4 = alfa(ich:ich) s(1) = 1. ! brtio = brch*kiso(i)/rsum brtio = brch CALL bratio(piso(i),brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) IF (wtopeflag.EQ.1) & write(10,22)'isomerisa° /',kiso(i)/rsum,'G',p(1) !IF (wtflag.EQ.1) write(6,*) '*bratio* : output :',p(1),' ',piso(i) !IF (wtflag.EQ.1) write(6,*) p(1),' ',piso(i) ! second product (if delocalisation allowed) np = 1 IF (fliso2(i).NE.0) THEN np = np + 1 s(1) = siso(i,1) s(np) = siso(i,2) CALL bratio(piso2(i),brtio,p(np), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ENDIF ! other products ! np = 1 DO j=1,mca IF (copiso(i,j)(1:1).NE.' ') THEN np = np + 1 IF (np.GT.mnp) then mesg = 'np > mnp (too many products in the reaction: 3)' CALL errexit(progname,mesg,rdct(lco+1:lcf)) ENDIF s(np) = 1. p(np) = copiso(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: 4)' CALL errexit(progname,mesg,rdct(lco+1:lcf)) ENDIF s(np) = rdtcopchem p(np) = copchem IF (wtopeflag.EQ.1) & WRITE(10,'(20X,A1,f5.3,1X,A1,A6)')'#',kiso(i)*s(np)/rsum, & 'G',p(np) ENDIF ! reactant a1 = '1' r(1) = rdct(1:lco) IF (kisomfg.EQ.2) THEN r(2) = 'ISOM ' ENDIF ar1 = ar1iso(i) ar2 = ar2iso(i) ar3 = ar3iso(i) f298 = ar1*(298.**ar2)*exp(-ar3/298.) fratio = f298/rsum298 CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) IF (kisomfg.EQ.2) THEN WRITE(17,'(A7,5(ES10.3,1x),A1)') ' ISOM/',(FSD(i,j),j=1,5),'/' ENDIF 300 CONTINUE ! write out decomposition ! ----------------------- DO 350 i=1,ndec IF (fldec(i).EQ.0) GO TO 350 CALL rxinit3(a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ich = ich + 1 IF (nich.GT.1) a4 = alfa(ich:ich) ! first product s(1) = 1. ! brtio = brch*kdec(i)/rsum brtio = brch CALL bratio(pdec(i,1),brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) np = 1 ! second product IF(INDEX(pdec(i,2),' ').GT.1)THEN np = np+1 s(np) = 1. CALL bratio(pdec(i,2),brtio,p(np), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) IF (wtopeflag.EQ.1) & write(10,21)'decomp / ',kdec(i)/rsum,'G',p(1),' + ','G',p(2) ELSE IF (wtopeflag.EQ.1) & write(10,21)'decomp / ',kdec(i)/rsum,'G',p(1) ENDIF ! third product (if delocalisation allowed) IF (fldec2(i).NE.0) THEN np = np+1 s(1) = sdec(i,1) s(np) = sdec(i,2) CALL bratio(pdec2(i),brtio,p(np), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ENDIF ! other coproducts ! np = 2 DO j=1,mca IF (copdec(i,j)(1:1).NE.' ') THEN np = np + 1 IF (np.GT.mnp) then mesg = 'np > mnp (too many products in the reaction: 5)' CALL errexit(progname,mesg,rdct(lco+1:lcf)) ENDIF s(np) = 1. p(np) = copdec(i,j) IF (wtopeflag.EQ.1) write(10,'(20X,A1,A6)') 'G',p(np) ENDIF ENDDO DO j=1,mca IF ((copdec2(i,j)(1:1).NE.' ').AND.(fldec2(i).NE.0)) THEN np = np + 1 IF (np.GT.mnp) then mesg = 'np > mnp (too many products in the reaction: 6)' CALL errexit(progname,mesg,rdct(lco+1:lcf)) ENDIF s(np) = 1. p(np) = copdec2(i,j) IF (wtopeflag.EQ.1) write(10,'(20X,A1,A6)') 'G',p(np) ENDIF ENDDO IF (rdtcopchem1.GT.0.) THEN np = np + 1 IF (np.GT.mnp) then mesg = 'np > mnp (too many products in the reaction: 7)' CALL errexit(progname,mesg,rdct(lco+1:lcf)) ENDIF s(np) = rdtcopchem1 p(np) = copchem1 IF (wtopeflag.EQ.1) & WRITE(10,'(23X,A1,f5.3,1X,A1,A6)')'#',s(np)*kdec(i)/rsum, & 'G',p(np) ENDIF IF (rdtcopchem2.GT.0.) THEN np = np + 1 IF (np.GT.mnp) then mesg = 'np > mnp (too many products in the reaction: 8)' CALL errexit(progname,mesg,rdct(lco+1:lcf)) ENDIF s(np) = rdtcopchem2 p(np) = copchem2 IF (wtopeflag.EQ.1) & WRITE(10,'(23X,A1,f5.3,1X,A1,A6)')'#',s(np)*kdec(i)/rsum, & 'G',p(np) ENDIF ! reactant a1 = '1' r(1) = rdct(1:lco) ar1 = ar1dec(i) ar2 = ar2dec(i) ar3 = ar3dec(i) f298 = ar1*(298.**ar2)*exp(-ar3/298.) fratio = f298/rsum298 CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) 350 CONTINUE IF (wtflag.EQ.1) write(6,*)'*end of ro*' IF (wtopeflag.EQ.1) write(6,*)'end of ro' IF (wtflag.EQ.1) write(10,*)'end' ! DEBUG ! ! PRINT*,"!!debug end ro!!" ! STOP ! DEBUG ! RETURN END