!**********************************************************************
!          
!     Parallel  Wolff cluster Monte Carlo algorithm for the   
!     2-dimensional Ising model. (HPF version)
!               
!     Uses the parallel Ants-in-the-labyrinth (Breadth-first-search)
!     algorithm.   
!
!     This program is modified from the program written by
!     John Apostolakis and Paul Coddington. ( Some routines for 
!     Swendsen-Wang algorithm are replaced with new ones.)
!                
!**********************************************************************


      MODULE Var
      IMPLICIT NONE
      INTEGER dim, dimm1, Nproc, MAXINT, Maxiter
      REAL*4  fdim, vol
      PARAMETER (dim     = 8)
      PARAMETER (fdim    = 8.0)
      PARAMETER (Nproc   = 8)
      PARAMETER (dimm1   = dim - 1)
      PARAMETER (vol     = fdim*fdim)
      PARAMETER (MAXINT  = 2147483647)
      PARAMETER (Maxiter = 10000000)

      integer, dimension(0:dimm1, 0:dimm1) :: spin
      real,    dimension(0:dimm1, 0:dimm1) :: ranum
      logical, dimension(0:dimm1, 0:dimm1) :: cluster, new_cluster
      logical, array(0:3, 0:dimm1, 0:dimm1) :: neighbor


      INTEGER therm_sweeps, sweeps, measurements
      INTEGER nits, size
      REAL*4  beta, p
      REAL*8  sum_energy, sum_enersq, sum_spins, sum_spinsq
      REAL*8  sum_nits, sum_nitsq, sum_size, sum_sizesq
      END MODULE Var

!
!

      MODULE Directions
      integer, parameter :: right  = 0
      integer, parameter :: left   = 1
      integer, parameter :: down   = 2
      integer, parameter :: up     = 3
      END MODULE Directions


!**********************************************************************
!
!      Main Program
!
!**********************************************************************

 


      PROGRAM Wolff_Cluster
      USE Var
      INTEGER m, dummy
      REAL    time, t1, t2, total_sweeps
      REAL*8  elapsed_time, persite_time, scale_time 


!
!  Read in parameters, initialize spins, etc.
!
      CALL initialize
 
      CALL random
      sum_size = 100

!
!  Thermalize the configuration.
!
!      CALL wolff_update(therm_sweeps)

!
!  Wolff update of the configuration, taking measurements
!  every 'sweeps' iterations.
!
      time = 0.

      DO m = 1, measurements
        t1 = dtime(dummy)
!        CALL wolff_update(sweeps)
        t2 = dtime(dummy)
        time = time + t2
!        CALL measure
      END DO
 
