************************************************************************ * MASTER MECHANISM - ROUTINE NAME : no3add_c1 * * * * * * PURPOSE : * * This subroutine computes the reaction rates for NO3 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 reaction products are returned in pchem for the main organic * * product and in coprod for coproducts of the reation. * * * * INPUT: * * - chem : chemical formula * * - group(i) : groups at position (carbon) i * * - bond(i,j) : carbon-carbon bond matrix of chem * * - ncd : number of "Cd" carbons in chem * * - conjug : =1 if conjugated Cd (C=C-C=C), otherwise =0 * * - cdtable(i) : carbon number bearing a "Cd" * * - cdsub(i) : number of -C- substitutents (including -CO-) bonded * * to the Cd corresponding to cdtable(i) * * * * INPUT/OUTPUT * * - nr : number of reaction channels associated with chem * * - flag(i) : flag that active the channel i * * - rno3(i,3) : reaction rate (298 K) for channel i * * - pchem(i) : the main organic product of reaction channel i. * * - coprod(i,j) : coproducts of reaction channel i * * * ************************************************************************ SUBROUTINE no3add_c1(chem,bond,group,ncd,conjug,cdtable,cdsub, & nr,flag,rno3,pchem,coprod,flag_del,pchem_del, & coprod_del,sc_del) 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) REAL rno3(mnr) CHARACTER(LEN=lco) coprod(mnr,mca) * internal INTEGER i,j,k,j1,nc,nb,ic1,ic2,nb1,nb2 INTEGER tbond(mca,mca),nring CHARACTER(LEN=lgr) tgroup(mca), pold, pnew CHARACTER(LEN=lfo) tempgr CHARACTER(LEN=lco) tprod(mca) REAL fract(4) REAL rtot REAL beta,delta INTEGER ialpha, ibeta, igamma, idelta 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) REAL arrhc(3) CHARACTER(LEN=lsb) :: progname= '*no3add_c1*' CHARACTER(LEN=ler) :: mesg * write info for finding bugs IF (wtflag.NE.0) WRITE(*,*) progname * ----------- * initialize * ----------- tgroup = group tbond = bond * 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 as follows : * CH2=CH2 : 2.05E-16 * C-CH=CH2 : rate const. for 1 butene (1.35E-14) * C-C(C)=CH2 : rate const. for 2 methyl propene (3.32E-13) * C-CH=CH-C : rate const. for 2 butene (3.7E-13) (mean value cis-trans) * C-C(C)=CH-C : rate const. for 2 methyl 2 butene (9.37E-12) * c-C(C)=C(C)-C : rate const. for 2,3 dimethyl 2 butene (5.72E-11) * Branching ratio is set assuming 100 % yield of the more stable radical * produced * check that there is no problem with the Cd case IF (conjug.ne.0) THEN WRITE(6,*) '-error in no3add_c1, conjug eq 0 expected ' WRITE(99,*) 'no3add_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 CALL raddno3(tbond,tgroup,cdtable(ic1),cdtable(ic2),arrhc) rtot = arrhc(1) * -------------------------------------- * set rate constant and branching ratio * -------------------------------------- * assign the rate constant c nb=cdsub(ic1)+cdsub(ic2) c IF (nb.eq.0) rtot=2.05E-16 c IF (nb.eq.1) rtot=1.35E-14 c IF (nb.eq.2) THEN c IF (cdsub(ic1).eq.1) THEN c rtot=3.70E-13 c ELSE c rtot=3.32E-13 c ENDIF c ENDIF c IF (nb.eq.3) rtot=9.37E-12 c IF (nb.eq.4) rtot=5.72E-11 * assign the branching ratio IF (cdsub(ic1).eq.cdsub(ic2)) THEN fract(ic1)=0.5 fract(ic2)=0.5 ELSE IF (cdsub(ic1).gt.cdsub(ic2)) THEN fract(ic2)=1. fract(ic1)=0. ELSE fract(ic2)=0. fract(ic1)=1. ENDIF * ---------------- * do the reaction * ---------------- DO i=ic1,ic2 IF (fract(i).gt.0.) THEN CALL addreac(nr,progname,chem,flag) rno3(nr) = rtot*fract(i) c totr = totr + rno3(nr) * convert to single carbon bonds: tbond(cdtable(ic1),cdtable(ic2)) = 1 tbond(cdtable(ic2),cdtable(ic1)) = 1 pold = 'Cd' pnew = 'C' tempgr=group(cdtable(ic1)) CALL swap(tempgr,pold,tgroup(cdtable(ic1)),pnew) tempgr=group(cdtable(ic2)) CALL swap(tempgr,pold,tgroup(cdtable(ic2)),pnew) * add (ONO2) to cdtable(I) carbon, nc = INDEX(tgroup(cdtable(i)),' ') tgroup(cdtable(i))(nc:nc+5) = '(ONO2)' * add radical dot to the other carbon: IF (i.eq.ic1) THEN j1=ic2 ELSE j1=ic1 ENDIF nc = INDEX(tgroup(cdtable(j1)),' ') tgroup(cdtable(j1))(nc:nc) = '.' * rebuild, check, and rename: CALL rebond(tbond,tgroup,tempgr,nring) c CALL radchk(tempgr,pchem(nr),tprod) CALL radchk(tempgr,rdckprod,rdcktprod,nip,sc) pchem(nr) = rdckprod(1) sc_del(nr,1) = sc(1) IF (nip.EQ.2) STOP CALL stdchm(pchem(nr)) coprod(nr,:) = rdcktprod(1,:) * reset groups and bonds: tgroup = group tbond = bond ENDIF ENDDO 40 CONTINUE 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 asigned * by considering first the alkyl group at terminal position. A correction * is then made to account for the contribution of the internal alkyl group. * Rate constants are assigned as follows : * CH2=CH-CH=CH2 : 1.00E-13 * C-CH=CH-CH=CH2 : 1.50E-12 (from CH3-CH=CH-CH=CH2) * C-CH=CH-CH=CH-C : 1.60E-11 (from CH3-CH=CH-CH=CH-CH3) * and then multiply by 6.8 if there is 1 internal alkyl group (ratio * isoprene / 1,3 butadiene) and 21 if there are 2 internal alkyl groups * (ratio 2,3 dimethyl / 1,3 butadiene and butadiene) * 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 : NO3 is assumed to add at the less substitued terminal * position. If the 2 terminal positions are equivalent, then the dominant * position is assumed to be the one yielding the highest substitued radical * If label 102 is reached, then the molecule has a C=C-C=C structure 102 CONTINUE * -------------------------------------- * set rate constant and branching ratio * -------------------------------------- * set the rate constant - contribution of external alkyl group c nb=cdsub(1)+cdsub(4) c IF (nb.eq.0) rtot=1.00E-13 c IF (nb.eq.1) rtot=1.50E-12 c IF (nb.gt.1) rtot=1.60E-11 * contribution of internal alkyl group (nb=3 means 1 internal alkyl group) c nb=cdsub(2)+cdsub(3) c IF (nb.eq.3) rtot=rtot*6.8 c IF (nb.gt.3) rtot=rtot*21. * assign the branching ratio for NO3 addition c IF (cdsub(1).eq.cdsub(4)) THEN c nb1=cdsub(2)+cdsub(4) c nb2=cdsub(1)+cdsub(3) c IF (nb1.eq.nb2) THEN c fract(1)=0.5 c fract(4)=0.5 c ELSE IF (nb1.gt.nb2) THEN c fract(1)=1. c fract(4)=0. c ELSE c fract(1)=0. c fract(4)=1. c ENDIF c ELSE IF (cdsub(1).gt.cdsub(4)) THEN c fract(4)=1. c fract(1)=0. c ELSE c fract(4)=0. c fract(1)=1. c ENDIF ! ric : the new SAR estimates rate constant for each double bond fract(1)=1. fract(4)=1. * ---------------- * do the reaction * ---------------- * NO3 adds at terminal position (ialpha position), radical dot is at * beta position (ibeta) or delta position (idelta). DO i=1,4,3 ialpha=i IF (ialpha.eq.1) THEN ibeta=2 igamma=3 idelta=4 ELSE ibeta=3 igamma=2 idelta=1 ENDIF ! ric : set the rate constant - kerdouci et al., 2010 CALL raddno3(tbond,tgroup,cdtable(ialpha),cdtable(ibeta), & arrhc) rtot = arrhc(1) IF (fract(ialpha).ne.0.) THEN IF (cdsub(ibeta).eq.cdsub(idelta)) THEN beta=0.5 delta=0.5 ELSE IF (cdsub(ibeta).gt.cdsub(idelta)) THEN beta=1. delta=0. ELSE beta=0. delta=1. ENDIF * ric : do the beta addition, delocalisation is now done in radchk beta=1. delta=0. * BETA ADDITION * ------------- IF (beta.ne.0.) THEN CALL addreac(nr,progname,chem,flag) rno3(nr) = rtot*fract(ialpha)*beta * convert to single carbon bonds: tbond(cdtable(ialpha),cdtable(ibeta)) = 1 tbond(cdtable(ibeta),cdtable(ialpha)) = 1 pold = 'Cd' pnew = 'C' tempgr=group(cdtable(ialpha)) CALL swap(tempgr,pold,tgroup(cdtable(ialpha)),pnew) tempgr=group(cdtable(ibeta)) CALL swap(tempgr,pold,tgroup(cdtable(ibeta)),pnew) * add (ONO2) to cdtable(i) carbon, nc = INDEX(tgroup(cdtable(ialpha)),' ') tgroup(cdtable(ialpha))(nc:nc+5) = '(ONO2)' * add radical dot to the other carbon: nc = INDEX(tgroup(cdtable(ibeta)),' ') tgroup(cdtable(ibeta))(nc:nc) = '.' * rebuild, check, and rename: CALL rebond(tbond,tgroup,tempgr,nring) c CALL radchk(tempgr,pchem(nr),tprod) CALL radchk(tempgr,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) CALL stdchm(pchem_del(nr)) coprod_del(nr,:) = rdcktprod(2,:) ENDIF CALL stdchm(pchem(nr)) coprod(nr,:) = rdcktprod(1,:) * reset groups and bonds: tgroup = group tbond = bond ENDIF * DELTA ADDITION * --------------- c IF (delta.ne.0.) THEN IF ((delta.ne.0.).AND.(delta.NE.0.5)) THEN !print*,'delta' !print*,i,':',ialpha,ibeta,igamma,idelta CALL addreac(nr,progname,chem,flag) rno3(nr) = rtot*fract(ialpha)*delta * convert to single carbon bonds: tbond(cdtable(ialpha),cdtable(ibeta)) = 1 tbond(cdtable(ibeta),cdtable(ialpha)) = 1 tbond(cdtable(igamma),cdtable(idelta)) = 1 tbond(cdtable(idelta),cdtable(igamma)) = 1 pold = 'Cd' pnew = 'C' tempgr=group(cdtable(ialpha)) CALL swap(tempgr,pold,tgroup(cdtable(ialpha)),pnew) tempgr=group(cdtable(idelta)) CALL swap(tempgr,pold,tgroup(cdtable(idelta)),pnew) * add (ONO2) to cdtable(i) carbon, nc = INDEX(tgroup(cdtable(ialpha)),' ') tgroup(cdtable(ialpha))(nc:nc+5) = '(ONO2)' * add radical dot to the other carbon: nc = INDEX(tgroup(cdtable(idelta)),' ') tgroup(cdtable(idelta))(nc:nc) = '.' * add the double bond between beta and gamma position tbond(cdtable(ibeta),cdtable(igamma)) = 2 tbond(cdtable(igamma),cdtable(ibeta)) = 2 * rebuild, check, and rename: CALL rebond(tbond,tgroup,tempgr,nring) c CALL radchk(tempgr,pchem(nr),tprod) CALL radchk(tempgr,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) CALL stdchm(pchem_del(nr)) coprod_del(nr,:) = rdcktprod(2,:) ENDIF CALL stdchm(pchem(nr)) coprod(nr,:) = rdcktprod(1,:) * reset groups and bonds: tgroup = group tbond = bond ENDIF ENDIF ENDDO !IF (wtflag.NE.0) WRITE(*,*) progname,': end :' RETURN END