SUBROUTINE hoadd_arom(chem,bond,group,nca,
     &             nr,flag,tarrhc,pchem,coprod)
      IMPLICIT NONE
      INCLUDE 'general.h'
      INCLUDE 'common.h'

! input:
      CHARACTER(LEN=lfo),INTENT(in) :: chem
      INTEGER,INTENT(in)         :: bond(mca,mca)
      CHARACTER(LEN=lgr),INTENT(in) :: group(mca)
      INTEGER,INTENT(in)         :: nca

! input/output
      INTEGER,INTENT(inout)      :: nr, flag(mnr)
      REAL,INTENT(inout)            :: tarrhc(mnr,3)
      CHARACTER(LEN=lfo),INTENT(inout) :: pchem(mnr)
      CHARACTER(LEN=lco),INTENT(inout) :: coprod(mnr,mca)

! internal
      INTEGER                    :: i,j,k,li,num_sub
      REAL                       :: karom(3),kipso(3)
      REAL                       :: kabs(3),kadd1(3),kadd2(3)
      REAL                       :: kabs_298,kadd1_298,kadd2_298
      INTEGER                    :: o_sub(2),m_sub(2),p_sub(2)
      INTEGER                    :: path(6),nring
      REAL                       :: frac(3),sum298,arrhc(3)
      INTEGER                    :: nb_alp,alp(4)
      INTEGER                    :: track(mco,mca)
      INTEGER                    :: trlen(mco)
      INTEGER                    :: ntr,nc

      CHARACTER(LEN=lfo) rdckprod(mca),pchem_del(mnr),tempkc
      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)

      INTEGER         :: tbond(mca,mca)
      CHARACTER(LEN=lgr) :: tgroup(mca)
      CHARACTER(LEN=lgr) :: pold,pnew

      INTEGER    rngflg       ! 0 = 'no', 1 = 'yes'
      INTEGER    ring(mca)    ! =1 if node participates in current ring

      CHARACTER(LEN=lsb) :: progname='*hoadd_arom*'
      CHARACTER(LEN=ler) :: mesg

! write info for finding bugs
      IF (wtflag.NE.0) WRITE(*,*) progname,': input : ',chem

!!!!
      karom(1)=0.378E-12
      karom(2)=0
      karom(3)=190

      kipso(1)=0.378E-12
      kipso(2)=0
      kipso(3)=89

!      kabs(1)=0.378
!      kabs(2)=0
!      kabs(3)=190
!
!      kadd1(1)=0.378
!      kadd1(2)=0
!      kadd1(3)=190
!
!      kadd2(1)=0.378
!      kadd2(2)=0
!      kadd2(3)=190

       frac(:)=0
       arrhc(:)=0
!!!!

compute addition rate constant on each carbon of the aromatic ring      
      DO i=1,nca
        IF (INDEX(group(i),'O').NE.0) CYCLE
        IF (group(i)(1:1).EQ.'c') THEN
          CALL arom_data(i,group,bond,nca,o_sub,m_sub,p_sub)
          
          nb_alp=0
          alp(:)=0
          DO j=1,nca
            IF ((group(j)(1:1).EQ.'c').AND.(bond(i,j).EQ.1)) THEN
              nb_alp=nb_alp+1
              alp(nb_alp)=j
            ENDIF
          ENDDO
          IF (INDEX(chem,'c(OH)').NE.0) THEN
            IF ((group(alp(1)).NE.'c(OH)').AND.
     &          (group(alp(2)).NE.'c(OH)')) CYCLE
          ENDIF





!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!     COMPUTE RATE CONSTANTS     !!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! addition on a non substitued carbon in the aromatic ring
          IF (group(i)(1:3).EQ.'cH ') THEN
            arrhc(1)=karom(1)
            arrhc(3)=karom(3)
          ELSE
