      
      
      PROGRAM ising_model
         IMPLICIT NONE
         INTEGER, PARAMETER :: lx= 100 , ly= 100
         INTEGER, PARAMETER :: no_e= 20 ! Number of replica
         INTEGER :: s(lx,ly)
         INTEGER :: i
         REAL :: t, m, x
         REAL :: mag(no_e)
         REAL :: start, finish
         
         OPEN ( UNIT= 1, FILE= 'm')
         !OPEN ( UNIT= 2, FILE= 'cv')
         !OPEN ( UNIT= 3, FILE= 'xi')
         
         CALL RANDOM_SEED
         CALL CPU_TIME(start)
         DO t= 3., 0.01, -0.1
            print *, t
            DO i= 1, no_e
               CALL init(lx, ly, s)
               !s= 1
               CALL monte_carlo(lx, ly, s, t, m)
               mag(i)= abs(m)
            ENDDO
            
            !x=(lx*ly)*(SUM(mag**2)/REAL(no_e)-(SUM(mag)/ REAL(no_e))**2)/t ! take variance from enesembles
            !print*, x, SUM(mag**2)/REAL(no_e)
            WRITE (1,*) t, SUM(mag)/ REAL(no_e)
         ENDDO
         CALL CPU_TIME(finish)
         print*,finish-start
      ENDPROGRAM ising_model
      !******************************************
      
      !******************************************
      SUBROUTINE init( lx, ly, s )
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: lx, ly
         INTEGER, INTENT(OUT) :: s(lx,ly)
         INTEGER :: i, j
         REAL :: x, y
         
         DO i= 1, lx
            DO j= 1, ly
               CALL RANDOM_NUMBER(x)
               s(i,j)= -1
               IF ( x .LT. 0.5) s(i,j)= 1
            ENDDO
         ENDDO
         RETURN
      ENDSUBROUTINE init
      !******************************************
      
      !******************************************
      SUBROUTINE monte_carlo(lx, ly, s, t, m)
         IMPLICIT NONE
         INTEGER, INTENT (IN) :: lx, ly
         INTEGER, INTENT (INOUT) :: s(lx,ly)
         REAL, INTENT(IN) :: t
         REAL, INTENT(OUT) :: m
         INTEGER*8:: n 
         REAL    :: tra_prob(-8:8)
         INTEGER :: nxp, nxm, nyp, nym, kx, ky, k
         INTEGER :: i, j
         REAL    :: x, y
         
         tra_prob= (/1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,&
         & 0.0, exp(-4.0/t), 0.0, 0.0, 0.0, exp(-8.0/t)/)
         !print*, tra_prob(-8), tra_prob(0), tra_prob(8) 
         
         n = lx* ly
         DO i= 1, n*n
            CALL RANDOM_NUMBER (x)
            CALL RANDOM_NUMBER (y)
            x= x*lx
            y= y*lx
            
            kx= FLOOR(x)+ 1 
            Ky= FLOOR(Y)+ 1

            
            CALL boundery (lx, ly, kx, ky, nxp, nxm, nyp, nym)
            
            k= s(nxp,ky)+ s(nxm,ky)+ s(kx,nyp)+ s(kx,nym)
            k= 2*k*s(kx,ky)
            
            CALL RANDOM_NUMBER(y)
            IF (y .LT. tra_prob(k)) s(kx,ky)= -s(kx,ky)
            
         ENDDO
         m= SUM(s)/REAL(lx*ly)

         
         RETURN
      ENDSUBROUTINE monte_carlo
      !******************************************
      
      !******************************************
      SUBROUTINE boundery(lx, ly, kx, ky, nxp, nxm, nyp, nym)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: lx, ly, kx, ky
         INTEGER, INTENT(OUT) :: nxp, nxm, nyp, nym
         INTEGER :: i, j
         REAL :: x, y
         !read*, kx, ky
         ! nxp ---> right neighbor
         ! nxm ---> left  neighbor
         ! nyp ---> up    neighbor
         ! nym ---> down  neighbor
         nxp= kx+ 1; nxm= kx-1
         nyp= ky+ 1; nym= ky-1
         
         IF (kx == 1 ) THEN
            nxp= 2
            nxm= lx
            
         ENDIF
         IF (kx == lx ) THEN
            nxp= 1
            nxm= lx-1
         ENDIF
         IF (ky == 1 ) THEN
            nyp= 2
            nym= ly
         ENDIF
         IF (ky == ly ) THEN
            nyp= 1
            nym= ly-1
         ENDIF
         !print*,nxp, nxm, nyp, nym
         RETURN
      ENDSUBROUTINE boundery
      !******************************************
