!********************************************************** ! SUBROUTINE : gromhe ! ! Purpose : This subroutine returns the effective Henry's law ! coefficient (Keff) for the species given as input (chem) ! ! The gromhe method to compute the henry's law coefficient ! is described in Raventos-Duran et al., ACPD, 2010 ! ! INPUT : ! - chem : formula of the input species ! - Keff : Intrinsic Henry's law coefficient ! !!!!!!!!! SUBROUTINE gromhe(chem,Keff,nhenrydat,henrydat,chemhenry, & nhydratdat,hydratdat,chemhydrat) IMPLICIT NONE INCLUDE 'general.h' ! INPUT CHARACTER(LEN=lfo) chem, chemhydrat(mrd) INTEGER nhydratdat REAL hydratdat(mrd) ! database INTEGER nhenrydat CHARACTER(LEN=lfo) chemhenry(mrd) REAL henrydat(mrd,3) ! OUTPUT REAL keff ! INTERNAL CHARACTER(LEN=lfo) tchem ! group contributors : REAL sigma INTEGER onitrofol,haloica,nogrp REAL caoxa,caoxb INTEGER noh16,noh15,nfcd,nfaro INTEGER ic,ih,in,io,ir,is,ifl,ibr,icl INTEGER nc,i,j,k,l,ipos,it,i1,i2,ncd,ngrp,inte,ierr CHARACTER(LEN=lgr) group(mca) INTEGER bond(mca,mca),node,dbflg,nring INTEGER rjg(mri,2) ! ring-join group pairs REAL kest ! estimated intrinsic Henry's law coef. CHARACTER*1 nodetype(mca) REAL alifun(20),cdfun(20),arofun(20) INTEGER funflg(mca) INTEGER tabester(4,2) ! 1= -O- side, 2= CO side REAL mapfun(mca,3,20) REAL alisig(20),sigester ! taft sigmas REAL mult,tsig,dist INTEGER nabcde(9), tabcde(9,mco,mca) ! provided for hydration INTEGER nwa CHARACTER(LEN=lfo) chemhyd(mhyd,mhiso) REAL yhyd(mhyd,mhiso) INTEGER nhyd(mhyd) REAL khydstar ! Tafta sigma values, alifatic compounds c DATA funnam/'ROH ','RNO2 ','RONO2','ROOH ','RF ', c & 'RCl ','RBr ','RI ','RCHO ','RCOR ', c & 'RCOOH','COOOH','PAN ','ROR ','RCOOR', c & 'HCOOR','RCOF ','RCOCl','RCOBr','RCOI '/ DATA alisig/ 0.62 , 1.47 , 1.38 , 0.62 , 1.10 , & 0.94 , 1.00 , 1.00 , 2.15 , 1.81 , & 2.08 , 2.08 , 2.00 , 1.81 , 2.56 , & 2.90 , 2.44 , 2.37 , 2.37 , 2.37 / ! sigma for ester depnd whether ester is connected from the -CO- side ! or the -O- side. Sigma=2.56 for -O- side and 2.00 for CO side DATA sigester / 2.00/ ! sigma for ester, CO side !hammet orto meta para ! DATA arometa/ 0.13 , 0.74 , 0.55 , 0.00 , 0.34 , ! & 0.37 , 0.39 , 0.35 , 0.36 , 0.36 , ! & 0.35 , 0.00 , 0.00 , 0.11 , 0.32 , ! & 0.00 , 0.55 , 0.53 , 0.53 , 0.53 / ! DATA aropara/ -0.38 , 0.78 , 0.70 , 0.00 , 0.06 , ! & 0.24 , 0.22 , 0.21 , 0.44 , 0.47 , ! & 0.44 , 0.00 , 0.00 , -0.28 , 0.39 , ! & 0.00 , 0.70 , 0.69 , 0.69 , 0.69 / ! benzoic as a reference for substituents in ortho ! DATA arorto/ 1.22 , 1.99 , 0.00 , 0.00 , 0.93 , ! & 1.28 , 1.35 , 1.34 , 0.00 , 0.07 , ! & 0.95 , 0.00 , 0.00 , 0.12 , 0.63 , ! & 0.00 , 0.00 , 0.00 , 0.00 , 0.00 / ! tchem(1:lfo)=chem(1:lfo) IF (chem(1:3).EQ.'#mm'.OR.chem(1:3).EQ.'#bb' ) THEN chem(1:)=tchem(4:) ENDIF !------------------------------------------------------------- !Check if henry's law constant is already known in the database !------------------------------------------------------------- DO i=1,nhenrydat IF (chem .eq. chemhenry(i)) THEN Keff = henrydat(i,1) GOTO 999 ENDIF ENDDO ! --------- ! Initialize ! --------- ! build the group and bond matrix for chem nc=index(chem,' ') CALL grbond(chem,nc,group,bond,dbflg,nring) DO i=1,mca bond(i,i)=0 ENDDO CALL rjgrm(nring,group,rjg) ! rm ring index and get ring closing nodes ! count the number of group in chem node=0 DO i=1,mca IF (group(i)(1:1).NE.' ') node = node + 1 ENDDO ierr = 0 ! ---------------------- ! Get hydration constant ! ---------------------- CALL get_hydrate(chem,nwa,nhyd,chemhyd,yhyd,khydstar, & nhydratdat,hydratdat,chemhydrat) ! --------------------------------- ! Count functionalities ! --------------------------------- ! get the number of atom - C and H counters ***MC CALL number(chem,nc,ic,ih,in,io,ir,is,ifl,ibr,icl) ***MC ! table of functions - Index of functionalities ! -------------------------- ! 1= -OH ; 2= -NO2 ; 3= -ONO2 ; 4= -OOH ; 5= -F ! 6= -Cl ; 7= -Br ; 8= -I ; 9= -CHO ; 10= -CO- ! 11= -COOH ; 12= -CO(OOH) ; 13= -PAN ; 14= -O- ; 15= R-COO-R ! 16 = HCO-O-R; 17= -CO(F) ; 18= -CO(Cl) ; 19= -CO(Br) ; 20= -CO(I) ! write(92,*) chem(1:80),ierr CALL chemmap(chem,node,group,bond,ngrp,nodetype, & alifun,cdfun,arofun,mapfun,funflg, & tabester,nfcd,nfaro,ierr) IF (ierr.ne.0) THEN WRITE(6,*) 'STOP in chemmap, using chem:' WRITE(6,*) chem(1:60) c STOP ! write(92,*) 'err',chem(1:80),ierr ENDIF ! --------------------------------- ! compute taft sigmas ! --------------------------------- sigma=0. ! 1-1 interaction ! ---------------- DO i=1,node IF (funflg(i).gt.1) THEN DO k=1,8 DO it=1,2 ! check for both saturated and unsaturated aliphatic ! self interaction IF (mapfun(i,it,k).ne.0) THEN IF (mapfun(i,it,k).gt.1) THEN write(98,*) 'alisig(k)=',alisig(k) sigma=sigma+(alisig(k)/0.4)* & mapfun(i,it,k)*(mapfun(i,it,k)-1) ENDIF ENDIF ! cross interaction DO l=1,8 IF ((mapfun(i,it,k).gt.0.).and.(l.ne.k)) THEN sigma=sigma+(alisig(k)/0.4)*mapfun(i,it,k)* & mapfun(i,it,l) ENDIF ENDDO ENDDO ENDDO ENDIF ENDDO ! end do loop for the 1-1 interaction ! multiple node interaction ! ------------------------- DO i=1,node ! start loop among node IF (funflg(i).ne.0) THEN ! manage node i CALL abcde_map(bond,i,node,nabcde,tabcde) ! ester functionality is on 2 nodes - remove the path that would lead ! to count some interaction twice. This is performed by "killing" the ! "second" node of the functionality. DO k=1,4 IF (tabester(k,1).ne.0) THEN DO ipos=2,9 DO it=1,nabcde(ipos) DO j=1,ipos-1 IF (tabcde(ipos,it,j).eq.tabester(k,1)) THEN IF (tabcde(ipos,it,j+1).eq.tabester(k,2)) THEN tabcde(ipos,it,j+1)=0 ! set the node to 0 ENDIF ENDIF IF (tabcde(ipos,it,j).eq.tabester(k,2)) THEN IF (tabcde(ipos,it,j+1).eq.tabester(k,1)) THEN tabcde(ipos,it,j+1)=0 ! set the node to 0 ENDIF ENDIF ENDDO ENDDO ENDDO ENDIF ENDDO ! end of loop to kill irrelevant ester nodes ! start loop for sigma effect DO j=1,20 IF ( (mapfun(i,1,j).ne.0).OR. & (mapfun(i,2,j).ne.0) ) THEN !start sigma - taft DO inte=2,9 DO k=1,nabcde(inte) ncd=0. IF (tabcde(inte,k,inte).eq. 0) CYCLE IF (funflg(tabcde(inte,k,inte)).eq. 0) CYCLE !CMV! rule out case when tabcde(int,k,int) .eq. 0 because !CMV! associated funflg(0) (invisible seg fault!) !CMV! can be .ne. 0, depending on memory issues DO l=1,inte-1 ! start the loop to see if Cd are in between i1=tabcde(inte,k,l) i2=tabcde(inte,k,l+1) IF (bond(i1,i2).eq.2) ncd=ncd+1 ENDDO mult=real(funflg(tabcde(inte,k,inte))) IF (tabcde(inte,k,2).eq.0) mult=0. ! kill second node on ester tsig=alisig(j) IF (j.eq.15) THEN ! check ester IF (group(i)(1:2).eq.'CO') tsig=sigester ! CO instead of -O- ENDIF dist=REAL(inte)-REAL(ncd)-2. IF (group(tabcde(inte,k,inte))(1:3).EQ.'-O-') THEN ! decrease the distance dist=dist-1. ENDIF IF (nodetype(i).eq.'o') dist=dist-1. sigma=sigma+(tsig*0.4**(dist)*mult) ENDDO ENDDO ENDIF ! end sigma - taft ENDDO ! end loop for sigma effect ENDIF ! manage node i ENDDO ! end loop among node ! lump group on Cd and saturated aliphatic DO i=1,20 IF (cdfun(i).gt.0) THEN alifun(i)=alifun(i)+cdfun(i) ENDIF ENDDO ! lump group on aromatic and saturated aliphatic DO i=1,20 IF (arofun(i).gt.0) THEN alifun(i)=alifun(i)+arofun(i) ENDIF ENDDO ! ------------------------------------- ! compute other group-group descriptor ! ------------------------------------- ! compute correcting factors nogrp=1 IF (ngrp.ne.0) nogrp=0 ! check for ortho nitro phenols CALL nitrofol(node,group,bond,onitrofol) ! check for halogen next to a carboxylic group CALL haloic(node,group,bond,haloica) ! check for H bounding from an alcohol noh15=0 noh16=0 IF(alifun(1).gt.0) CALL hydol(node,group,bond,noh15,noh16) ! check for -CO-C(X) and -CO-C-C(X) structure CALL caox(node,group,bond,caoxa,caoxb) ! species having more than 1 ring where excluded during the ! development of gromhe ! IF (nring.le.1) THEN ! nring=0 ! ELSE ! nring=1 ! ENDIF ! IF (nring.ge.1) GOTO 999 ! species having no H atmo where excluded during the ! development of Grohme ! IF (ih.EQ.0) GOTO 999 cccccccccccccc ALL DATA kest = -1.51807 & - 0.13764*sigma & - 2.66308*onitrofol & + 0.973418*haloica & - 1.76731*caoxa & - 1.09538*caoxb & - 1.02636*noh16 & - 0.596601*noh15 & + 0.499757*ic & - 0.31011*ih & - 0.27637*nogrp & - 0.524482*nfcd & - 1.12707*nfaro & + 4.56593*alifun(1) & + 3.01781*alifun(2) & + 2.03719*alifun(3) & + 4.8687*alifun(4) & + 0.59747*alifun(5) & + 0.870225*alifun(6) & + 1.05894*alifun(7) & + 1.22831*alifun(8) & + 2.59268*alifun(9) & + 3.16535*alifun(10) & + 5.09321*alifun(11) & + 4.683*alifun(12) & + 1.93179*alifun(13) & + 2.4008*alifun(14) & + 2.78457*alifun(15) & + 2.35694*alifun(16) ! keff=10**kest ! keff=keff*(1+khydstar) ! keff=log10(keff) ! Keff=kest! + log10(1+khydstar) Keff=kest + log10(1+khydstar) 999 CONTINUE chem(1:)=tchem(1:) END