! addition on a substitued carbon in the aromatic ring
            arrhc(1)=kipso(1)
            arrhc(3)=kipso(3)
          ENDIF

          num_sub=0
          DO j=1,2
            IF (o_sub(j).NE.0) num_sub=num_sub+1
            IF (m_sub(j).NE.0) num_sub=num_sub+1
            IF (p_sub(j).NE.0) num_sub=num_sub+1
          ENDDO

          IF (num_sub.EQ.1) THEN
            IF (m_sub(1).NE.0) THEN
              arrhc(1)=arrhc(1)*0.7
              arrhc(3)=arrhc(3)-207
            ELSE
              arrhc(1)=arrhc(1)*0.8
              arrhc(3)=arrhc(3)-659
            ENDIF
          ELSE IF (num_sub.EQ.2) THEN
            IF ((m_sub(1).NE.0).AND.(m_sub(2).NE.0)) THEN
              arrhc(1)=arrhc(1)*2.6
              arrhc(3)=arrhc(3)-416
            ELSE IF (((o_sub(1).NE.0).AND.(o_sub(2).NE.0))
     &               .OR.((o_sub(1).NE.0).AND.(p_sub(1).NE.0))) THEN
              arrhc(1)=arrhc(1)*0.8
              arrhc(3)=arrhc(3)-659
            ELSE
              arrhc(1)=arrhc(1)*0.6
              arrhc(3)=arrhc(3)-1203
            ENDIF
          ELSE IF (num_sub.EQ.3) THEN
            IF ( (o_sub(1).NE.0).AND.(o_sub(2).NE.0).AND.
     &           (p_sub(1).NE.0) )THEN
              arrhc(1)=arrhc(1)*1.9
              arrhc(3)=arrhc(3)-409
            ELSE IF ( ((o_sub(1).NE.0).AND.(o_sub(2).NE.0).AND.
     &                 (m_sub(1).NE.0)) .OR. ((o_sub(1).NE.0).AND.
     &                 (p_sub(1).NE.0).AND.(m_sub(1).NE.0)) ) THEN
              arrhc(1)=arrhc(1)*6.8
              arrhc(3)=arrhc(3)-760
            ELSE
              arrhc(1)=arrhc(1)*0.5
              arrhc(3)=arrhc(3)-1200
            ENDIF
          ELSE IF (num_sub.EQ.4) THEN
            arrhc(1)=arrhc(1)*3.5
            arrhc(3)=arrhc(3)-341
          ELSE IF (num_sub.EQ.5) THEN
            arrhc(1)=arrhc(1)*2.0
            arrhc(3)=arrhc(3)-998
          ENDIF

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!            REACTIONS           !!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! to this point, OH has been added on the ring, forming a delocalized
! radical. This radical can follow three reactions with O2 :
! 1- H abstraction if OH was added on a cH group
! 2 and 3 : addition of O2 in alpha position of the OH bearing carbon
          
! CASE 1
          kabs(:)=0.
          IF (group(i)(1:3).EQ.'cH ') THEN
            kabs(1)=1.75E-10
            kabs(2)=0
            kabs(3)=1500
          ENDIF
! CASE 2-3: look for neighbours where O2 will be added
!          nb_alp=0
!          alp(:)=0
!          DO j=1,nca
!            IF ((group(j)(1:1).EQ.'c').AND.(bond(i,j).EQ.1)) THEN
!              nb_alp=nb_alp+1
!              alp(nb_alp)=j
!            ENDIF
!          ENDDO

! CASE 2:
          CALL gettrack(bond,alp(1),nca,ntr,track,trlen)
          CALL arom_path(i,ntr,track,group,bond,path)

          kadd1(1)=1.50E-10
          kadd1(2)=0
          kadd1(3)=1700
          IF (group(path(2))(1:3).NE.'cH ') kadd1(3)=kadd1(3)-207
          IF (group(path(3))(1:3).NE.'cH ') kadd1(3)=kadd1(3)-620
          IF (group(path(4))(1:3).NE.'cH ') kadd1(3)=kadd1(3)-207
          IF (group(path(5))(1:3).NE.'cH ') kadd1(3)=kadd1(3)-558

! CASE 3:
          CALL gettrack(bond,alp(2),nca,ntr,track,trlen)
          CALL arom_path(i,ntr,track,group,bond,path)

          kadd2(1)=1.50E-10
          kadd2(2)=0
          kadd2(3)=1700
          IF (group(path(2))(1:3).NE.'cH ') kadd2(3)=kadd2(3)-207
          IF (group(path(3))(1:3).NE.'cH ') kadd2(3)=kadd2(3)-620
          IF (group(path(4))(1:3).NE.'cH ') kadd2(3)=kadd2(3)-207
          IF (group(path(5))(1:3).NE.'cH ') kadd2(3)=kadd2(3)-558

