PROGRAM ising
*     Monte Carlo simulation of the Ising model on the square lattice
*     using the Metropolis algorithm
      IMPLICIT NONE
      REAL*8 Ce(0:20), Cm(0:20), E, M, T, accept, accum(10),
     +       esave(100), msave(100), w(-8:8)
      INTEGER IR_, L, N, i, mcs, nsave, nequil, pass, spin(32,32)
      CALL GWopen(IR_, 0)
      CALL initial(N,L,T,spin,mcs,nequil,w,E,M)
      CALL plotspins(L,spin)
      DO i = 1, nequil
         CALL Metropolis(N,L,spin,E,M,w,accept)
      END DO
      CALL initialize(accum,accept)
      CALL initialize_correl(Ce,Cm,esave,msave,nsave)
      DO pass = 1, mcs
         CALL Metropolis(N,L,spin,E,M,w,accept)
         CALL mydata(E,M,accum)
         CALL correl(Ce,Cm,E,M,esave,msave,pass,nsave)
      END DO
      CALL output(T,N,mcs,accum,accept)
      CALL c_output(Ce,Cm,accum,mcs,nsave)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(N,L,T,spin,mcs,nequil,w,E,M)
      IMPLICIT NONE
      REAL*8 sum, T, w(-8:8), E, M
      REAL dmy, rnd
      INTEGER N, L, mcs, nequil, dE, right, up, x, y, spin(32,32)
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'linear dimension of lattice = '
      READ(*,*) L
      N = L*L                               ! number of spins
*     temperature measured in units of J/k
      WRITE(*,'(A,$)') 'temperature = '
      READ(*,*) T
      WRITE(*,'(A,$)') '# MC steps per spin for equilibrium = '
      READ(*,*) nequil
      WRITE(*,'(A,$)') '# MC steps per spin for data = '
      READ(*,*) mcs
      M = 0
      DO y = 1, L
         DO x = 1, L
            IF(rnd(0) .LT. 0.5D0) THEN
               spin(x,y) = 1                ! spin up
            ELSE
               spin(x,y) = -1
            END IF
            M = M + spin(x,y)               ! total magnetization
         END DO
      END DO
      E = 0
      DO y = 1, L
         IF(y .EQ. L) THEN
            up = 1                          ! periodic boundary conditions
         ELSE
            up = y + 1
         END IF
         DO x = 1, L
            IF(x .EQ. L) THEN
               right = 1
            ELSE
               right = x + 1
            END IF
            sum = spin(x,up) + spin(right,y)
            E = E - spin(x,y)*sum           ! total energy
         END DO
      END DO
*     compute Boltzmann probability ratios
      DO dE = -8, 8, 4
         w(dE) = exp(-dE/T)
      END DO
      END
C
      SUBROUTINE plotspins(L,spin)
      IMPLICIT NONE
      REAL xwin, ywin
      INTEGER IR_, i, j, L, spin(32,32)
      CALL GWclear(IR_, -1)
      CALL compute_aspect_ratio(REAL(L+1),xwin,ywin)
      CALL GWindow(IR_, 0.0, 0.0, xwin, ywin)
      DO i = 1, L
         DO j = 1, L
            CALL GWerase(IR_, (i-1)*L+j, 0);
            CALL GWsetogn(IR_, 0, (i-1)*L+j);
            IF(spin(i,j) .EQ. 1) THEN
               CALL GWsrect(IR_, REAL(i), REAL(j), i+1.0, j+1.0, 13)
            ELSE
               CALL GWsrect(IR_, REAL(i), REAL(j), i+1.0, j+1.0, 14)
            END IF
         END DO
      END DO
      END
C
      SUBROUTINE initialize_correl(Ce,Cm,esave,msave,nsave)
      IMPLICIT NONE
      REAL*8 Ce(0:20), Cm(0:20), esave(100), msave(100)
      INTEGER i, nsave
      nsave = 10
      DO i = 1, nsave
         esave(i) = 0
         msave(i) = 0
      END DO
      DO i = 1, nsave
         Ce(i) = 0
         Cm(i) = 0
      END DO
      END
