
C----------------------------------------------------------------------C
C                                                                      C
C       Metropolis Monte Carlo update of the Ising model.              C
C      (Metropolis et al., J. Chem. Phys. 21 (1953) 1087).             C
C                                                                      C
C    Based on the code in the book "Statistical Field Theory",         C
C    G. Parisi, (Addison-Wesley, 1988).                                C
C                                                                      C
C       Paul Coddington, August 1992, June 1993.                       C
C                                                                      C
C----------------------------------------------------------------------C

	PROGRAM metropolis
	INCLUDE "common.h"
	INTEGER i,seed,start

C
C  Initialize.
C
C ----  Initialize some Express stuff.
	CALL KXINIT

C ----  Read in the parameters.
	CALL input(seed,start)

C ----  Do domain decomposition of the lattice over the nodes.
	CALL decompose
 
C ----  Initialize the random number generator.
	CALL randinit(seed)
 
C ----  Initialize the spins.
	CALL init_spins(start)
 
C ----  Initialize the measurement data.
	CALL init_measure

        WRITE(6,*) 'finished initialization'
	CALL KFLUSH(6)
C
C  Thermalize the system.
C
	CALL metupdate(therm_sweeps)

        WRITE(6,*) 'finished thermalization'
	CALL KFLUSH(6)

C
C  Metropolis update, taking measurements every 'sweeps' iterations.
C
	DO i=1,measurements
	  CALL metupdate(sweeps)
	  CALL measure
	ENDDO
 
C
C  Print out the results.
C
	CALL output

	STOP
	END



C----------------------------------------------------------------------C
C                                                                      C
C                                                                      C
C    Metropolis Monte Carlo update of the Ising model.                 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    Do black/white checkerboard update, so can pass edge values       C
C    (blocked data) which will be more efficient than passing          C
C    single spins at a time (unblocked data).                          C
C                                                                      C
C----------------------------------------------------------------------C

	SUBROUTINE metupdate(iterations)
	INCLUDE "common.h"
	INTEGER iterations, black_white
	INTEGER n, i,j, jstart
        INTEGER new_spin, newE,oldE,deltaE
        REAL rand
	
	DO n=1,iterations

	  accept = 0

C ------  Checkerboard update.
C ------  Stride of inner loop is 2, starting point of inner loop
C ------  depends on whether we are updating black or white sites.

	  DO black_white = 1,2

C --------  Pass edge values to get neighboring spins.
	    CALL all_edges
             
	    jstart = black_white
	    jstart = 3 - jstart
	    DO i=1,isiz
C --------- Alternate between jstart=1 and jstart=2
	    jstart = 3 - jstart
	    DO j=jstart,jsiz,2
C
C  New trial spin (= - old spin)
C
              new_spin = - spin(i,j)

C
C  Calculate old and new energy (new energy = - old energy).
C  (May need to get some neighboring spin values values from 
C   other processors).
C
	      oldE = - spin(i,j) * 
     &                ( spin(i-1,j) + spin(i+1,j) +
     &                  spin(i,j-1) + spin(i,j+1) )
              newE = - oldE
              deltaE = newE - oldE

C
C  Metropolis accept/reject 
C 
              IF ( ( deltaE.LE.0 ) .OR.
     &             ( EXP(-beta*deltaE).GT.rand() ) ) THEN
                spin(i,j) = new_spin
                accept = accept + 1
              ENDIF

	    ENDDO
	    ENDDO

C ------  End of loop over black/white
	  ENDDO

	ENDDO

C
C       Pass edges to neighboring processors.
C
        CALL spin_edges

	RETURN 
	END


