******************************************************************************
C                            PIMCSSF_GRSOL.F
C Calculate the energy and superfluid density (using winding numbers).
******************************************************************************
C Subroutines and functions:
      include 'metrop_sqr3sc.f'
      include 'getdatadu_sqr3sc.f'
c      include 'frozen_he_layer.f'
      include 'frozen_he_II.f'
      include 'initialize_sqr3sc.f'
      include 'permid_sqr3sc.f'
      include 'mgauss_sqr3sc.f'
      include 'samdist_sqr3sc.f'
      include 'ssfobsv_sqr3sc.f'
      include 'dstpbc_sqr3sc.f'
      include 'ratint.f'
      include 'gasdev.f'
      include 'ran1.f'
      include 'ran2.f'
      include 'pimc_sqr3sc.dcl'
******************************************************************************

      dimension x_sol(nfrozen),y_sol(nfrozen),z_sol(nfrozen)
      dimension x(2*m,nat),y(2*m,nat),z(2*m,nat)
      dimension xp(2*m,nat),yp(2*m,nat),zp(2*m,nat)
      dimension r1(nr1),r2(nr1),u00(nr1)
      dimension dutmov(lev,nr1),ddutmov(lev,nr1)
      dimension dutfix(lev,nr1),ddutfix(lev,nr1)
      dimension u(npol,nr1)
      dimension rdudt1(nr1),rdudt2(nr1),dudt1(nr1),dudt2(npol,nr1)
      dimension rdudl1(nr1),rdudl2(nr1),dudl1(nr1),dudl2(npol,nr1)
      dimension rgrlat(nr1),ugrlat(nr1),ugrdtlat(nr1),ugrdllat(nr1)
      dimension ssfx(nssfx),ssfy(nssfy)

      integer im(npol),in(npol)
      integer image(2*m,nat)

      data im/0,1,1,0,2/
      data in/1,0,1,2,0/

      call getdata(r1,r2,rdudt1,rdudt2,rdudl1,rdudl2,rgrlat,dutmov,
     x     ddutmov,dutfix,ddutfix,u00,u,dudt1,dudt2,dudl1,dudl2,ugrlat,
     x     ugrdtlat,ugrdllat,nd,nu,n1,n2,n3,n4,ngrlat)

      open(99,file='data/pthermsqr3sct40m80d208n32.dat',status='old')
c      open(100,file='data/pthermsqr3sct40m80d208n32.dat',
c     x     status='unknown')
      open(101,file='ssfxt40m80d208n32.dat',status='unknown')
      open(102,file='ssfyt40m80d208n32.dat',status='unknown')

      itherm=1             ! 0 to start anew, 1 to start from old data

*//////////////////////////////////////////////////////////////////////////////
C Read positions of the frozen helium sites.
*//////////////////////////////////////////////////////////////////////////////

      call slab(x_sol,y_sol,z_sol)

*//////////////////////////////////////////////////////////////////////////////
C Initialize particle positions.  For ITHERM=0, use subroutine SOLID.F.
C For ITHERM=1, read positions from old data data file.
*//////////////////////////////////////////////////////////////////////////////
      
      if(itherm.eq.0) then
         call initialize(x,y,z,xp,yp,zp)
         do i=1,2*m
            do j=1,nat
               image(i,j)=j
            end do
         end do
         idum=-abs(int(x(1,1)*1.e5))
      else
         do i=1,2*m
            do j=1,nat
               read(99,*) x(i,j),y(i,j),z(i,j),image(i,j)
               xp(i,j)=x(i,j)
               yp(i,j)=y(i,j)
               zp(i,j)=z(i,j)
            end do
         end do
         idum=-abs(int(x(1,1)*1.e5))
      end if
      
*//////////////////////////////////////////////////////////////////////////////
      
      icount2=0
      icount3=0
      icount4=0
      nacc=0
      sum=0.
      fac=1./(4.*float(m)*tau*float(nat+nfrozen)*alambda)
      do iter=1,niter
         do ibin=1,nbin
            call metrop(x,y,z,xp,yp,zp,x_sol,y_sol,z_sol,r1,r2,rgrlat,
     x           dutmov,ddutmov,dutfix,ddutfix,u00,u,ugrlat,image,im,
     x           in,idum,nd,nu,ngrlat,icount2,icount3,icount4,nacc)
         end do
         call ssfobsv(x,y,z,ssfx,ssfy)
         print*,ssfx()
      end do

      pi=4.*atan(1.)
      qx=2.*pi/xdim
      qy=2.*pi/ydim
      qxy=sqrt(qx*qx+qy*qy)

      fac=float(niter)*float(nbin)*float(m)*float(nat)
      do i=1,nssfx
         write(101,*) qx*i,1.+ssfx(i)/fac
      end do

      do i=1,nssfy
         write(102,*) qy*i,1.+ssfy(i)/fac
      end do

c      do i=1,2*m
c         do j=1,nat
c            write(100,*) x(i,j),y(i,j),z(i,j),image(i,j)
c         end do
c      end do
      end
