!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MASTER MECHANISM - ROUTINE NAME : isom_rate                             !
!                                                                         !
!                                                                         !
! PURPOSE :                                                               !
!  Find rate constants for alkoxy radical isomerization based on          !
!  Vereecken and Peeters 2010.                                            !
!                                                                         !
! INPUT:                                                                  !
! - tgroup(i)     : groups at position (carbon) i                         !
! - tbond(i,j)    : carbon-carbon bond matrix of chem                     !
! - ig            : group bearing the leaving H                           !
! - span          : it is the span of the migration :                     !
!                   - 4 for a 1-4 migration                               !
!                   - 5 for a 1-5 migration                               !
!                   - 6 for a 1-6 migration                               !
!                   - 7 for a 1-7 migration                               !
!                                                                         !
! - Fspan,tunnel = a_span,tunnel * exp(b_span_tunnel/T)                   !
! - Fsubs,tunnel = a_subs,tunnel * exp(b_subs_tunnel/T)                   !
!                                                                         !
!                                                                         !
!                                                                         !
!                                                                         !
!                                                                         !
!                                                                         !
! OUTPUT                                                                  !
! - ar1, ar2, ar3 : arrhenius coefficients 
!                                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE isom_rate(bond,group,ig,span,ar1,ar2,ar3,FSD)
      IMPLICIT NONE
      INCLUDE 'general.h'

! input
      INTEGER,INTENT(in)          :: bond(mca,mca)
      CHARACTER(LEN=lgr),INTENT(in)  :: group(mca)
      INTEGER,INTENT(in)          :: ig,span

! output
      REAL,INTENT(out)            :: ar1,ar2,ar3,FSD(5)
! internal:
      INTEGER          :: tbond(mca,mca)
      CHARACTER(LEN=lgr)  :: tgroup(mca)
      INTEGER          :: i,j,k,irad,ipath,log_fg
      INTEGER          :: nca                  ! number of nodes
      INTEGER          :: nc
      INTEGER          :: nabcde(9), tabcde(9,mco,mca)
      REAL             :: ktemp, temp
! coefficient of the polynome
      REAL             ::  a1,b1,c1,d1,a2,b2,c2,d2

     
! -----------
! initialize
! -----------
      log_fg=1

      nca = 0
      tgroup = group
      tbond = bond
      FSD(1:4) = 0
      FSD(5) = 1
      ipath=1

      ktemp = 0
      temp = 298

      a1 = 0
      b1 = 1
      c1 = 0
      d1 = 1
      a2 = 0
      b2 = 1
      c2 = 0
      d2 = 1


! check that the span given as an input is authorized
      IF ((span.NE.4).AND.(span.NE.5).AND.
     &    (span.NE.6).AND.(span.NE.7).AND.(span.NE.8)) THEN
        WRITE(6,*) 'ERROR in isomrate.f :'
        WRITE(6,*) 'The span given as an input is not possible'
        STOP
      ENDIF

! count the number of nodes
      DO i=1,mca
        IF (tgroup(i)(1:1).EQ.' ') EXIT
        IF (tgroup(i)(1:1).NE.' ') nca=nca+1
        IF (INDEX(tgroup(i),'(O.)').NE.0) irad = i
        IF (log_fg.EQ.1) THEN 
          WRITE(46,*) i, tgroup(i)
        ENDIF
      ENDDO

! get the track between the O. and the leaving H
      CALL abcde_map(tbond,irad,nca,nabcde,tabcde)
! get the right pathway to reach the leaving H
      DO i=1,nabcde(span-1)
        IF (tabcde(span-1,i,span-1).EQ.ig) THEN
          ipath=i
          EXIT
        ENDIF
      ENDDO

      IF (log_fg.EQ.1) WRITE(46,*) 'span=',span
      IF (log_fg.EQ.1) WRITE(46,*) 'ig=',ig
      IF (log_fg.EQ.1) WRITE(46,*) 'irad=',irad

! reference rate coefficients and correction for F_span_tunnel:
      IF (tgroup(ig)(1:3).EQ.'CH3') THEN
        ar1 = 3*4E10
        ar2 = 0.
        ar3 = 3825
        IF (span.EQ.6) THEN
          ar1=ar1*1.13
          ar3=ar3+59
        ELSE IF (span.EQ.7) THEN
          ar1=ar1*1.09
          ar3=ar3+49
        ELSE IF (span.EQ.8) THEN
          ar1=ar1*0.85
          ar3=ar3-65
        ELSE IF (span.EQ.4) THEN
          ar1=ar1*0.19
          ar3=ar3-678
        ENDIF
      ELSE IF (tgroup(ig)(1:3).EQ.'CH2') THEN
        ar1 = 2*4E10
        ar2 = 0.
        ar3 = 3010
        IF (span.EQ.6) THEN
          ar1=ar1*1.26
          ar3=ar3+117
        ELSE IF (span.EQ.7) THEN
          ar1=ar1*1.24
          ar3=ar3+116
        ELSE IF (span.EQ.8) THEN
          ar1=ar1*0.94
          ar3=ar3-16
        ELSE IF (span.EQ.4) THEN
          ar1=ar1*0.28
          ar3=ar3-556
        ENDIF
      ELSE IF (tgroup(ig)(1:2).EQ.'CH') THEN
        ar1 = 4E10
        ar2 = 0.
        ar3 = 2440
        IF (span.EQ.6) THEN
          ar1=ar1*1.25
          ar3=ar3+123
        ELSE IF (span.EQ.7) THEN
          ar1=ar1*0.97
          ar3=ar3-9
        ELSE IF (span.EQ.8) THEN
          ar1=ar1*0.96
          ar3=ar3-15
        ELSE IF (span.EQ.4) THEN
          ar1=ar1*0.36
          ar3=ar3-473
        ENDIF
      ENDIF