! compute fraction for three cases
          kabs_298=kabs(1)*EXP(-kabs(3)/298)
          kadd1_298=kadd1(1)*EXP(-kadd1(3)/298)
          kadd2_298=kadd2(1)*EXP(-kadd2(3)/298)

          sum298=kabs_298+kadd1_298+kadd2_298

          frac(1)=kabs_298/sum298
          frac(2)=kadd1_298/sum298
          frac(3)=kadd2_298/sum298

!          arrhc(1:3)=tarrhc(nr,1:3)

! WRITE NEW SPECIES
! CASE 1  
          tgroup(:)=group(:)
          tbond(:,:)=bond(:,:)

          IF (tgroup(i)(1:3).EQ.'cH ') THEN
            nr=nr+1
            flag(nr)=1
            tarrhc(nr,1)=arrhc(1)*frac(1)
            tarrhc(nr,3)=arrhc(3)

            CALL gettrack(bond,alp(1),nca,ntr,track,trlen)
            CALL arom_path(i,ntr,track,group,bond,path)
            tgroup(i)='c(OH)'
!            pold='c'
!            pnew='Cd'
!            DO j=1,5
!              CALL swap(group(path(j)),pold,tgroup(path(j)),pnew)
!            ENDDO
!            CALL setbond(tbond,i,path(1),2)
!            CALL setbond(tbond,path(1),path(2),1)
!            CALL setbond(tbond,path(2),path(3),2)
!            CALL setbond(tbond,path(3),path(4),1)
!            CALL setbond(tbond,path(4),path(5),2)
!            CALL setbond(tbond,path(5),i,1)

            CALL rebond(tbond,tgroup,pchem(nr),nring)
            CALL stdchm(pchem(nr))

            tgroup(:)=group(:)
            tbond(:,:)=bond(:,:)
          ENDIF

! CASE 2

          DO k=1,2

! if there is an alcohol moeity on the aromatic ring, only addition OH
! in alpha position of the alcohol, (from phenol OH oxidation, Xu and
! Wang., 2013)
            IF (INDEX(chem,'c(OH)').NE.0) THEN     
              IF (group(alp(k))(1:5).NE.'c(OH)') CYCLE
            ENDIF

            nr=nr+1
            flag(nr)=1
            tarrhc(nr,1)=arrhc(1)*frac(k+1)
            tarrhc(nr,3)=arrhc(3)
            IF (tgroup(i)(1:3).EQ.'cH ') THEN 
              tgroup(i)='CH(OH)'
            ELSE
!              tgroup(i)='C(OH)'
              nc=INDEX(tgroup(i),' ')
              tgroup(i)(nc:nc+3)='(OH)'
              tgroup(i)(1:1)='C'
            ENDIF
            CALL gettrack(bond,alp(k),nca,ntr,track,trlen)
            CALL arom_path(i,ntr,track,group,bond,path)
!            IF (tgroup(path(1))(1:3).EQ.'cH ') THEN
!              tgroup(path(1))='CH(OO.) '
!            ELSE
!              tgroup(path(1))='C(OO.) '
!            ENDIF
            nc=INDEX(tgroup(path(1)),' ')
            tgroup(path(1))(nc:nc+4)='(OO.)'
            tgroup(path(1))(1:1)='C'
            pold='c'
            pnew='Cd'
            DO j=2,5
              CALL swap(group(path(j)),pold,tgroup(path(j)),pnew)
            ENDDO
 
            CALL setbond(tbond,i,path(1),1)
            CALL setbond(tbond,path(1),path(2),1)
            CALL setbond(tbond,path(2),path(3),2)
            CALL setbond(tbond,path(3),path(4),1)
            CALL setbond(tbond,path(4),path(5),2)
            CALL setbond(tbond,path(5),i,1)
            tempkc=' '
            CALL rebond(tbond,tgroup,tempkc,nring)
            CALL radchk(tempkc,rdckprod,rdcktprod,nip,sc)
            pchem(nr) = rdckprod(1)
            CALL stdchm(pchem(nr))
            
            tgroup(:)=group(:)
            tbond(:,:)=bond(:,:)

          ENDDO

        ENDIF
      ENDDO









      END