!*********************************************************
! SUBROUTINE  : get_hydrate                              
!                                                        
! PURPOSE     : get all hydrates that can be made from a
!               given formula and related hydration constants
! INPUT :
! - chem : formula of the input species
! 
! OUTPUT :
! - nwa : number of water molecule that can be added to chem.
!         (e.g. if the molecule bear 3 ketones, the nwa is return as 3).
! - nhyd(i) : number of distinct hydrate that can be made for i molecule
!             of water added.
! - chemhyd(j,i) : formula of the jth hydrate bearing ith added water 
!                  molecule
! - yhyd(j,i) : hydration constant (with respect to the fully 
!               non-hydrated molecule) of the jth hydrate bearing 
!               ith added water molecules (log scale)
! - khydstar : hydration constant taking all possible hydrate into
!              account (apparent constant K_star). NOT A LOG SCALE !
!
! INTERNAL 
! - ndat : number of distinct molecule (hydrate) that can be made
!          from chem (again 1 water molecule added only !)
! - chemdat(i) : formula of the i'th hydrate
! - ydat(i) : hydration constant of the i'th hydrate (with respect 
!             to the fully non hydrated molecule !)
!
!*********************************************************
      SUBROUTINE get_hydrate(chem,nwa,nhyd,chemhyd,yhyd,khydstar,
     &          nhydratdat,hydratdat,chemhydrat)
      IMPLICIT NONE
      INCLUDE 'general.h'

! input
      CHARACTER(LEN=lfo) chem, chemhydrat(mrd)
      INTEGER         nhydratdat
      REAL            hydratdat(mrd)
      
! output
      INTEGER         nwa
      CHARACTER(LEN=lfo) chemhyd(mhyd,mhiso)
      REAL            yhyd(mhyd,mhiso)
      INTEGER         nhyd(mhyd)
      REAL            khydstar

! one dimension local table
      INTEGER  nc,i,j,k,l
      CHARACTER(LEN=lfo) t1chemhyd(mhiso)
      REAL            t1yhyd(mhiso)
      INTEGER         t1nhyd

      REAL            yield,numy,sumy
      CHARACTER(LEN=lfo) chemdat(mhyd)
      REAL            ydat(mhyd)
      INTEGER         ndat
      INTEGER         jp
      CHARACTER(LEN=lfo) tempkc
      INTEGER         hlev
      

! 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)

! -----------
! initialize
! -----------
      nwa=0
      khydstar=0.
      DO i=1,mhyd
        nhyd(i)=0
        DO j=1,mhiso
          chemhyd(i,j)=' '
          yhyd(i,j)=0.
        ENDDO
      ENDDO

      ndat = 0
      ydat(:) = 0.
      chemdat(:)=' '

! -----------
! monohydrate
! -----------
      t1nhyd=0
      DO i=1,mhiso
        t1chemhyd(i)=' '
        t1yhyd(i)=0.
      ENDDO

      yield=0. ! parent compound
      CALL khydration(chem,yield,ndat,chemdat,ydat)
      IF (ndat.eq.0) RETURN   ! no hydrate can be made from chem
      t1nhyd=ndat
      DO i=1,ndat
        t1chemhyd(i)=chemdat(i)
        t1yhyd(i)=ydat(i)
      ENDDO

! collapse identical formula
      DO i=1,t1nhyd-1
        numy=1
        sumy=t1yhyd(i)
        DO j=i+1,t1nhyd
          IF (t1chemhyd(i).eq.t1chemhyd(j)) THEN
            numy=numy+1
            sumy=sumy+t1yhyd(j)
            t1chemhyd(j)=' '
            t1yhyd(j)=0.
          ENDIF
        ENDDO
        IF (numy.gt.1) THEN
          t1yhyd(i)=sumy/real(numy)
        ENDIF
      ENDDO      

! write the table for monohydrate
      DO i=1,t1nhyd
        IF (t1chemhyd(i)(1:1).ne.' ') THEN
          nhyd(1)=nhyd(1)+1
          chemhyd(1,nhyd(1))=t1chemhyd(i)
          yhyd(1,nhyd(1))=t1yhyd(i)
        ENDIF
      ENDDO

! -----------
! multi-hydrate
! -----------
      DO hlev=2,mhyd
        IF (nhyd(hlev-1).eq.0) THEN  ! all hydrate were found
          nwa=hlev-1
          GOTO 300
        ENDIF

        t1nhyd=0
        DO i=1,mhiso
          t1chemhyd(i)=' '
          t1yhyd(i)=0.
        ENDDO
      
        DO jp=1,nhyd(hlev-1)
          yield=yhyd(hlev-1,jp) ! parent compound
          tempkc=chemhyd(hlev-1,jp)
          CALL khydration(tempkc,yield,ndat,chemdat,ydat)
          DO i=1,ndat
            t1chemhyd(t1nhyd+i)=chemdat(i)
            t1yhyd(t1nhyd+i)=ydat(i)
          ENDDO
          t1nhyd=t1nhyd+ndat
        ENDDO

! collapse identical formula
        DO i=1,t1nhyd-1
          numy=1
          sumy=t1yhyd(i)
          DO j=i+1,t1nhyd
            IF (t1chemhyd(i).eq.t1chemhyd(j)) THEN
              numy=numy+1
              sumy=sumy+t1yhyd(j)
              t1chemhyd(j)=' '
              t1yhyd(j)=0.
            ENDIF
          ENDDO
          IF (numy.gt.1) THEN
            t1yhyd(i)=sumy/real(numy)
          ENDIF
        ENDDO      

! write the table for monohydrate
        DO i=1,t1nhyd
          IF (t1chemhyd(i)(1:1).ne.' ') THEN
            nhyd(hlev)=nhyd(hlev)+1
            chemhyd(hlev,nhyd(hlev))=t1chemhyd(i)
            yhyd(hlev,nhyd(hlev))=t1yhyd(i)
          ENDIF
        ENDDO

      ENDDO

300   CONTINUE
! exit of the DO loop usually occurs for the formula. Must
! adjust nwa if exit did not occur.
      IF (hlev.ge.mhyd) nwa=mhyd   

! compute Kstar
      DO i=1,nwa
        DO j=1,nhyd(i)
          khydstar=khydstar+10**(yhyd(i,j))
        ENDDO
      ENDDO

!-------------------------------------------------------------
!Check if hydration constant is already known in the database
!-------------------------------------------------------------
      DO i=1,nhydratdat
        IF (chem .eq. chemhydrat(i)) THEN
          khydstar = 10**hydratdat(i)
        ENDIF
      ENDDO

      END