      subroutine hydrat(did,iclust,idim,iwat,ihydtp,iat,ihyd,itot,
     &           elemt,watat,ix)            
c ----------------------------------------------------------------------HYDR002
c     This program module is part of the author's X-ray Library         HYDR003
c     and may not be sold or reproduced for commercial purposes.        HYDR004
c     Although it has been tested carefully, there is no warranty       HYDR005
c     regarding accuracy or freedom from side-effects.                  HYDR006
c                   (C) 1986-1994 Bernhard Rupp                         HYDR007
c ----------------------------------------------------------------------HYDR008
c     HYDRAT estimates hydration shells for water atoms                 HYDR009
c ----------------------------------------------------------------------HYDR010
      
      character elemt*2                                                                HYDR011
      dimension did(iclust,idim),ihydtp(iat,2),ihyd(20,2),watat(iat,4)               HYDR012
      dimension elemt(iat)

c---  search all waters around water atom --
      ix=int(did(1,5))
      watat(ix,1)=ix
      iwats=0
      do 0001 i=2,itot
         if (elemt(did(i,5)).eq.'OH') then
            watat(ix,2)=did(i,4)
            iwats=iwats+1
         else
            watat(ix,3)=did(i,4)
            goto 0002
         end if
 0001 continue
 0002 continue      
      watat(ix,4)=iwats

c --- in case the closest is a nonwater, check for any waters within 4A
      if (watat(ix,4).lt.0.001)then
         do 0003 i=2,itot
            if(did(i,4).le.4.0) then
              if (elemt(did(i,5)).eq.'OH') then
                 watat(ix,2)=0.5555
              end if   
            else
              goto 0004
            end if
 0003    continue
 0004    continue      
      end if

c      write(*,*) 'watat(ix,i)=',(watat(ix,j),j=1,4)

c --- check for level of hydration ---                                  HYDR014
      iwat=iwat+1                                                       HYDR015
      ihydtp(iwat,1)=int(did(1,5))                                      HYDR016
      ihydtp(iwat,2)=1                                                  HYDR017
c --- first, isolated water atom ---                                    HYDR018
      if (ihyd(1,1).eq.0) then                                          HYDR019
         ihydtp(iwat,2)=0                                               HYDR020
      else                                                              HYDR021
         j=1                                                            HYDR022
 0102    continue                                                       HYDR023
         if ((ihyd(j,1).gt.0).and.(ihyd(j,2).eq.0)) then                HYDR024
            ihydtp(iwat,2)=ihydtp(iwat,2)+1                             HYDR025
            j=j+1                                                       HYDR026
            goto 0102                                                   HYDR027
         end if                                                         HYDR028
      end if                                                            HYDR029
      return                                                            HYDR030
      end                                                               HYDR031