C
      SUBROUTINE initialize(accum,accept)
*     use array to save accumulated values of magnetization and
*     energy. Array used to make it easier to add other quantities      
      IMPLICIT NONE
      REAL*8 accum(10), accept
      INTEGER i
      DO i = 1, 5
         accum(i) = 0
      END DO
      accept = 0
      END
C
      SUBROUTINE Metropolis(N,L,spin,E,M,w,accept)
*     one Monte Carlo step per spin
      IMPLICIT NONE
      REAL*8 E, M, w(-8:8), accept
      REAL rnd
      INTEGER dE, ispin, x, y, N, L, spin(32,32), DeltaE
      DO ispin = 1, N
         x = int(L*rnd(0) + 1)              ! random x coordinate
         y = int(L*rnd(0) + 1)              ! random y coordinate
         dE = DeltaE(x,y,L,spin)            ! compute change in energy
         IF(rnd(0) .LE. w(dE)) THEN
            spin(x,y) = -spin(x,y)          ! flip spin
            accept = accept + 1
            M = M + 2*spin(x,y)
            E = E + dE
            CALL show_change(x,y,L,spin)
         END IF
      END DO
      END
C
      FUNCTION DeltaE(x,y,L,spin)
*     periodic boundary conditions
      IMPLICIT NONE
      INTEGER DeltaE, down, left, right, up, x, y, L, spin(32,32)
      IF(x .EQ. 1) THEN
         left = spin(L,y)
      ELSE
         left = spin(x - 1,y)
      END IF
      IF(x .EQ. L) THEN
         right = spin(1,y)
      ELSE
         right = spin(x + 1,y)
      END IF
      IF(y .EQ. 1) THEN
         down = spin(x,L)
      ELSE
         down = spin(x,y - 1)
      END IF
      IF(y .EQ. L) THEN
         up = spin(x,1)
      ELSE
         up = spin(x,y + 1)
      END IF
      DeltaE = 2*spin(x,y)*(left + right + up + down)
      END
C
      SUBROUTINE show_change(i,j,L,spin)
      IMPLICIT NONE
      REAL X, Y
      INTEGER IR_, RBTN, i, j, L, spin(32,32)
      LOGICAL show
      SAVE show
      DATA show/.TRUE./
      CALL GWmouse(IR_, RBTN, X, Y)
      IF(IR_ .NE. 0) THEN
         IF(show) THEN
            show = .FALSE.
         ELSE
            show = .TRUE.
            CALL plotspins(L,spin)
         END IF
         DO WHILE(IR_ .NE. 0)
            CALL GWmouse(IR_, RBTN, X, Y)
         END DO
      ENDIF
      IF(.NOT.show) RETURN
      CALL GWerase(IR_, (i-1)*L+j, 0)
      CALL GWsetogn(IR_, 0, (i-1)*L+j)
      IF(spin(i,j) .EQ. 1) THEN
         CALL GWsrect(IR_, REAL(i), REAL(j), i+1.0, j+1.0, 13)
      ELSE
         CALL GWsrect(IR_, REAL(i), REAL(j), i+1.0, j+1.0, 14)
      END IF
      END
C
      SUBROUTINE mydata(E,M,accum)
*     accumulate data after every Monte Carlo step per spin
      IMPLICIT NONE
      REAL*8 E, M, accum(10)
      accum(1) = accum(1) + E
      accum(2) = accum(2) + E*E
      accum(3) = accum(3) + M
      accum(4) = accum(4) + M*M
      accum(5) = accum(5) + abs(M)
      END
