!C----------------------------------------------------------------------C
!C                                                                      C
!C    Parallel Metropolis algorithm for 2-dimensional Ising model.      C
!C                                                                      C
!C    HPF version 0.1   October 7, 1994                                 C
!C    Paul Coddington, paulc@npac.syr.edu                               C 
!C                                                                      C
!C    Main difference between parallel and sequential Metropolis        C
!C    program is that we need to do a checkerboard update, doing        C
!C   black sites first, then white sites. We're not allowed to         C
!C   update neighboring sites at the same time, since they are         C
!C   dependent on each other. Hence the black_white mask in the        C
!C   update subroutine.                                                C
!C                                                                     C
!C----------------------------------------------------------------------C

      PROGRAM metropolis
      USE global_variables
      INTEGER m,n, seed, start 
      REAL time1(2), time2(2) 

!C
!C  Initialize.
!C
!C --- Read in parameters.
      CALL input(seed,start)
 
!C --- Initialize the random number generator.
      CALL HPF_RANDOM_SEED(seed)
 
!C --- Initialize the spins.
      CALL init_spins(start)
 
!C --- Initialize the measurement data.
      CALL init_measure
 
!C --- Set up black/white checkerboard mask.
      CALL checkerboard
 
!C --- Initialize timer. 
!C      CALL dtime(time1)

!c      CALL picture
!c      CALL print_rand

!C
!C  Thermalize the system.
!C
      DO n = 1, therm_sweeps
        CALL HPF_RANDOM_NUMBER(ranum)
        CALL metupdate(black_white)
        black_white = .NOT. black_white
        CALL metupdate(black_white)
      END DO
!c      PRINT *
!c      PRINT *, 'Completed thermalization'

!C
!C  Metropolis update, taking measurements every 'sweeps' iterations.
!C
      DO m = 1, measurements
        DO n = 1, sweeps
          accept = 0
          CALL HPF_RANDOM_NUMBER(ranum)
          CALL metupdate(black_white)
          black_white = .NOT. black_white
          CALL metupdate(black_white)
        END DO
        acceptance = FLOAT( SUM(accept) ) / volume
        CALL measure
      END DO
 
!C --- Get time taken. 
!C      CALL dtime(time2)

!C
!C  Print out the results.
!C
      CALL output
!C     CALL picture

      PRINT *
!C      PRINT *, 'Time taken = ', time2(1) - time1(1)

      STOP
      END PROGRAM




!C----------------------------------------------------------------------C
!C                                                                     C
!C   The Metropolis Monte Carlo update on black or white sites.        C
!C                                                                     C
!C   For the Ising model new_spin = - old_spin, so                     C
!C   the old energy   oldE = - old_spin * (sum of neighboring spins)   C
!C   the new energy   newE = - new_spin * (sum of neighboring spins)   C
!C                         = - oldE                                    C
!C                                                                     C
!C   The Metropolis algorithm says that a lattice site should be       C
!C   updated if the change in energy deltaE = newE - oldE is less      C
!C   than or equal to zero, or with a probability                      C
!C      p = exp{-beta * deltaE}                                        C
!C   if the energy is increased. So an update mask is set up which     C
!C   specifies the sites for which the update is to be done, i.e.      C
!C      update = black .AND. ( (deltaE .LE. 0) .OR. (p .GT. ranum) )   C
!C   where ranum is an array of random numbers.                        C
!C   This update mask is set up in stages to avoid nested WHEREs and   C
!C   doing more calculations than necessary.                           C
!C                                                                     C
!C----------------------------------------------------------------------C
 
      SUBROUTINE metupdate(black)
      USE global_variables
      INTEGER i, j
      LOGICAL, DIMENSION(0:dimm1, 0:dimm1) :: black, update
      REAL, DIMENSION(0:dimm1, 0:dimm1) :: p
      INTEGER, DIMENSION(0:dimm1, 0:dimm1) :: new_spin
      INTEGER, DIMENSION(0:dimm1,0:dimm1) :: oldE, newE, deltaE

!HPF$ ALIGN update(:,:)   WITH LATTICE(:,:)
!HPF$ ALIGN new_spin(:,:) WITH LATTICE(:,:)
!HPF$ ALIGN oldE(:,:)     WITH LATTICE(:,:)
!HPF$ ALIGN newE(:,:)     WITH LATTICE(:,:)
!HPF$ ALIGN deltaE(:,:)   WITH LATTICE(:,:)
!HPF$ ALIGN p(:,:)        WITH LATTICE(:,:)

!C
!C  New trial spin (= - old spin)
!C
      new_spin = -spin
!C
!C  Calculate old and new energy (new energy = - old energy)
!C  and initialize the update mask.
!C
      update = .FALSE.
      WHERE (black)
        oldE = - spin *  &
     &         ( CSHIFT(spin,1,1) + CSHIFT(spin,-1,1) + &
     &           CSHIFT(spin,1,2) + CSHIFT(spin,-1,2) )
        newE = -oldE
        deltaE = newE - oldE
        update = deltaE .LE. 0 
      END WHERE
!C
!C  Calculate the update mask (don't want to calculate the 
!C  exponential unless we have to, so do it separately here).
!C 
      WHERE ( black .AND. (deltaE .GT. 0) )
        p = EXP(-beta*deltaE)
	update = p .GT. ranum 
      END WHERE
!C
!C  Metropolis accept/reject 
!C 
      WHERE ( update )
          spin = new_spin
	  accept = 1
      END WHERE

      RETURN
      END SUBROUTINE


