************************************************************************ * MASTER MECHANISM - ROUTINE NAME : cdcase * * * * * * PURPOSE : 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 returns the "case" to which * * the species belongs. Six 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 : is for structures containing the >C=C-C=O structure * * but no C=C-C=C (or C=C=O) * * CASE 3 : is for the -CO-C=C-C=C-C=O structure only (i.e. containing * * carbonyl at both sides of the conjugated C=C-C=C) * * CASE 4 : Two double bonds non conjugated (i.e. not C=C-C=C) but * * containing at least one C=C-C=O * * CASE 5 : is for the -CO-C=C-C=C< structure (i.e. containing * * carbonyl at only one side of the conjugated C=C-C=C) * * CASE 6 : -C=C=O e.g. ketene (only case for this group) * * CASE 7 : -C=C-O- => vinyl ether chemistry * * * * Nov 20, 2015, J.Lee-Taylor, CIRES: added check for ring structures. * * * * INPUT: * * - chem : chemical formula * * - group(i) : groups at position (carbon) i * * - bond(i,j) : carbon-carbon bond matrix of chem * * - rxnflg : = 0 if calling subroutine requries identification * * only, = 1 if chemistry also required * * - nring : number of distinct rings in the molecule * * * * OUTPUT: * * - cdcase : case to which chem belongs * * - 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", ordered by bonds * * - tcdcase(i) : carbon number at terminal position in C=C-C=C * * - cdsub(i) : number of -C- substitutents (including -CO-) bonded * * to the Cd corresponding to cdtable(i) * * - cdcarbo(i,j) : number of -CO- substitutents bonded to the Cd * * corresponding to cdtable(i) * * * * EDITED * * Mar 6, 2012, J.Lee-Taylor, NCAR: carbons are now listed in cdtable * * in order of their connection (path). Allows for Cd on side branches. * ************************************************************************ SUBROUTINE cdcase(chem,bond,group,rxnflg,nring, & ncd,ncdcase,conjug, & cdtable,tcdtable,cdsub,cdcarbo,cdeth) IMPLICIT NONE INCLUDE 'general.h' INCLUDE 'common.h' INCLUDE 'organic.h' * input: CHARACTER(LEN=lfo),INTENT(in) :: chem INTEGER,INTENT(in) :: bond(mca,mca) CHARACTER(LEN=lgr),INTENT(in) :: group(mca) INTEGER,INTENT(in) :: rxnflg, nring * output : INTEGER,INTENT(out) :: ncd INTEGER,INTENT(out) :: ncdcase, cdtable(4),tcdtable(4), cdsub(4) INTEGER,INTENT(out) :: cdcarbo(4,2) INTEGER,INTENT(out) :: cdeth(4,2) * internal INTEGER :: conjug, l, i , j, k, nb,o INTEGER :: clngth(mca,mca) INTEGER :: path(mca,mca,mca), pathtmp(mca) INTEGER :: top,sec INTEGER :: cnum,onum,nc,nca INTEGER :: dbflg INTEGER :: len INTEGER :: rjg(mri,2) ! ring-join group pairs CHARACTER(LEN=lco) :: prod CHARACTER(LEN=lsb) :: progname='*cdcase* ' CHARACTER(LEN=ler) :: mesg * write info for finding bugs IF (wtflag.NE.0) WRITE(*,*) progname 12 CONTINUE * ----------- * initialize * ----------- ncdcase=0 ncd=0 l=0 cdtable=0 cdsub=0 tcdtable=0 cdcarbo=0 cdeth=0 dbflg = 1 clngth = 0 len = 0 pathtmp = 0 * get the number of >C< and >c< (nca) in the molecule nc = INDEX(chem,' ') - 1 nca = cnum(chem,nc)+onum(chem,nc) * IF RINGS EXIST remove ring-join characters from groups IF (nring.gt.0) THEN CALL rjgrm(nring,group,rjg) ENDIF * count number of Cd in the molecule and put Cd number in cdtable DO i=1,mca IF (INDEX(group(i),'Cd').NE.0) THEN ncd=ncd+1 cdtable(ncd)=i ENDIF ENDDO ! troubleshooting ! print*,'---cdtable-initial-----------' ! print*, (cdtable(j),j=1,ncd) ! print*,'---cdtable bond matrix-------------' ! DO i=1,ncd ! print*,cdtable(i),(bond(cdtable(i),cdtable(j)),j=1,ncd) ! ENDDO ! print*,'----------------' ! use lntreeCd to find longest path starting at a Cd and ! containing all the double bonds DO i=1,ncd top = cdtable(i) DO j=1,ncd IF(bond(top,j).EQ.2) THEN sec = j CALL lntreeCd(bond,dbflg,top,sec,nca,clngth,path) IF(clngth(top,sec).GT.len)THEN len = clngth(top,sec) pathtmp(:) = path(top,sec,:) ! print*,'top, sec, clngth ',top,sec,clngth(top,sec) ! print*,'pathtmp', pathtmp(1:len) ENDIF ENDIF ENDDO ENDDO ! replace cdtable with appropriate carbons from 'pathtmp' (reorder ! cdtable). i=1 DO j=1,len IF (INDEX(group(pathtmp(j)),'Cd').NE.0) THEN cdtable(i) = pathtmp(j) i=i+1 ENDIF ENDDO ! troubleshooting ! print*,'---cdtable-revised-----------' ! print*, (cdtable(j),j=1,ncd) ! print*,'---cdtable bond matrix-------------' ! DO i=1,ncd ! print*,cdtable(i),(bond(cdtable(i),cdtable(j)),j=1,ncd) ! ENDDO ! print*,'----------------' * check that the molecule is allowed. * ----------------------------------- IF (ncd.eq.0) RETURN * check for >C=C=C< structure IF (ncd.eq.3) THEN mesg = '>C=C=C< structure not allowed' CALL errexit(progname,mesg,chem) ENDIF * more than 2 double bonds are not treated here IF (ncd.gt.4) THEN WRITE(6,'(a4,i1)') 'ncd=',ncd mesg = 'number of double bonded carbons > 4' CALL errexit(progname,mesg,chem) ENDIF * check for >C=CR-OH (enol: should have already tautomerised to ketone * in subroutine alkcheck) and >C=CR-ONO2 (not available yet) IF(rxnflg.NE.0)THEN DO i=1,4 IF(cdtable(i).GT.0)THEN IF (INDEX(group(cdtable(i)),'(ONO2)').NE.0)THEN mesg = '>C=CR-ONO2 structure not allowed' CALL errexit(progname,mesg,chem) ENDIF IF (INDEX(group(cdtable(i)),'(OH)').NE.0)THEN mesg = '>C=CR-OH structure not allowed: pass to alkcheck' ! CALL errexit(progname,mesg,chem) CALL alkcheck(chem,prod) IF (INDEX(chem,'=').EQ.0) RETURN nc = INDEX(chem,' ') - 1 CALL grbond(chem,nc,group,bond,dbflg,nring) GOTO 12 ENDIF ENDIF ENDDO ENDIF * check if the molecule is a >C=C-C=C< structure. At the end, tcdtable * should only contain the "terminal Cd's". This is necessary for * further checks conjug=0 IF (ncd.eq.4) THEN tcdtable(1:4)=cdtable(1:4) DO i=1,3 DO j=i+1,4 IF (bond(cdtable(i),cdtable(j)).eq.1) THEN conjug=1 tcdtable(i)=0 tcdtable(j)=0 GOTO 457 ENDIF ENDDO ENDDO ENDIF 457 CONTINUE ! troubleshooting ! print*,'---tcdtable------------' ! print*,tcdtable * count the number of carbons or -O- bonded to each Cd, except the carbon * involved in the C=C bond (for which bond(i,j)=2) DO i = 1, 4 IF (cdtable(i).ne.0) THEN nb=0 DO j=1,mca IF ((bond(cdtable(i),j).eq.1).OR. & (bond(cdtable(i),j).eq.3)) nb=nb+1 ENDDO cdsub(i)=nb * JMLT: if C=C=O present, case 6 IF (INDEX(group(cdtable(i)),'CdO').NE.0) THEN ncdcase=6 GOTO 700 ENDIF ! Ludo: check the number of ether functions on C=C o=0 DO j=1,mca IF (bond(cdtable(i),j).EQ.3) THEN IF (o.LT.1) THEN o=o+1 cdeth(i,o)=j ELSE mesg = 'more than one ether function on a Cd' CALL errexit(progname,mesg,chem) ENDIF ENDIF ENDDO ENDIF ENDDO ! Ludo: if C=C-O- structure, case 7 DO i=1,3,2 IF (((cdeth(i,1)).NE.0).AND.((cdeth(i+1,1)).NE.0)) THEN mesg = 'ether function found both sides of a Cd=Cd' CALL errexit(progname,mesg,chem) ELSE IF ((((cdeth(i,1)).NE.0).AND.((cdeth(i+1,1)).EQ.0)).OR. & (((cdeth(i,1)).EQ.0).AND.((cdeth(i+1,1)).NE.0))) THEN ncdcase=7 WRITE(32,*) 'case 7',chem(1:100) GOTO 700 ENDIF ENDDO * count the number of carbonyls conjugated to each Cd DO i=1,4 l=0 IF (cdtable(i).ne.0) THEN DO j=1,mca IF (bond(cdtable(i),j).eq.1) THEN ! IF (group(j).eq.'CHO') THEN IF (group(j).eq.'CHO' .OR. & group(j)(1:2).eq.'CO') THEN l=l+1 cdcarbo(i,l)=j ENDIF ENDIF ENDDO ENDIF IF (l.gt.2) THEN mesg = 'more than 2 carbonyls bonded to a Cd' CALL errexit(progname,mesg,chem) ENDIF ENDDO * check if the molecule has a -CO-C=C-C=C-CO- structure (case 3) IF (conjug.eq.1) THEN l=0 DO 345 i=1,4 IF (tcdtable(i).ne.0) THEN DO j=1,mca IF (bond(tcdtable(i),j).eq.1) THEN IF (group(j)(1:3).EQ.'CHO') THEN l=l+1 c Ci=tcdtable(i) GOTO 345 ENDIF IF (group(j)(1:2).EQ.'CO') THEN l=l+1 c Ci=tcdtable(i) GOTO 345 ENDIF ENDIF ENDDO ENDIF 345 CONTINUE IF (l.ge.2) THEN ncdcase=3 GOTO 700 ENDIF ENDIF * -C=C-C=C-CO- is substituted with at least one carbonyl (case 5) IF (conjug.eq.1) THEN DO i=1,4 IF (cdtable(i).ne.0) THEN DO j=1,mca IF (bond(cdtable(i),j).eq.1) THEN IF (group(j)(1:3).EQ.'CHO') THEN IF ((cdcarbo(2,1).NE.0).AND. & (cdcarbo(3,1).NE.0)) THEN ncdcase=3 GOTO 700 ELSE ncdcase=5 GOTO 700 ENDIF ENDIF IF (group(j)(1:2).EQ.'CO') THEN IF ((cdcarbo(2,1).NE.0).AND. & (cdcarbo(3,1).NE.0)) THEN ncdcase=3 GOTO 700 ELSE ncdcase=5 GOTO 700 ENDIF ENDIF ENDIF ENDDO ENDIF ENDDO ENDIF * check for C=C-C=O structure (case 2 or 4) DO i=1,4 IF (cdtable(i).ne.0) THEN DO j=1,mca IF (bond(cdtable(i),j).eq.1) THEN IF (group(j)(1:3).EQ.'CHO') THEN IF (ncd.eq.2) THEN ncdcase=2 GOTO 700 ELSE IF (ncd.eq.4) THEN ncdcase=4 GOTO 700 ENDIF ENDIF IF (group(j)(1:2).EQ.'CO') THEN IF (ncd.eq.2) THEN ncdcase=2 GOTO 700 ELSE IF (ncd.eq.4) THEN ncdcase=4 GOTO 700 ENDIF ENDIF ENDIF ENDDO ENDIF ENDDO * otherwise case 1 (if this point is reached) ncdcase=1 700 CONTINUE IF (wtflag.NE.0) WRITE(*,*) progname, ': end : case =',ncdcase END