C
      SUBROUTINE output(T,N,mcs,accum,accept)
      IMPLICIT NONE
      REAL*8 abs_mave, e2ave, eave, m2ave, mave, norm, T,
     +       accum(10), accept
      INTEGER N, mcs
      norm = 1.0D0/(mcs*N)                      ! averages per spin
      accept = accept*norm
      eave = accum(1)*norm
      e2ave = accum(2)*norm
      mave = accum(3)*norm
      m2ave = accum(4)*norm
      abs_mave = accum(5)*norm
      WRITE(*,*) 'temperature = ', T
      WRITE(*,*) 'acceptance probability = ', accept
      WRITE(*,*) 'mean energy per spin = ', eave
      WRITE(*,*) 'mean squared energy per spin = ', e2ave
      WRITE(*,*) 'mean magnetization per spin = ', mave
      WRITE(*,*) 'mean of absolute magnetization per spin = ', abs_mave
      WRITE(*,*) 'mean squared magnetization per spin = ', m2ave
      END
C
      SUBROUTINE save_config(L,T,spin)
      IMPLICIT NONE
      REAL*8 T
      INTEGER IO2_, L, x, y, spin(32,32)
      CHARACTER file*80
      DATA IO2_/12/
      WRITE(*,'(A,$)') 'name of file for last configuration = '
      READ(*,'(A)') file
      OPEN(IO2_, file = file)
      WRITE(*,*) T
      DO y = 1, L
         DO x = 1, L
            WRITE(*,*) spin(x,y)
         END DO
      END DO
      CLOSE(IO2_)
      END
C
      SUBROUTINE read_config(L,T,spin)
      IMPLICIT NONE
      REAL*8 T
      INTEGER IO1_, L, x, y, spin(32,32)
      CHARACTER file*80
      DATA IO1_/11/
      WRITE(*,'(A,$)') 'filename?'
      READ(*,'(A)') file
      OPEN(IO1_, file = file)
      READ(IO1_,*) T
      DO y = 1, L
         DO x = 1, L
            READ(IO1_,*) spin(x,y)
         END DO
      END DO
      CLOSE(IO1_)
      END
C
      SUBROUTINE correl(Ce,Cm,E,M,esave,msave,pass,nsave)
*     accumulate data for time correlation functions
*     save last nsave values of M and E
      IMPLICIT NONE
      REAL*8 Ce(0:20), Cm(0:20), E, M, esave(100), msave(100)
      INTEGER index, index0, tdiff, pass, nsave
*     index0 = array index for earliest saved time
      index0 = mod(pass-1,nsave) + 1
      IF(pass .GT. nsave) THEN
*        compute Ce and Cm after nsave values are saved
         index = index0
         DO tdiff = nsave, 1, -1
            Ce(tdiff) = Ce(tdiff) + E*esave(index)
            Cm(tdiff) = Cm(tdiff) + M*msave(index)
            index = index + 1
            IF(index .GT. nsave) index = 1
         END DO
      END IF
*     save latest value in array position of earliest value
      esave(index0) = E
      msave(index0) = M
      END
C
      SUBROUTINE c_output(Ce,Cm,accum,mcs,nsave)
      IMPLICIT NONE
      REAL*8 e2bar, ebar, m2bar, mbar, norm, Ce(0:20), Cm(0:20),
     +       accum(10)
      INTEGER tdiff, mcs, nsave
*     compute time correlation functions
      ebar = accum(1)/mcs
      e2bar = accum(2)/mcs
      Ce(0) = e2bar - ebar*ebar
      mbar = accum(3)/mcs
      m2bar = accum(4)/mcs
      Cm(0) = m2bar - mbar*mbar
      norm = 1.0D0/(mcs - nsave)
      WRITE(*,*)
      WRITE(*,'(3A13)') 't', 'Ce(t)', 'Cm(t)'
      WRITE(*,*)
      DO tdiff = 1, nsave
*        correlation functions defined so that C(t=0) = 1
*        and C(infinity) = 0
         Ce(tdiff) = (Ce(tdiff)*norm - ebar*ebar)/Ce(0)
         Cm(tdiff) = (Cm(tdiff)*norm - mbar*mbar)/Cm(0)
         WRITE(*,'(I13,2F13.6)') tdiff, Ce(tdiff), Cm(tdiff)
      END DO
      END