***************************************************************************
* SUBROUTINE DELOC                                                        *
*                                                                         *
* this subroutine checks if there is a delocalisation is possible in the  *
* case of a C=C double bond in alpha position of the alkyl radical :      *
*   -C.-C=C<  <->  -C=C-C.<                                               *
* the program produces the most stable radical. If the two possibles      *
* radicals have the same substitution degree, the output is a table with  *
* the differents formula.                                                 *
*                                                                         *
* INPUT  : -chem                                                          *
*                                                                         *
* OUTPUT : -tchem(mca) : formula of the different species.                *
*          -nchem : number of different species.                          *
*                   (currently max = 2)                                   *
*                                                                         *
***************************************************************************

      SUBROUTINE deloc(chem,tchem,nchem,sc)
      IMPLICIT NONE
      INCLUDE 'general.h'
      INCLUDE 'common.h'

* input
      CHARACTER(LEN=lfo)  chem

* output
      CHARACTER(LEN=lfo)  tchem(mca),temp
      INTEGER          nchem

*internal
      INTEGER          i,j,k,n1,n3,nring,dbflg
      INTEGER          j1,j2,j3,nc,nca
      CHARACTER(LEN=lgr)  group(mca)
      INTEGER          bond(mca,mca)
      INTEGER          rjg(mri,2)
      REAL             sc(mca),sctmp

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

      IF (wtflag.NE.0) WRITE(6,*)progname,': input : ',chem
      !IF (wtflag.NE.0) WRITE(6,*)progname

!-------------
* initialize:
!-------------
      nchem = 0
      j1=0
      j2=0
      j3=0
      n1 = 0
      n3 = 0
      nca = 0
      tchem = ' '
      tchem(1) = chem
      nchem = 2
      sc = 0.
      sc(1) = 1.


* find groups and bond matrix from chem      

      nc = INDEX(chem,' ') - 1
      CALL grbond(chem,nc,group,bond,dbflg,nring)
      IF (nring.GT.0) CALL rjgrm(nring,group,rjg)

      DO i=1,mca
        IF (group(i)(1:1).NE.' ') THEN
          nca=nca+1
        ENDIF
      ENDDO

* j1 = radical group (if none, escape)

      DO i=1,nca
        IF (INDEX(group(i),'.').NE.0) j1 = i
      ENDDO
      IF (j1.EQ.0) RETURN

* only treat uncomplicated radical groups
! JMLT: allow CH(OH). structure (but not any other ').' structures)

      IF ((INDEX(group(j1),'.').NE.0).AND.
     &    (INDEX(group(j1),'O.').EQ.0).AND.
     &    (INDEX(group(j1),'Cd').EQ.0))  THEN

       IF (INDEX(group(j1),').').EQ.0) THEN
c       IF((INDEX(group(j1),').').EQ.0).OR.
c     &    (INDEX(group(j1),'OH).').NE.0)) THEN

* j2 = double-bonded group attached to radical group (if none, escape)

        DO i=1,nca
          IF ((bond(i,j1).NE.0).AND.group(i)(1:2).EQ.'Cd') j2 = i
        ENDDO
        IF (j2.EQ.0) RETURN

* j3 = group at other end of double bond

        DO i=1,nca
          IF (bond(i,j2).EQ.2) j3=i 
        ENDDO

* n1 = number of bonds to radical group j1
* n3 = number of bonds to double-bonded group j3

        DO i=1,nca
          IF (bond(i,j1).NE.0) n1=n1+1 
          IF (bond(i,j3).NE.0) n3=n3+1 
        ENDDO

* move double bond one carbon over, to radical center 

        !bond(j2,j3)=1
        !bond(j3,j2)=1
        !bond(j1,j2)=2
        !bond(j2,j1)=2

        CALL setbond(bond,j1,j2,2)
        CALL setbond(bond,j2,j3,1)

* convert radical group j1 to double-bonded group
        
	IF (group(j1)(1:2).EQ.'C.')   group(j1)(1:2)='Cd'
        IF (group(j1)(1:3).EQ.'CH.')  group(j1)(1:3)='CdH'
        IF (group(j1)(1:4).EQ.'CH2.') group(j1)(1:4)='CdH2'
! JMLT: allow CH(OH). & C(OH). structures
        IF (group(j1)(1:6).EQ.'C(OH).') group(j1)(1:6)='Cd(OH)'
        IF (group(j1)(1:7).EQ.'CH(OH).') group(j1)(1:7)='CdH(OH)'

* convert 'double' group j3 to radical group

        IF (group(j3)(1:4).EQ.'CdH2') THEN
          group(j3)(1:4)='CH2.'
        ELSE IF (group(j3)(1:3).EQ.'CdH') THEN
          group(j3)(1:3)='CH.'
        ELSE IF (group(j3)(1:2).EQ.'Cd') THEN
          group(j3)(1:2)='C.'
        ENDIF

* write (non-standard) formula for new product

        CALL rebond(bond,group,tchem(2),nring)

* prioritise possible products
        IF (n1.GT.n3) THEN
	  sc(1) = 1.
	  sc(2) = 0.
c          nchem=1
        ELSE IF (n3.GT.n1) THEN
	  sc(1) = 0.
	  sc(2) = 1.
c          nchem=1
        ELSE IF (n3.EQ.n1) THEN
	  sc(1) = 0.5
	  sc(2) = 0.5
c          nchem=2
        ENDIF

! JMLT ! allow OH to attract electron density, as per Paulot ('09)
        IF(INDEX(group(j1),'OH').NE.0
     &    .OR.INDEX(group(j3),'OH').NE.0) THEN
          print*,'!Paulot case! ',group(j1),group(j3)
          sc(1) = 0.65
          sc(2) = 0.35
        ENDIF
! JMLT ! end additional section

        IF (sc(1).LT.sc(2)) THEN
          temp=tchem(1)
          sctmp=sc(1)
          tchem(1)=tchem(2)
          sc(1)=sc(2)
          tchem(2)=temp
          sc(2)=sctmp
        ENDIF

        WRITE(42,*) '------------'
        WRITE(42,*) 'deloc done:'
        WRITE(42,*) 'chem :',chem(1:50)	       
        WRITE(42,*) 'first :',sc(1),'  ',tchem(1)(1:50)	       
        WRITE(42,*) 'second:',sc(2),'  ',tchem(2)(1:50)	       

        IF (wtflag.NE.0) THEN
          print*,progname,': output (1) :',sc(1),tchem(1)
          print*,progname,': output (2) :',sc(2),tchem(2)
        ENDIF

       ENDIF
      ENDIF

      !IF (wtflag.NE.0) print*,progname,': end : '
      END