*****************************************************************************
C                            ADSORB_SQR3.F
C Returns the location of the graphite adsorbtion sites.
C Also determines the sublattice (A=1, B=2, C=3) that each particle 
C belongs to.  This is stored in the array isub.
*****************************************************************************
      subroutine adsorb(x_ads,y_ads,isub)
      include 'pimc_sqr3sc.dcl'
      dimension x_ads(2*ntile_x*3*ntile_y)
      dimension y_ads(2*ntile_x*3*ntile_y)
      dimension rx1(2),ry1(2)
      integer isub(2*ntile_x*3*ntile_y)
      integer icyclic(3)

      a1=a_sqr3sc    ! Use 4.26
      a2=a_gr      ! Use 2.46

*---------------------------------------------------------------------------

      rx1(1)=0.0     
      ry1(1)=0.0
      rx1(2)=a1/2.
      ry1(2)=a2/2.

*---------------------------------------------------------------------------
C Place the frozen helium atoms at hcp lattice sites.
*---------------------------------------------------------------------------
      icyclic(1)=3
      icyclic(2)=1
      icyclic(3)=2

      iads=0
      isublat=0
      do mm=1,ntile_x  
         isublat=1
         cell_x=real(mm-1)*a1   
         do nn=1,3*ntile_y
            cell_y=real(nn-1)*a2
            do igr=1,2
               iads=iads+1
               x_ads(iads)=cell_x+rx1(igr)
               y_ads(iads)=cell_y+ry1(igr)
               if(igr.eq.1) then
                  isub(iads)=isublat
               else
                  isub(iads)=icyclic(isublat)
               end if
c               print*,x_ads(iads),y_ads(iads)
            end do
            isublat=isublat+1
            if(isublat.gt.3) isublat=1
         end do
      end do

c      do mm=1,ntile_x
c         cell_x=real(mm-1)*a1
c         x=cell_x+rx1(1)
c         y=ntile_y*a2+ry1(1)
c         print*,x,y
c      end do
c      do mm=1,ntile_y
c         cell_y=real(mm-1)*a2
c         y=cell_y+ry1(1)
c         x=ntile_grx*a1+rx1(1)
c         print*,x,y
c      end do
c      x=ntile_x*a1+rx1(1)
c      y=3*ntile_y*a2+ry1(1)
c      print*,x,y
      end
