************************************************************************ * MASTER MECHANISM - ROUTINE NAME : rozone2 * * * * * * PURPOSE : * * perform reaction of O3 with VOC, i.e. : * * --> Add O3 to unsaturated compounds, make carbonyl and criegee * * * * NOTE : hot criegee (and as well as thermalized "stable" criegee) * * are not included as a species in the mechanism and are * * substituted by their decomposition or reaction products. The * * substitution is achieved by the xcrieg subroutine, which is * * called at the end of the subroutine, just before writing the * * product. * * * * 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 rings present * * - 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 O3 (see o3_rate.dat) * * arrhenius parameter (see o3_rate.dat) * * - no3dat : total number of known NO3 rate constant * * - o3chem(i) : chemical formula of ith species having a known NO3 * * rate constant * * - o3rate(i,3) : arrhenius parameter of the known NO3 rate constant * * corresponding to the ith species * * * * 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 * * * ************************************************************************ SUBROUTINE rozone2(chem,rdct,bond,group,nring,brch, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & no3dat,o3chem,o3rate,nfn,namfn,chemfn, & nkwo3,nkwo3pd,kwo3rct,kwo3pd,kwo3copd,kwo3yld) IMPLICIT NONE INCLUDE 'general.h' INCLUDE 'organic.h' INCLUDE 'common.h' * input: CHARACTER(LEN=lcf) rdct CHARACTER(LEN=lgr) group(mca) INTEGER bond(mca,mca) REAL brch REAL cut_off CHARACTER(LEN=lfo) o3chem(mrd) REAL o3rate(mrd,3) INTEGER no3dat INTEGER nkwo3,nkwo3pd(mkr) CHARACTER(LEN=lfo) kwo3rct(mkr), kwo3pd(mkr,mnr) CHARACTER(LEN=lco) kwo3copd(mkr,mnr,mcp) REAL kwo3yld(mkr,mnr) * 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: allow for up to "mnr" reaction channels CHARACTER(LEN=lfo) chem,tempkc CHARACTER(LEN=lfo) pchem(mnr,2) CHARACTER(LEN=lgr) tgroup(mca), pold, pnew, tempkg CHARACTER(LEN=lco) pp(mnr),coprod(mnr,mca) REAL ss(mnr) INTEGER tbond(mca,mca), flag(mnr) INTEGER nr,np,nch,ich,i,j,k,l,j1,kk,ii,nc,ncd REAL brtio, totr, sum, rate(mnr) REAL rtot REAL ratio(mnr) REAL tarrhc(mnr,3),rarrhc(3) REAL rk298 INTEGER cdtable(4),tcdtable(4), cdsub(4),conjug, ncdcase INTEGER cdcarbo(4,2),cdeth(4,2) INTEGER cpcdtable(4) CHARACTER*1 a1,a2,a3,a4 CHARACTER(LEN=lco) r(3), p(mnp) REAL s(mnp), ar1,ar2,ar3,f298,fratio INTEGER idreac, nlabel REAL xlabel,folow(3),fotroe(4) CHARACTER(LEN=lco) copchem REAL rdtcopchem,A,E,nt INTEGER nring,ring(mca),rjg(mri,2) INTEGER rxnflg CHARACTER(LEN=lfo) rdckprod(mca) CHARACTER(LEN=lco) rdcktprod(mca,mca) INTEGER nip REAL sc(mca) INTEGER known_species,no3prod INTEGER :: parent_group(mnr) REAL :: arrhc(mnr,3) REAL :: totr1 CHARACTER(LEN=lsb) :: progname='*rozone2* ' CHARACTER(LEN=ler) :: mesg * write info for finding bugs IF (wtflag.NE.0) WRITE(*,*) progname,'----------: input : ', chem !------------- * initialize: !------------- copchem=' ' rdtcopchem=0. totr = 0. nr = 0 rk298 = 0. !rk298 = 1.0E-20 no3prod = 0 known_species = 0 nch = 0 tempkc = ' ' coprod = ' ' pp = ' ' pchem = ' ' flag = 0 rxnflg = 1 rate = 0. tarrhc=0. arrhc=0. rarrhc=0. chem=' ' i = lco + INDEX(rdct(lco+1:lcf),' ') chem=rdct(lco+1:i) ! nc is the size of the carbon skeleton ! i.e. number of groups in species nc = 0 DO i = 1, mca IF(group(i)(1:1) .ne. ' ') nc = nc + 1 ENDDO *If rings exist remove ring-join characters from groups IF (nring.gt.0) THEN CALL rjgrm(nring,group,rjg) ENDIF tgroup = group tbond = bond * check that calling "ozone reaction" make sense IF (INDEX(chem,'Cd').EQ.0) THEN mesg = 'the routine was called with no C=C bond' CALL errexit(progname,mesg,chem) ENDIF * ------------------------------------------------------------ * check if the species is known in the rate constant data base * ------------------------------------------------------------ DO i=1,no3dat IF (chem.EQ.o3chem(i)) THEN rk298=o3rate(i,1)*298**o3rate(i,2)*exp(-o3rate(i,3)/298.) * if database rate constant < 1.0E-20, no reaction IF (rk298.LT.1.0E-20) THEN IF (wtflag.EQ.1) WRITE(10,*)'returning from',progname RETURN ENDIF rarrhc(:)=o3rate(i,:) ENDIF ENDDO * ----------------------------------------------------------- * Check file o3_prod.dat for any known products & branching ratios * * WARNING: for those reactions, yields < 5% are not automatically * * eliminated * * ----------------------------------------------------------- IF (rk298.NE.0.) THEN DO i=1,nkwo3 IF (chem.EQ.kwo3rct(i)) THEN known_species = 1 no3prod = nkwo3pd(i) * loop through known products DO j=1,no3prod CALL addreac(nr,progname,rdct(lco+1:lcf),flag) tempkc = kwo3pd(i,j) ratio(j) = kwo3yld(i,j) * check the products IF (index(tempkc,'.').ne.0) THEN CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc) pchem(j,1)=rdckprod(1) coprod(j,1:mca) = rdcktprod(1,1:mca) ELSE pchem(j,1)= tempkc ENDIF CALL stdchm(pchem(j,1)) * add known coproduct(s) if necessary DO k=1,mcp ! loop possible coproducts IF (kwo3copd(i,j,k).NE.' ') THEN DO 455 l=1,mca IF (coprod(j,l).ne.' ') GOTO 455 coprod(j,l) = kwo3copd(i,j,k) WRITE(6,*) 'pp(',l,')=',pp(l) 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 !!!!!!! NITRO ALKENES !!!!!! ! nitro-alkene compounds can be formed during the aromatic oxidation ! processes. We assume that the O3 reactivity is slow compared to the OH ! reactivity, then we ignore ozonolysis of nitro alkenes for now. There ! is still a lack of data concerning nitro alkenes, and this is ! something that may be updated in the future if new informations arise. IF (INDEX(chem,'Cd(NO2)').NE.0) RETURN IF (INDEX(chem,'CdH(NO2)').NE.0) RETURN * ----------------------------------------------------------- * send to the reaction subroutine * ----------------------------------------------------------- * * The parameterization used to assign a rate constant for O3 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 -C=C-C=O type). Subroutine cdcase * returns the "case" to which the species belongs. 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 : >C=C-C=O structure but no C=C-C=C * CASE 3 : -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 : -CO-C=C-C=C< structure (i.e. containing carbonyl * at only one side of the conjugated C=C-C=C * CASE 6 : O=C=C-R reactions not considered, even if R is C=C-R ! Ludo: ! CASE 7 : reactions of ozone with vinyl ether (-O-C=C) and dihydrofurans * find the position of the Cd carbon in the molecule and check if the * molecule is allowed (i.e. no >C=C=C<) 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 * modified with new SAR * ------- parent_group = 0 IF (ncdcase.eq.1) THEN IF (criegee_fg .eq. 2) THEN CALL o3add_c1_cmv(chem,bond,group,ncd,conjug,cdtable,cdsub, & nr,flag,rate,pchem,arrhc, parent_group) ELSE CALL o3add_c1(chem,bond,group,ncd,conjug,cdtable,cdsub, & nr,flag,rate,pchem,arrhc) ENDIF * case 2 : C=C-C=O bonds, only 1 C=C (i.e no conjugated Cds) * ------- ELSE IF (ncdcase.eq.2) THEN CALL o3add_c2(chem,bond,group, & ncd,conjug,cdtable,cdsub,cdcarbo, & nr,flag,rate,pchem, parent_group) * case 3 : -CO-C=C-C=C-C=O- structure : no reaction (much too slow) * ------- ELSE IF (ncdcase.eq.3) THEN IF (wtflag.EQ.1) WRITE(10,*)'returning from',progname RETURN * 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) 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 o3add_c2(chem,bond,group, & ncd,conjug,cdtable,cdsub,cdcarbo, & nr,flag,rate,pchem, parent_group) ELSE IF (cdcarbo(i+1,1).ne.0) THEN CALL o3add_c2(chem,bond,group, & ncd,conjug,cdtable,cdsub,cdcarbo, & nr,flag,rate,pchem, parent_group) ELSE IF (criegee_fg .eq. 2) THEN CALL o3add_c1_cmv(chem,bond,group,ncd,conjug,cdtable,cdsub, & nr,flag,rate,pchem,arrhc, parent_group) ELSE CALL o3add_c1(chem,bond,group,ncd,conjug,cdtable,cdsub, & nr,flag,rate,pchem,arrhc) ENDIF ENDIF ENDDO * case 5 : -CO-C=C-C=C< * ------- ELSE IF (ncdcase.eq.5) THEN CALL o3add_c5(chem,bond,group, & ncd,conjug,cdtable,cdsub,cdcarbo, & nr,flag,rate,pchem) * case 6 : -C=C=O * ------- ! ozone reaction extremely slow: assume no reaction ! case 7 : -C=C-O- ! ---------------- ELSE IF (ncdcase.EQ.7) THEN CALL o3add_c7(chem,bond,group, & ncd,conjug,cdtable,cdsub, & cdeth,nr,flag,rate,pchem) ENDIF * ----------------------------------------------------------- * end reaction * ----------------------------------------------------------- * compute total rate constant totr=0. totr1=0. DO i=1,nr totr = totr + rate(i)*EXP(-arrhc(i,3)/298.) totr1 = totr1 + rate(i) ENDDO * if total rate constant less than 1.0E-20, then no reaction IF (totr.LT.1.0E-20) THEN WRITE(99,'(a)') '--warning--' WRITE(99,'(a)') 'O3 rate constant estimated less than 1E-20' WRITE(99,'(a)') 'reaction with O3 is not considered for:' WRITE(99,'(a)') rdct(lco+1:lcf) ENDIF * flag down for duplicate products IF (nr.gt.1) THEN DO i=1,nr-1 DO j= i+1,nr IF ( ( (pchem(i,1).EQ.pchem(j,1)) .AND. & (pchem(i,2).EQ.pchem(j,2)) ) .OR. & ( (pchem(i,1).EQ.pchem(j,2)) .AND. & (pchem(i,2).EQ.pchem(j,1)) ) ) THEN flag(j) = 0 rate(i) = rate(i) + rate(j) ENDIF ENDDO ENDDO ENDIF * active channel counter DO i=1,nr ! IF (rate(i).LT.cut_off*totr) flag(i) = 0 IF (rate(i)*EXP(-arrhc(i,3)/298.).LT.cut_off*totr) flag(i) = 0 ENDDO * set branching ratio for each channel nch = 0 sum = 0. DO i=1,nr IF (flag(i).EQ.1) THEN nch=nch+1 sum = sum + rate(i) ENDIF ENDDO DO i=1,nr IF (flag(i).EQ.1) THEN ratio(i)=rate(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 else rescale tarrhc(i,1) IF (rk298.NE.0.) THEN IF(rarrhc(1).GT.0.) WRITE(70,*) rdct(1:6),' ',rarrhc(:) 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 (wtflag.NE.0) WRITE(77,*) rk298,sum,' ',chem ELSE IF(totr1.GT.0.) WRITE(70,*) rdct(1:6),' ',totr1,rarrhc(2:3) DO i=1,nr IF (flag(i).EQ.1) THEN ! tarrhc(i,1)=totr*ratio(i) tarrhc(i,1)=totr1*ratio(i) tarrhc(i,2)=0. ! tarrhc(i,3)=0. tarrhc(i,3)=arrhc(i,3) ENDIF ENDDO ENDIF * -------------------- * output * -------------------- ich = 10 nt =0. E = 0. 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,f5.1,1X,f7.0)') & '****INIT GO3 +','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=totr ENDIF fratio=ratio(i) r(1) = rdct(1:lco) r(2) = 'O3 ' * separate hot-criegee from carbonyl products (kk=criegge, ii=carbonyl) IF (INDEX(pchem(i,1),hot_criegee).NE.0) THEN kk = 1 ii = 2 ELSE kk = 2 ii = 1 ENDIF * load non-criegee products: np = 1 s(1) = 1. brtio = brch * ratio(i) * CALL bratio(pchem(i,ii),brtio,p(1),copchem,rdtcopchem, CALL bratio(pchem(i,ii),brtio,p(1), & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) IF(wtopeflag.EQ.1)WRITE(10,'(f5.3,2X,A1,A6)')ratio(i),'G',p(1) IF (known_species.EQ.0) THEN * find products from hot Criegee: ! correct nc if species is CH2CH2 IF(criegee_fg .eq. 2) THEN IF(chem(1:10) .eq. 'CdH2=CdH2 ') nc = 1 CALL xcrieg_cmv(pchem(i,kk), & nc, bond, parent_group(i), & brtio,ss,pp, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ELSE CALL xcrieg(pchem(i,kk), & brtio,ss,pp, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ENDIF * load to output arrays: DO j=1,mnp IF (pp(j)(1:1).NE.' ') THEN CALL addprod(np,progname,chem) s(np) = ss(j) p(np) = pp(j) IF (wtopeflag.EQ.1) & WRITE(10,'(f5.3,2X,A1,A6)') s(np)*ratio(i),'G',p(np) ENDIF ENDDO IF (rdtcopchem.GT.0.) THEN CALL addprod(np,progname,chem) 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 ELSEIF (known_species.EQ.1) THEN 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)' 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 ENDIF * write out: reaction with O3 are thermal reaction (idreac=0 and does * not require labels) CALL rxwrit3(17,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ! DEBUG ! ! CALL rxwrit3(6,a1,a2,a3,a4,r,s,p,ar1,ar2,ar3, ! & f298,fratio,idreac,nlabel,xlabel,folow,fotroe) ! DEBUG ! 300 CONTINUE * end of rozone IF (wtopeflag.EQ.1) WRITE(10,*)'end',progname !IF (wtflag.NE.0) WRITE(*,*) progname,': end : ' ! DEBUG ! ! PRINT*,"rozone2 debug return stop" ! STOP ! DEBUG ! RETURN END