
!-----------------------------------------------------------------------
!
!  Q State 2d Potts Model in Fortran 90 / HPF
!  K.A.Hawick, NPAC, August 1994.
!
!-----------------------------------------------------------------------
      program potts2d
      integer, parameter :: m = 64   ! Lattice dimensions
      integer, parameter :: n = 64

!HPF$ PROCESSORS PROCS(8)
!HPF$ TEMPLATE, DISTRIBUTE(*,BLOCK) ONTO PROCS :: RECTANGLE(m,n)

      integer, parameter :: q = 4     ! Number of allowable States for Degrees of Freedom
      real, parameter :: cup = -4.0   ! Coupling in kT units, -ve for ferromagnetic

      integer, dimension( 1:m, 1:n ) :: latt  ! The lattice of DOF's
!HPF$ ALIGN latt(:,:) WITH RECTANGLE(:,:)

      integer, dimension( 1:m, 1:n ) :: itmp  ! Trial choices for DOF's
!HPF$ ALIGN itmp(:,:) WITH RECTANGLE(:,:)

      real, dimension( 1:m, 1:n ) :: tmp, p   ! TMP for RNG and Probability 
!HPF$ ALIGN tmp(:,:) WITH RECTANGLE(:,:)
!HPF$ ALIGN p(:,:)   WITH RECTANGLE(:,:)

      logical, dimension( 1:m, 1:n ) :: mask  ! mask used in ensuring detailed balance prevails
!HPF$ ALIGN mask(:,:) WITH RECTANGLE(:,:)

      integer :: nmeasure = 20                ! number of measurements
      integer :: niter = 10                   ! number of sweeps per measurement

      call setup_random( 12345 )

      call random_number( tmp )    ! initialise between 0 and q-1 inclusive
      latt = int( tmp * q )

      do imeasure=1,nmeasure

! Measurements:

        itmp =        merge( 1, 0, latt == cshift(latt,  1, 1 ) )   ! Present bonds
        itmp = itmp + merge( 1, 0, latt == cshift(latt, -1, 1 ) )
        itmp = itmp + merge( 1, 0, latt == cshift(latt,  1, 2 ) )
        itmp = itmp + merge( 1, 0, latt == cshift(latt, -1, 2 ) )
        write(6,*) '# Bonds = ', sum( itmp )
!        do iq=1,q
!          write(6,*) '#Q= ',iq,'   is ', count( latt == iq )
!        enddo

      t1 = etime(x)
      do iter=1,niter

        call random_number(tmp)    ! generate "new" values for the Potts variables
        itmp = mod( (latt + int( tmp * q-1 ) + 1 ), q )  ! using the Coddington recipe


        do i0=1,2    ! use four-colour ordering to preserve detailed-balance
        do j0=1,2

          call random_number( tmp )  ! have to get too many cos of this stupid non-function RNG

          mask = .false.
          forall(i=i0:m:2, j=j0:n:2 ) mask(i,j) = .true.
          where( mask )
            p =     merge( cup, 0.0, latt == cshift(latt,  1, 1 ) )   ! Present bonds
            p = p + merge( cup, 0.0, latt == cshift(latt, -1, 1 ) )
            p = p + merge( cup, 0.0, latt == cshift(latt,  1, 2 ) )
            p = p + merge( cup, 0.0, latt == cshift(latt, -1, 2 ) )

            p = p - merge( cup, 0.0, latt == cshift(itmp,  1, 1 ) )   ! future bonds
            p = p - merge( cup, 0.0, latt == cshift(itmp, -1, 1 ) )   ! if dof is flipped
            p = p - merge( cup, 0.0, latt == cshift(itmp,  1, 2 ) )
            p = p - merge( cup, 0.0, latt == cshift(itmp, -1, 2 ) )

            p = exp( p )
            latt = merge( itmp, latt, tmp < p )
          end where

        end do  !
        end do  ! end of checkerboard loop

      end do  ! end of iter loop
      t2 = etime(x)
      write(6,*) 'Time: ', t2 - t1
      end do  ! end of measure loop


      write(9,100)  ! write out picture of configuration as a PPM file
  100 Format('P6' )
      write(9,101)
  101 Format('# 2d Potts Model')
      write(9,*) m, n
      write(9,*) 127
      do i=1,m
      do j=1,m
        if( latt(i,j) == 0 )then
          write(9,'(3A1)',ADVANCE='NO')achar(127), achar(0),achar(0)
        else if( latt(i,j) == 1 )then
          write(9,'(3A1)',ADVANCE='NO')achar(0), achar(127),achar(0)
        else if( latt(i,j) == 2 )then
          write(9,'(3A1)',ADVANCE='NO')achar(0), achar(0),achar(127)
        else if( latt(i,j) == 3 )then
          write(9,'(3A1)',ADVANCE='NO')achar(127), achar(127),achar(0)
        else if( latt(i,j) == 4 )then
          write(9,'(3A1)',ADVANCE='NO')achar(127), achar(0),achar(127)
        else if( latt(i,j) == 5 )then
          write(9,'(3A1)',ADVANCE='NO')achar(0), achar(127),achar(127)
        else if( latt(i,j) == 6 )then
          write(9,'(3A1)',ADVANCE='NO')achar(127), achar(127),achar(127)
        else
          write(9,'(3A1)',ADVANCE='NO')achar(0), achar(0),achar(0)
        endif
      enddo
      enddo




        
      end program 

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

      subroutine setup_random(iseed)  ! put a new seed in the bloody awful system supplied RNG
      integer :: n
      integer, dimension(:), allocatable ::  seed
      call random_seed(size=n)
      allocate( seed(n) )
      forall(i=1:n) seed(i) = mod(i*4096+150889+iseed,714025 )
      call random_seed( put=seed(1:n) )
      deallocate( seed )
      end subroutine


