************************************************************************ * MASTER MECHANISM - ROUTINE NAME : o3add_c1 * * * * * * PURPOSE : * * This subroutine computes the reaction rate for O3 addition to * * >C=C< bond (case 1). Species having a C=O group conjugated with * * the C=C are not taken into account here * * * * Part I performs the reaction for simple >C=C< bond ; Part II for * * the conjugated >C=C-C=C< bonds. The method used for each case is * * described in the body of the program. * * * * The 2 reaction products (carbonyl and criegge) are returned as * * pchem(i,1) and pchem(i,2). * * * * INPUT: * * - chem : chemical formula * * - group(i) : groups at position (carbon) i * * - bond(i,j) : carbon-carbon bond matrix of chem * * - ncd : number of "Cd" carbon in chem * * * * INPUT/OUTPUT * * - nr : number of reaction channels associated with chem * * - flag(i) : flag for active channel i * * - tarrhc(i,3) : arrhenius coefficient for channel i * * - pchem(i,2) : the 2 main organic products of reaction channel i * * - rate(i) : reaction rate constant * * * ************************************************************************* SUBROUTINE o3add_c1(chem,bond,group,ncd,conjug,cdtable,cdsub, & nr,flag,rate,pchem,tarrhc) IMPLICIT NONE INCLUDE 'general.h' INCLUDE 'common.h' * input: CHARACTER(LEN=lfo) chem INTEGER bond(mca,mca) CHARACTER(LEN=lgr) group(mca) INTEGER ncd, conjug INTEGER cdtable(4),cdsub(4) * input/output INTEGER nr, flag(mnr) CHARACTER(LEN=lfo) pchem(mnr,2) REAL rate(mnr) ! output REAL,INTENT(out) :: tarrhc(mnr,3) * internal INTEGER i,j,k,l,j1,nc,nb,ic1,ic2 INTEGER tbond(mca,mca) CHARACTER(LEN=lgr) tgroup(mca), pold, pnew CHARACTER(LEN=lfo) tempkg CHARACTER(LEN=lco) tprod(mca) REAL fract(4) REAL rtot INTEGER nring, ndis INTEGER nca,m,n, pos1, pos2 INTEGER rngflg INTEGER ring(mca) INTEGER last,ncx REAL :: arrhc(3) !INTEGER,PARAMETER :: magnify_fg = 1 CHARACTER(LEN=lsb) :: progname='*o3add_c1* ' CHARACTER(LEN=ler) :: mesg * write info for finding bugs IF (wtflag.NE.0) WRITE(*,*) progname * Initialize * ----------- tgroup(:)=group(:) tbond(:,:)=bond(:,:) tarrhc(:,:) = 0 nca=0 last=0 DO i=1,mca IF (group(i)(1:1).NE.' ') THEN nca=nca+1 last=i ENDIF ENDDO ncx = 0 DO i=1,last DO j=i+1,last IF(bond(i,j).GT.0) ncx = ncx + 1 ENDDO ENDDO nring=ncx-nca+1 * Send to Part II if the molecule has conjugated Cds IF (conjug.eq.1) GOTO 102 * ================================================================== * PART I : non conjugated Cds (simple >C=C<) structures * ================================================================== * * Rate constants are assigned for the C=C bond, not for each Cd. * Evaluation is made using group rate constants provided in the SAPRC99 * mechanism. Values are mean values based on the rate constants provided * in Atkinson, 1997, J. Phys. Chem. Ref. Data. * CH2=CH2 : 1.68E-18 * C-CH=CH2 : 1.01E-17 * C-C(C)=CH2 : 1.18E-17 * C-CH=CH-C : 1.15E-16 * C-C(C)=CH-C : 3.48E-16 * C-C(C)=C(C)-C : 1.13E-15 * Branching ratio is set assuming 50 % for both channels for a symmetrical * C=C bond (in the sense of the group defined above). For a non-symmetrical * C=C bond, the following is used (which is almost identical to yields * provided in the SAPRC99 chemical scheme): * C-CH=CH2 : 50 % for each side (see data in Atkinson, 1997) * C-C(C)=CH2 : 65 % (HCHO + >C.(OO.)) , 35 % (>C=O + CH2.(OO.)) * C-C(C)=CH-C : 30 % (>C=O + R-CH.(OO.)) , 70 % (R-CHO + >C.(OO.)) * check that there is no problem with the Cd case IF (conjug.ne.0) THEN WRITE(6,*) '-error in O3add_c1, conjug eq 0 expected ' WRITE(99,*) 'o3add_c1',chem !STOP ENDIF * Loop over the 2 possible double bonds. In cdtable and other related * tables, table elements 1-2 and 3-4 are double bonded. If only 1 C=C * bond, then exit after treatment. DO 40 ic1=1,3,2 ic2=ic1+1 IF (cdtable(ic1).eq.0) GOTO 40 * assign the rate constant and branching ratio * -------------------------------------------- nb=cdsub(ic1)+cdsub(ic2) * ethene IF (nb.eq.0) THEN rtot=1.68E-18 fract(ic1)=0.5 fract(ic2)=0.5 * -CH=CH2 ELSE IF (nb.eq.1) THEN rtot=1.01E-17 fract(ic1)=0.5 fract(ic2)=0.5 * -CH=CH- or >C=CH2 ELSE IF (nb.eq.2) THEN IF (cdsub(ic1).eq.1) THEN rtot=1.15E-16 fract(ic1)=0.5 fract(ic2)=0.5 ELSE rtot=1.18E-17 IF ( (cdsub(ic1).eq.0).AND.(cdsub(ic2).eq.2) ) THEN fract(ic1)=0.65 fract(ic2)=0.35 ELSE IF ( (cdsub(ic1).eq.2).AND.(cdsub(ic2).eq.0) ) THEN fract(ic1)=0.35 fract(ic2)=0.65 ENDIF ENDIF * >C=CH- ELSE IF (nb.eq.3) THEN rtot=3.48E-16 IF ( (cdsub(ic1).eq.1).AND.(cdsub(ic2).eq.2) ) THEN fract(ic1)=0.70 fract(ic2)=0.30 ELSE IF ( (cdsub(ic1).eq.2).AND.(cdsub(ic2).eq.1) ) THEN fract(ic1)=0.30 fract(ic2)=0.70 ELSE WRITE(6,*) '-error in O3add_c1, no ratio found ' WRITE(99,*) 'o3add_c1',chem !STOP ENDIF * >C=C< ELSE IF (nb .eq.4) THEN rtot=1.13E-15 fract(ic1)=0.5 fract(ic2)=0.5 ELSE WRITE(6,*) '-error1-- in O3add_c1, no rate constant found ' WRITE(99,*) 'o3add_c1',chem STOP ENDIF IF (hoadd_c1_fg.EQ.3) THEN CALL o3rate(bond,group,nca,conjug,cdtable,cdsub, & ic1,ic2,arrhc) rtot=arrhc(1) ENDIF * perform the reaction * -------------------- DO i=ic1,ic2 IF (fract(i).gt.0.) THEN nr = nr + 1 IF (nr.GT.mnr) THEN WRITE(6,'(a)') '--error--' WRITE(6,'(a)') 'from MASTER MECHANISM ROUTINE:rozone' WRITE(6,'(a)') 'too many reactions created for:' WRITE(6,'(a)') chem WRITE(99,*) 'o3add_c1',chem !STOP ENDIF flag(nr) = 1 rate(nr) = rtot*fract(i) tarrhc(nr,1) = arrhc(1) *fract(i) tarrhc(nr,3) = arrhc(3) * check if bond that will be brocken belong to a cycle. pos1=cdtable(ic1) pos2=cdtable(ic2) CALL findring(pos1,pos2,nca,tbond,rngflg,ring) * break the C=C bond and change "Cd" to "C" tbond(cdtable(ic1),cdtable(ic2)) = 0 tbond(cdtable(ic2),cdtable(ic1)) = 0 pold = 'Cd' pnew = 'C' tempkg=group(cdtable(ic1)) CALL swap(tempkg,pold,tgroup(cdtable(ic1)),pnew) tempkg=group(cdtable(ic2)) CALL swap(tempkg,pold,tgroup(cdtable(ic2)),pnew) * change cdtable(i) carbon to carbonyl IF (tgroup(cdtable(i))(1:3).EQ.'CH2') THEN pold = 'CH2' pnew = 'CH2O' ELSE IF (tgroup(cdtable(i))(1:2).EQ.'CH') THEN pold = 'CH' pnew = 'CHO' ELSE IF (tgroup(cdtable(i))(1:1).EQ.'C') THEN pold = 'C' pnew = 'CO' ENDIF CALL swap(tgroup(cdtable(i)),pold,tempkg,pnew) tgroup(cdtable(i)) = tempkg * change the other carbon to hot criegee: IF (i.eq.ic1) THEN j1=ic2 ELSE j1=ic1 ENDIF nc = INDEX(tgroup(cdtable(j1)),' ') tgroup(cdtable(j1))(nc:nc+6) = '.(OO.)*' * Check the number of distinct molecule (could be 1 because of * ring structure, eg terpenes). If 2 then fragments the species into 2 parts. c CALL chk_nmol(tbond,tgroup,ndis) IF (rngflg.eq.1) THEN CALL rejoin(nca,pos1,pos2,m,n,tbond,tgroup) CALL rebond(tbond,tgroup,tempkg,nring) pchem(nr,1)=tempkg pchem(nr,2)=' ' CALL stdchm(pchem(nr,1)) ELSE IF (rngflg.eq.0) THEN CALL fragm(tbond,tgroup,pchem(nr,1),pchem(nr,2)) CALL stdchm(pchem(nr,1)) CALL stdchm(pchem(nr,2)) ELSE WRITE(6,*) '--error--, number of distinct molecule is' WRITE(6,*) 'not equal to 1 or 2 in o3add_c1 (see rngflg)' STOP ENDIF * reset groups and bonds: tgroup(:) = group(:) tbond(:,:) = bond(:,:) ENDIF ENDDO 40 CONTINUE ! DEBUG ! ! print*,"mid o3add_c1 debug stop" ! STOP ! DEBUG ! RETURN * ================================================================= * PART II : conjugated Cds (>C=C-C=C<) structures * ================================================================= * Alkyl group at terminal position (i.e. alpha and delta) seems to have * a greater impact on rate constant than alkyl group at internal * position (i.e beta and gamma). Therefore, rate constant is assigned * considering first the alkyl group at terminal position. A correction * is then made to take into account the contribution of the internal * alkyl group. Rate constants are assigned as follows : * CH2=CH-CH=CH2 : 6.3E-18 * C-CH=CH-CH=CH2 : 3.5E-17 (from CH3-CH=CH-CH=CH2 mean value E/Z) * C-C(C)=CH-CH=CH2: 3.9E-17 (10 % increase from above, equal to the * difference between propene and isobutene) * C-CH=CH-CH=CH-C : 3.4E-16 (from CH3-CH=CH-CH=CH-CH3) * >CH=CH-CH=CH-C : 3.4E-16 (from CH3-CH=CH-CH=CH-CH3) * >CH=CH-CH=CH< : 3.1E-15 (from CH3-C(CH3)=CH-CH=C(CH3)CH3 Lewin 01) * Then multiply by 2.03 for each internal alkyl group (i.e. the ratio of * isoprene/butadiene). Note: agrees well with (2,3-dimethyl-1,3-butadiene * versus isoprene ratio). ALL THIS IS OF COURSE A VERY CRUDE ESTIMATE. * (more work should be done in the future to better estimate these rate * constants). * For branching ratio : from branching ratio for simple alkene (and data * coming from isoprene study), it can be expected TO VERY FIRST * APPROXIMATION that all possible reactions have a similar yield. Here, * we assume that each reaction occurs with 25 % yield. ALL THIS IS OF * COURSE A VERY CRUDE ESTIMATE. (more work should be done in the future * to better estimate these yields). * If label 102 is reached, then the molecule has a C=C-C=C structure 102 CONTINUE * set the rate constant - contribution of external alkyl groups nb=cdsub(1)+cdsub(4) IF (nb.eq.0) THEN rtot=6.3E-18 ELSE IF (nb.eq.1) THEN rtot=3.5E-17 ELSE IF (nb.eq.2) THEN IF (cdsub(1).eq.1) THEN rtot=3.4E-16 ELSE rtot=3.9E-17 ENDIF ELSE IF (nb.eq.3) THEN rtot=3.4E-16 ELSE IF (nb.eq.4) THEN rtot=3.1E-15 ELSE WRITE(6,*) '-error2-- in O3add_c1, rate constant not found ' WRITE(99,*) 'o3add_c1',chem !STOP ENDIF * internal alkyl group contribution (nb=3 means 1 internal alkyl group) nb=cdsub(2)+cdsub(3) IF (nb.eq.3) rtot=rtot*2.03 IF (nb.gt.3) rtot=rtot*2.03*2.03 * set the branching ratio - fract(i) is the branching ratio making the * carbonyl at position i (and the criegge at the other side of the C=C) fract(1:4)=0.25 IF (hoadd_c1_fg.EQ.3) THEN CALL o3rate(bond,group,nca,conjug,cdtable,cdsub, & ic1,ic2,arrhc) rtot=arrhc(1) ENDIF * perform the reactions DO i=1,4 nr = nr + 1 IF (nr.GT.mnr) THEN WRITE(6,'(a)') '--error--' WRITE(6,'(a)') 'from MASTER MECHANISM ROUTINE:3add_c1' WRITE(6,'(a)') 'too many reactions created for:' WRITE(6,'(a)') chem WRITE(99,*) 'o3add_c1',chem !STOP ENDIF flag(nr) = 1 rate(nr) = rtot*fract(i) tarrhc(nr,1)=arrhc(1)*fract(i) tarrhc(nr,3)=arrhc(3) * find the other carbon bonded to the cdtable(i) carbon IF (i.eq.1) j1=2 IF (i.eq.2) j1=1 IF (i.eq.3) j1=4 IF (i.eq.4) j1=3 * check if bond that will be brocken belong to a cycle. pos1=cdtable(i) pos2=cdtable(j1) CALL findring(pos1,pos2,nca,tbond,rngflg,ring) * break the C=C bond and change "Cd" to "C" tbond(cdtable(i),cdtable(j1)) = 0 tbond(cdtable(j1),cdtable(i)) = 0 pold = 'Cd' pnew = 'C' tempkg=group(cdtable(i)) CALL swap(tempkg,pold,tgroup(cdtable(i)),pnew) tempkg=group(cdtable(j1)) CALL swap(tempkg,pold,tgroup(cdtable(j1)),pnew) * change cdtable(i) carbon to carbonyl IF (tgroup(cdtable(i))(1:3).EQ.'CH2') THEN pold = 'CH2' pnew = 'CH2O' ELSE IF (tgroup(cdtable(i))(1:2).EQ.'CH') THEN pold = 'CH' pnew = 'CHO' ELSE IF (tgroup(cdtable(i))(1:1).EQ.'C') THEN pold = 'C' pnew = 'CO' ENDIF CALL swap(tgroup(cdtable(i)),pold,tempkg,pnew) tgroup(cdtable(i)) = tempkg * change the other carbon to hot criegee: nc = INDEX(tgroup(cdtable(j1)),' ') tgroup(cdtable(j1))(nc:nc+6) = '.(OO.)*' * Check the number of distinct molecule (could be 1 because of * ring structure, eg terpenes). If 2 then fragments the species into 2 parts. IF (rngflg.eq.1) THEN CALL rejoin(nca,pos1,pos2,m,n,tbond,tgroup) CALL rebond(tbond,tgroup,tempkg,nring) pchem(nr,1)=tempkg pchem(nr,2)=' ' CALL stdchm(pchem(nr,1)) ELSE IF (rngflg.eq.0) THEN CALL fragm(tbond,tgroup,pchem(nr,1),pchem(nr,2)) CALL stdchm(pchem(nr,1)) CALL stdchm(pchem(nr,2)) ELSE WRITE(6,*) '--error--, number of distinct molecule is' WRITE(6,*) 'not equal to 1 or 2 in o3add_c1 (see rngflg)' STOP ENDIF * reset groups and bonds: DO j=1,mca tgroup(j) = group(j) DO k=1,mca tbond(j,k) = bond(j,k) ENDDO ENDDO ENDDO * exit the conjugated case ! DEBUG ! ! print*,"end o3add_c1 debug stop" ! STOP ! DEBUG ! RETURN END