!
!     ...1.........2.........3.........4.........5.........6.........6.........7.........8
      program main
      implicit none
      integer(8),parameter :: n1x=50
      integer(8),parameter :: n2x=50
      real(8)   ,parameter :: jint=1.d0  ! spin-spin interaction
      real(8)   ,parameter :: hext=-0.d-2  ! external field
      real(8)   ,parameter :: kbt=2.3d0
      integer(8),parameter :: nstep=10**8  ! number of mc steps
      integer(8),parameter :: nsteppersample=nstep/1000
      integer(8)           :: state(n1x,n2x)
      integer(8)           :: i1,i2,istep
      integer(8)           :: isite,i1p,i1m,i2p,i2m,nn
      real(8)              :: mag,magdirect,magav ! magnetization/per site
      real(8)              :: mag2av,magfluct
      real(8)              :: etot,etotdirect,etotav
      real(8)              :: weight
      real(8)              :: deltae
      real(8)              :: ran
      real(8)              :: acceptanceratio
      real(8)              :: svar              !support variable
      logical              :: tflip
      character(n2x)       :: string
!     ************************************************************************************
!     ==  state=0 site is vacant
!     ==  state=1 site is occupied
      weight=1.d0/real(n1x*n2x)
!     ====================================================================================
!     ==  prepare initial state                                                         ==
!     ====================================================================================
      do i1=1,n1x
        do i2=1,n2x
          call random_number(ran)
          state(i1,i2)=nint(ran)
        enddo
      enddo
!
!     ====================================================================================
!     ==  calculate initial properties                                                  ==
!     ====================================================================================
      call totalenergy(n1x,n2x,hext,jint,state,etot,mag)
      print*,'magnetization',mag
      acceptanceratio=0.d0
      etotav=0.d0
      magav=0.d0
!
!     ====================================================================================
!     ==  monte carlo loop                                                              ==
!     ====================================================================================
      open(10,file='ising2d.dat')
      open(11,file='ising2d.out')
      do istep=1,nstep
!
!       ==================================================================================
!       == select a random site                                                         ==
!       ==================================================================================
        call random_number(ran)
        i1=nint(0.5d0+ran*real(n1x))
        call random_number(ran)
        i2=nint(0.5d0+ran*real(n2x))
!
!       ==================================================================================
!       == determine enery difference for the step                                      ==       
!       ==================================================================================
        i1p=1+modulo(i1,n1x)
        i1m=1+modulo(i1-2,n1x)
        i2p=1+modulo(i2,n2x)
        i2m=1+modulo(i2-2,n2x)
        nn=state(i1p,i2)+state(i1m,i2)+state(i1,i2p)+state(i1,i2m)
!       == calculate energy change for a spin flip ======================================
        deltae=-2.d0*real(2*state(i1,i2)-1)*(-hext-jint*real(2*nn-4))
!
!       ==================================================================================
!       == choose acceptance and flip spin                                              ==
!       ==================================================================================
        if(deltae.le.0.d0) then
          tflip=.true.
        else
          call random_number(ran)
          tflip=ran.lt.exp(-deltae/kbt)
        end if
        if(tflip) then
          state(i1,i2)=1-state(i1,i2)
          acceptanceratio=acceptanceratio+1
          etot=etot+deltae
          mag=mag+real(2*(2*state(i1,i2)-1))*weight
        end if
        magav=magav+mag
        mag2av=mag2av+mag**2
        etotav=etotav+etot
!        
!       ==================================================================================
!       == analysis                                                                     ==
!       ==================================================================================
        if(modulo(istep,nsteppersample).eq.0) then
!         == form time average ===========================================================
          svar=1.d0/real(nsteppersample)
          acceptanceratio=acceptanceratio*svar
          magav=magav*svar
          mag2av=mag2av*svar
          etotav=etotav*svar
!         == print result
!          magfluct=sqrt(mag2av-mag**2)
!         call totalenergy(n1x,n2x,hext,jint,state,etotdirect,magdirect)
!         write(10,*)istep,etot*weight,etotdirect*weight,acceptanceratio,mag,magdirect
          write(10,*)istep,etotav*weight,magav,acceptanceratio
          write(11,*)istep,etot*weight,mag
          call plotstate(11,n1x,n2x,state)
          acceptanceratio=0
          magav=0.d0
          etotav=0.d0
        end if
      enddo
      close(10)
!
!     ====================================================================================
!     ==  analyze result                                                                ==
!     ====================================================================================
      print*,'magnetization',mag
      close(11)
      stop
      end
!
!     ...1.........2.........3.........4.........5.........6.........6.........7.........8
      subroutine totalenergy(n1x,n2x,hext,jint,state,etot,mag)
      implicit none
      integer(8),intent(in) :: n1x
      integer(8),intent(in) :: n2x
      real(8)   ,intent(in) :: hext
      real(8)   ,intent(in) :: jint
      integer(8),intent(in) :: state(n1x,n2x)
      real(8)   ,intent(out):: etot
      real(8)   ,intent(out):: mag
      integer(8)            :: i1,i2
      integer(8)            :: i1p,i1m,i2p,i2m
      integer(8)            :: nn
      real(8)               :: hint
!     ************************************************************************************
      etot=0.d0 
      mag=0.d0
      do i1=1,n1x
        do i2=1,n2x
          i1p=1+modulo(i1,n1x)
          i1m=1+modulo(i1-2,n1x)
          i2p=1+modulo(i2,n2x)
          i2m=1+modulo(i2-2,n2x)
          nn=state(i1p,i2)+state(i1m,i2)+state(i1,i2p)+state(i1,i2m)
          hint=jint*real(2*nn-4)    ! sum_k jint*sigma_k
          etot=etot-(hext+0.5d0*hint)*real(2*state(i1,i2)-1)
          mag=mag+real(2*state(i1,i2)-1)
        enddo
      enddo
      mag=mag/real(n1x*n2x)
      return
      end
!
!     ...1.........2.........3.........4.........5.........6.........6.........7.........8
      subroutine plotstate(nfil,n1x,n2x,state)
!     ************************************************************************************
      implicit none
      integer  ,intent(in) :: nfil
      integer(8),intent(in) :: n1x
      integer(8),intent(in) :: n2x
      integer(8),intent(in) :: state(n1x,n2x)
      integer(8)            :: i1,i2
      character(n2x)        :: string
!     ************************************************************************************
      do i2=1,n2x
        string(i2:i2)='-'
      enddo
      write(nfil,*)'|'//string//'|'
      do i1=1,n1x
        string=''
        do i2=1,n2x
          if(state(i1,i2).eq.1)string(i2:i2)='o'
        enddo
        write(nfil,*)'|'//string//'|'
      enddo
      do i2=1,n2x
        string(i2:i2)='-'
      enddo
      write(nfil,*)'|'//string//'|'
      return
      end
