
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    As long as isiz and jsiz are both >=2, we don't have to do a      C
C    black/white checkerboard update, since sites being updated        C
C    simultaneously on different processors will not be neighbors.     C
C    However a checkerboard update is necessary if we want to do       C
C    efficient blocked data communications.                            C
C                                                                      C
C----------------------------------------------------------------------C

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

          accept = 0

          DO i=1,isiz
          DO j=1,jsiz

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) * 
     &              ( get_spin(i-1,j) + get_spin(i+1,j) +
     &                get_spin(i,j-1) + get_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

        ENDDO

C
C       Pass edges to neighboring processors.
C
        CALL spin_edges

	RETURN 
	END