!
!  Print out the results.
!
      CALL output

        total_sweeps = (sweeps*measurements)
        
        WRITE(6,10) '================================================'
        WRITE(6,*) '     Time for', sweeps*measurements,'sweeps'
        WRITE(6,*) '------------------------------------------------'
        WRITE(6,20) ' Time          =', time, 'sec'
        WRITE(6,10) '================================================'

        elapsed_time = (time / total_sweeps)*1000
        persite_time = elapsed_time / sum_size
        scale_time   = (elapsed_time/sum_size) * (dim*dim)
        
        WRITE(6,10) '*************************************************'
        WRITE(6,*) '           Timing Data for One Sweep'
        WRITE(6,*) '*************************************************'
        WRITE(6,20) ' Sweep Time    =', elapsed_time, 'msec'
        WRITE(6,20) ' Time per site =', persite_time, 'msec'
        WRITE(6,20) ' Scale Time    =', scale_time,   'msec'
        WRITE(6,*) '*************************************************' 
        WRITE(6,*) ' '
        WRITE(6,*) ' '
 10     FORMAT(//,1X,A)
 20     FORMAT(4X,A,F12.3,2X,A) 
 

!      call print_spin
 
      END program Wolff_Cluster



!********************************************************************
!
!  Initialize parameters and spins.
!
!********************************************************************

      SUBROUTINE initialize
      USE Var
      INTEGER i, j, start, seed
      LOGICAL ordered
      integer mode

!--- Distribution directive
!HPF$ PROCESSORS PROCS(Nproc)
!!HPF$ PROCESSORS PROCS(NUMBER_OF_PROCESSORS())
!HPF$ TEMPLATE, DISTRIBUTE(*,block) ONTO PROCS ::   &
                RECTANGLE(0:dimm1, 0:dimm1)
!HPF$ ALIGN spin WITH RECTANGLE
!HPF$ ALIGN ranum WITH RECTANGLE
!HPF$ ALIGN cluster WITH RECTANGLE
!HPF$ ALIGN new_cluster WITH RECTANGLE

!   !HPF$ ALIGN neighbor WITH RECTANGLE
!HPF$ ALIGN neighbor(*,J,K) WITH RECTANGLE(J,K)
!!HPF$ DISTRIBUTE neighbor(*,*,BLOCK)



!
!  Read in parameters.
!
      print *, ' ' 
      PRINT *, ' input the number of thermalization sweeps'
      READ *,  therm_sweeps
      PRINT *, ' input the number of sweeps per measurement'
      READ *,  sweeps
      PRINT *, ' input the number of measurements'
      READ *,  measurements
      PRINT *, ' input the beta value'
      READ *,  beta
      PRINT *, ' ordered start (0) or random start (1)?'
      READ *,  start
      PRINT *, ' input seed for random number generator'
      READ *,  seed

      p = 1 - EXP(-2.0*beta)

      ordered = start.EQ.0

!
!  Initialize random number generator.

      WRITE(6,*) ' Random seed; before            = ', seed
!      CALL random_seed(seed)
!      CALL setup_random(seed)
      CALL set_init_seed(seed, Nproc)
      WRITE(6,*) ' Random seed; after             = ', seed
!
!  Initialize spins.
!
      IF (ordered) THEN
        spin = 1
      ELSE
        CALL random
        WHERE (ranum .GE. 0.5)
          spin = 1
         ELSEWHERE
          spin = -1
        END WHERE
      END IF
!
!  Initialize sums.
!
      sum_energy = 0.0
      sum_spins  = 0.0
      sum_enersq = 0.0
      sum_spinsq = 0.0
      sum_nits   = 0.0
     
      WRITE(6,*) ' '
      WRITE(6,*) ' '
      WRITE(6,*) '****************************************************'
      WRITE(6,*) '         Wolff Cluster (HPF version) '
      WRITE(6,*) '----------------------------------------------------'
      WRITE(6,*) '                  input  data     '
      WRITE(6,*) '****************************************************'
      WRITE(6,22) ' Domain size               = ', dim, 'x', dim,  &
                                            '=',  dim*dim
!      WRITE(6,*) ' Number of processors      = ', nproc
      WRITE(6,*) ' Thermalization sweeps     = ', therm_sweeps
      WRITE(6,*) ' Sweeps per measurements   = ', sweeps
      WRITE(6,*) ' Number of measurments     = ', measurements
      WRITE(6,*) ' Beta value                = ', beta
      WRITE(6,*) ' Wolff probability         = ', p
      WRITE(6,*) ' Starting state            = ', start
      WRITE(6,*) ' Random seed               = ', seed
      WRITE(6,*) '****************************************************'
      WRITE(6,*) ' '
      WRITE(6,*) ' '
 22   FORMAT (1X,A,1X,I4,1X,A,1X,I4,1X,A,1X,I7)

      RETURN
      END



!********************************************************************
!
!
      SUBROUTINE measure
      USE Var
      INTEGER i, j
      REAL*4 ene, mag
      INTEGER e(0:dimm1, 0:dimm1)
!HPF$ ALIGN e WITH spin

!
!  Energy
!


!      print *, 'in measure'

      e = -spin * ( CSHIFT(spin,1,1) + CSHIFT(spin,2,1 ) )
      ene = sum(e)/vol
      sum_energy = sum_energy + ene
      sum_enersq = sum_enersq + (ene * ene)
!
!  Magnetization

      mag = sum(spin)/vol
      IF (mag.LT.0.0) mag = -mag
      sum_spins = sum_spins + mag
      sum_spinsq = sum_spinsq + (mag * mag)
!
!  Number of labeling iterations
!
      sum_nits = sum_nits + nits
      sum_nitsq = sum_nitsq + (nits * nits)

!
!  Cluster Size
!
      sum_size = sum_size + size
      sum_sizesq = sum_sizesq + ( size*size )

!     PRINT *, 'Ene ', ene,'    Mag ', mag

      RETURN
      END





!
!  Take averages, find errors and variances, and print out the results.
!

      SUBROUTINE output
      USE Var
      REAL*8 err_ener, err_spin, specific_heat, susceptibility
      REAL*8 err_nits, variance, err_size, var_size
!
!  Take averages.
!

!      print *, 'in output'

      sum_energy = sum_energy / measurements
      sum_enersq = sum_enersq / measurements
      sum_spins  = sum_spins  / measurements
      sum_spinsq = sum_spinsq / measurements
      sum_nits   = sum_nits  / measurements
      sum_nitsq  = sum_nitsq / measurements
      sum_size   = sum_size / measurements
      sum_sizesq = sum_sizesq / measurements
!
!  Work out errors and variances.
!
      specific_heat  = sum_enersq - (sum_energy * sum_energy)
      susceptibility = sum_spinsq - (sum_spins  * sum_spins)
      variance       = sum_nitsq  - (sum_nits   * sum_nits)
      var_size       = sum_sizesq  - (sum_size   * sum_size)
      err_ener = SQRT( abs(specific_heat  / measurements ))
      err_spin = SQRT( abs(susceptibility / measurements ))
      err_nits = SQRT( abs(variance       / measurements ))
      err_size = SQRT( abs(var_size       / measurements ))
      specific_heat  = specific_heat  * vol
      susceptibility = susceptibility * vol
!
!  Print out results.
!
      PRINT *, ' '
      PRINT 10, '  Average energy is ', sum_energy,' +/-',err_ener
      PRINT 10, '  Average magnet is ', sum_spins, ' +/-',err_spin
      PRINT *, ' '
      PRINT 20, '  Specific heat  is ', specific_heat
      PRINT 20, '  Susceptibility is ', susceptibility
      PRINT *, ' '
      PRINT 10, '  Average iters is ', sum_nits,' +/-',err_nits
      print 30, '  Average cluster size is ', sum_size, ' +/-',err_size
  10  FORMAT(A,F20.6,A,F12.6)
  20  FORMAT(A,F12.5)
 30   format(a,f20.6,a,f15.6//)
      RETURN
      END



!*********************************************************************
!
!     Wolff cluster update.
!
!*********************************************************************

      subroutine Wolff_update(iterations)
      USE Var
      integer iter, iterations

!      print *, 'in wolff_update'



      do iter = 1, iterations

!---  Set up neighbors for each site
         call setup_neighbors()

!---  Make a single cluster  
      
         call make_cluster()

!--- Check the cluster

         call check_cluster()

!---  Update the spins in the cluster  


         size = count(cluster)

         where (cluster)
           spin = -1 * spin
         end where

!D        print *, ' Sweep = ', iter

      end do
      

      end subroutine Wolff_update



!**********************************************************************
! 
!    Set up the four neighbors for each site.  
!
!     If spin[i, j] is connected with spin[i, j+1]           
!        then neighbor[right, i, j+1] = true and           
!             neighbor[left , i, j  ] = true.
!
!      *---------*
!     [i,j]     [i,j+1]
!
!     If spin[i, j] is connected with s[i+1, j]
!        then neighbor[down, i+1, j] = true and
!             neighbor[up  , i, j  ] = true.
!
!         *  [i,j]
!         |
!         |
!         *  [i+1,j]
!
!**********************************************************************

      subroutine setup_neighbors()
      USE Var
      USE Directions

!      print *, 'in setup'

      call random

!      call print_rand
!      print *, '** 1 ** '
!      print *, ' '

      neighbor(right, :, :) = ( cshift(spin, -1, 2) == spin )  &
                        .and. ( ranum < p )

      neighbor(left, :, :)  = cshift( neighbor(right, :, :), 1, 2 )

      call random

!      call print_rand
!      print *, '** 2 ** '   
!      print *, ' '



      neighbor(down, :, :)  = ( cshift(spin, -1, 1) == spin )  &
                         .and. ( ranum < p )

      neighbor(up, :, :)    = cshift( neighbor(down, :, :), 1, 1 )
       
  
      end subroutine setup_neighbors

       

!***********************************************************************
!
!   Using the parallel Ants-in-the-labyrinth (Breadth-first-search)
!   algorithm, make a single cluster
!
!***********************************************************************

      subroutine make_cluster()
      USE Var
      USE Directions

      logical exit
      integer counter, x, y
      integer, array(0:1) :: rand_site

!--- Initialize variables

      exit        = .false.
      counter     = 0
      cluster     = .false.
      new_cluster = .false.

!--- Choose a random site, so that it becomes the first element 
!--- of the cluster. It also belongs to the new cluster.

!?????????????????????????????????????????????
!      call random_number(rand_site, dim)
      call random

!????????????????????????????????????????????????


      rand_site(0) = ranum(0,0)*dim
      rand_site(1) = ranum(0,1)*dim
      new_cluster( rand_site(0), rand_site(1) ) = .true.
!      new_cluster( x, y ) = .true.

!      print *, 'rand site:', rand_site(0),rand_site(1)

!--- Repeat until the cluster does not change.

      do while ( .not. exit )

         counter = counter + 1

         cluster = cluster .or. new_cluster   ! update the cluster

!       Find neighbors of each spin which belongs to newcluster.
!       Those neighbors will belong to the new cluster        

         new_cluster =                                               &
          (                                                          &
           (neighbor(right, :, :) .and. cshift(new_cluster, -1, 2))  &
          .or.                                                       &
          (neighbor(left  , :, :) .and. cshift(new_cluster, 1, 2))   &
          .or.                                                       &
          (neighbor(down  , :, :) .and. cshift(new_cluster, -1, 1))  &
          .or.                                                       &
          (neighbor(up    , :, :) .and. cshift(new_cluster, 1,  1))  &
          )                                                          &
        .and.                                                        &
          ( .not. cluster )

!       If new cluster is empty, or counter > Maxcount, then stop

         if ( all(.not. new_cluster) .or. (counter > Maxiter))  & 
           exit = .true.

      end do

      nits = counter


      end subroutine make_cluster



!***********************************************************************
!
!   Check that the cluster is consistent.
!
!***********************************************************************

      subroutine check_cluster()
      USE Var
      USE Directions
      logical, array(0:dimm1, 0:dimm1) :: check, not_neighbor
      integer numbad
      integer i, j
      logical ok
!HPF$ ALIGN check WITH spin
!HPF$ ALIGN not_neighbor WITH spin

!      print *, 'in check'

      not_neighbor = .not. neighbor(left, :, :)
      check = (cluster == cshift(cluster, 1, 2)) .or.  &
             not_neighbor
      not_neighbor = .not. neighbor(up, :, :)
      check = ( (cluster == cshift(cluster, 1, 1)) .or. &
               not_neighbor )   .and.   check

      ok = all(check)

      if ( .not. ok ) then
         if ( nits > Maxiter ) then
            print 10, Maxiter + 1
 10         format(1x,'//Too many iteration to continue! ',  &
                     'Current iteration = ', I5)
         else
            numbad = sum(merge(0, 1, check))
            print *, ' Check of cluster fails!! Sum = ', numbad
         end if
      end if     
      
      end subroutine check_cluster

!-----------------------------------------------------------------

      subroutine print_spin
      USE Var
      integer i,j

      do i=0, dimm1
      WRITE(6,33) (spin(i,j), j=0,dimm1)
      enddo
 33   FORMAT(32 I3)

      end subroutine print_spin

!----------------------------------------------------------------
      subroutine print_rand
      USE Var
      integer i,j

      do i=0, dimm1
 !     WRITE(6,33) (ranum(i,j), j=0,dimm1)
         print *, (ranum(i,j), j=0,dimm1)
      enddo
 33   FORMAT(32 F5.2)
      
      print *, ' '
      print *, ' '

      end subroutine print_rand




!**************************************************************
!--------------------------------------------------------------
  