! Cstrain and Fspan(SD) Correction factors
      IF (span.EQ.4) THEN
        ar3=ar3+5440
        a1 = 2.68E-3
        b1 = 1.46
        c1 = 2.4E-3
        d1 = 0.498
      ELSE IF (span.EQ.6) THEN
        ar3=ar3+136
        a1 = -4.88E-4
        b1 = 0.6
        c1 = -6.12E-4
        d1 = 0.854
      ELSE IF (span.EQ.7) THEN
        ar3=ar3+1243 
        a1 = -3.09E-4
        b1 = 0.37
        c1 = -2.51E-4
        d1 = 0.562
      ELSE IF (span.EQ.8) THEN
        ar3=ar3+2179
        a1 = -2.48E-4
        b1 = 0.19
        c1 = 3.54E-4
        d1 = 0.603
      ENDIF

! Fspan(SD) and Fsubst(SD) : these coefficent are expressed in the
! following form : F(SD)=Frig*Fconf = (aT + b)*(cT + d) = acT² + (ad+bc)T + bd
! which give for Fspan(SD)*Fsubst(SD) something like :
! aa'cc'T⁴ + [(ac(a'd'+b'c')+a'c'(ad+bc)]*T³ +
! [(ad+bc)(a'd'+b'c')+(acb'd')+(bda'c')]*T² +
! [b'd'(ad+bc)bd(a'd'+b'c')]*T + bdb'd'
!
! The coefficients of the polynome have to be given in output to compute
! the rate in the boxmodel, using the temperature of the simulations
! coefficients of the polynome will be stored in FSD(1) to FSD(5)

      IF (log_fg.EQ.1) WRITE(46,*) 'substition correction factors :'
      
! Substition correction factor
! If the leaving H is in an aldehyde
      IF (group(ig)(1:4).EQ.'CHO ') THEN
        IF (log_fg.EQ.1) WRITE(46,*) ' CHO '
        IF (span.EQ.5) THEN
          ar1=ar1*1.5
          ar3=ar3+327+223
          a2 = 7.74E-3
          b2 = 2.34
          c2 = -1.03E-3
          d2 = 2.72
        ELSE IF (span.EQ.6) THEN
          ar1=ar1*1.27
          ar3=ar3+574+134
          a2 = 4.49E-4
          b2 = 2.23
          c2 = -8.32E-4
          d2 = 0.802
        ELSE IF (span.EQ.7) THEN
          ar1=ar1*1.64
          ar3=ar3+287+267
          a2 = 1.55E-3
          b2 = 4.92
          c2 = -1.03E-3
          d2 = 2.72
        ENDIF
      ENDIF

! endo B oxo : i.e. R(O.)-...-CO-CH2-...
      IF (tgroup(tabcde(span-1,ipath,span-2)).EQ.'CO') THEN
        IF (log_fg.EQ.1) WRITE(46,*) ' endo B oxo '
        IF (tgroup(ig)(1:3).EQ.'CH3') THEN
          IF (span.EQ.5) THEN
            ar1=ar1*0.57
            ar3=ar3-237+2159
            b2=0.74
          ELSE IF (span.EQ.6) THEN
            ar1=ar1*0.83
            ar3=ar3-76+1283
            a2=1.82E-3
            b2=1.15
          ENDIF
        ELSEIF (tgroup(ig)(1:2).EQ.'CH') THEN
          IF (span.EQ.5) THEN
            ar1=ar1*0.69
            ar3=ar3-165+1651
            b2=0.74
          ELSE IF (span.EQ.6) THEN
            ar1=ar1*0.92
            ar3=ar3-26+1057
            a2=1.82E-3
            b2=1.15
          ENDIF
        ENDIF
      ENDIF
      
! search for generic endo -oxo : R(O.)-CX-..-CO-..-CX-CHX-...
      IF (span.GE.5) THEN
        DO i=2,span-3
          IF (tgroup(tabcde(span-1,ipath,i)).EQ.'CO') THEN
            IF (log_fg.EQ.1) WRITE(46,*) ' generic endo-oxo '
            IF (span.EQ.5) THEN
              ar1=ar1*0.92
              ar3=ar3-34+354
              a2=4.02E-3
              b2=1.1
            ELSE IF (span.EQ.6) THEN
              ar1=ar1*0.92
              ar3=ar3-34+851
              a2=4.02E-3
              b2=1.1
            ENDIF
          ENDIF
        ENDDO
      ENDIF

! and acyl-oxi R(O.)-CO-R-CHX- ... subst correction factor
      IF (tgroup(irad).EQ.'CO(O.)') THEN
        IF (log_fg.EQ.1) WRITE(46,*) ' acyloxy '
        IF (tgroup(ig)(1:3).EQ.'CH3') THEN
          IF (span.EQ.5) THEN
            ar1=ar1*1.26
            ar3=ar3+114+667
            a2=7.3E-4
            b2=0.13
          ELSE IF (span.EQ.6) THEN
            ar1=ar1*1.26+549
            ar3=ar3+114
            a2=7.3E-4
            b2=0.13
          ENDIF
        ELSEIF (tgroup(ig)(1:2).EQ.'CH') THEN
          IF (span.EQ.5) THEN
            ar1=ar1*1.43+667
            ar3=ar3+187
            a2=4.27E-4
            b2=6.03E-2
          ELSE IF (span.EQ.6) THEN
            ar1=ar1*1.43+549
            ar3=ar3+187
            a2=4.27E-4
            b2=6.03E-2
          ENDIF
        ENDIF
      ENDIF

! alpha-OH subst correction factor
      IF (INDEX(tgroup(ig),'(OH)').NE.0) THEN
        IF (log_fg.EQ.1) WRITE(46,*) ' alpha-OH '
        IF (span.EQ.5) THEN
          ar1=ar1*1.62
          ar3=ar3+258-310
          a2=-9.21E-3
          b2=6.00
        ELSE IF (span.EQ.6) THEN
          ar1 = ar1*1.62
          ar3 = ar3+258-539
          a2 = -3.98E-2
          b2 = 18.2
        ENDIF
      ENDIF

! R(OH)(O.) subst correction factor
      IF (INDEX(tgroup(irad),'(OH)').NE.0) THEN
        IF (log_fg.EQ.1) WRITE(46,*) ' RC(OH)(O.) '
        ar1 = ar1*1.17
        ar3 = ar3+68-956
        a2 = -5.61E-4
        b2 = 0.633
        c2 = 1.16E-4
        d2 = 0.313
      ENDIF

! endo beta OH
      IF (INDEX(tgroup(tabcde(span-1,ipath,span-2)),'(OH)').NE.0) THEN
        IF (log_fg.EQ.1) WRITE(46,*) ' endo B OH '
        IF (span.EQ.5) THEN
          ar1=ar1*0.89
          ar3=ar3-47+974
          a2 = -3.64E-3
          b2 = 6.70
        ELSE IF (span.EQ.6) THEN
          ar1=ar1*0.93
          ar3=ar3-35-30
        ENDIF
      ENDIF

! exo beta OH
      DO i=1,nca
        IF ((tbond(i,ig).EQ.1).AND.(i.NE.tabcde(span-1,ipath,span-2))
     &  .AND.(INDEX(tgroup(i),'(OH)').NE.0)) THEN
          IF (log_fg.EQ.1) WRITE(46,*) ' exo B OH '
          IF (span.EQ.5) THEN
            ar3=ar3-252
            a2 = -4.53E-3
            b2 = 2.12
          ELSE IF (span.EQ.6) THEN
            ar1=ar1*1.18
            ar3=ar3+83-392
            a2 = -4.53E-3
            b2 = 2.12
          ENDIF
          EXIT
        ENDIF
      ENDDO

!!!!!!!!!!!!!!!!!!

      FSD(1) = a1*a2*c1*c2
      FSD(2) = (a1*c1*(a2*d2+b2*c2)+a2*c2*(a1*d1+b1*c1))
      FSD(3) = ((a1*d1+b1*c1)*(a2*d2+b2*c2)+(a1*c1*b2*d2)+(a2*c2*b1*d1))
      FSD(4) = b2*d2*(a1*d1+b1*c1) + b1*d1*(a2*d2+b2*c2)
      FSD(5) = b1*d1*b2*d2

!! Write the parameters in the log file
      IF (log_fg.EQ.1) THEN 
        WRITE(46,*) ' --- '
        WRITE(46,*)  'Fspan rig = ',b1,'+',a1,'T'
        WRITE(46,*)  'Fspan conf = ',d1,'+',c1,'T'
        WRITE(46,*)  'Fsubs rig = ',b2,'+',a2,'T'
        WRITE(46,*)  'Fsubs conf = ',d2,'+',c2,'T'
        WRITE(46,*) 'FSD(1:5)',FSD
        WRITE(46,*) 'arrhenius coef',ar1,ar2,ar3
        ktemp=ar1 *(FSD(1)*temp**4 + FSD(2)*temp**3 + FSD(3)*temp**2 + 
     &            FSD(4)*temp + FSD(5))* temp**ar2 * exp(-ar3/temp)
        WRITE(46,*) 'k_isom vereecken à ',temp,'K = ',ktemp
        WRITE(46,*) ' '
        WRITE(46,*) ' -------------------------------------- '
      ENDIF

      END