!
!     ...1.........2.........3.........4.........5.........6.........6.........7.........8
      program main
      implicit none
      integer(8),parameter :: n1x=50
      integer(8),parameter :: n2x=50
      integer(8),parameter :: nclasses=10
      real(8)   ,parameter :: jint=1.d0  ! spin-spin interaction
      real(8)   ,parameter :: hext=0.d0 ! external field
      real(8)   ,parameter :: kbt=2.3d0   
      integer(8),parameter :: nstep=10**7 ! number of mc steps
      real(8)   ,parameter :: tslice=1.d+5
      integer(8)           :: nmem(nclasses,n2x)
      integer(8)           :: nmem1
      integer(8)           :: state(n1x,n2x)
      integer(8)           :: class(n1x,n2x)
      integer(8)           :: i,i1,i2,istep,icount
      integer(8)           :: i1p,i1m,i2p,i2m,nn,imem,isvar
      integer(8)           :: iclass,thisclass
      integer(8)           :: nslice
      integer(8)           :: isteplast
      real(8)              :: x(nclasses)
      real(8)              :: p1(nclasses)
      real(8)              :: b
      real(8)              :: deltae(nclasses)
      real(8)              :: svar
      real(8)              :: acceptanceratio
      real(8)              :: etot,etotav
      real(8)              :: mag,magav
      real(8)              :: ran
      real(8)              :: time
      real(8)              :: deltat
      real(8)              :: profile(nclasses)
      real(8)              :: weight
!     ************************************************************************************
      weight=1.d0/real(n1x*n2x)
!
!     ====================================================================================
!     ==  set up process table                                                          ==
!     ==  each process class refers to a spin flip in a given environment               ==
!     ==  each process class is fully defined by the central state and the              ==
!     ==    average orientation of the nearest neigbors                                 ==
!     ====================================================================================
!     ==  class= 1:   central state 0; all neighbors 0                                  ==
!     ==  class= 2:   central state 0; 1 neighbor with state 1                          ==
!     ==  class= 3:   central state 0; 2 neighbors with state 1                         ==
!     ==  class= 4:   central state 0; 3 neighbors with state 1                         ==
!     ==  class= 5:   central state 0; 4 neighbors with state 1                         ==
!     ==  class= 6:   central state 1; all neighbors 0                                  ==
!     ==  class= 7:   central state 1; 1 neighbor with state 1                          ==
!     ==  class= 8:   central state 1; 2 neighbors with state 1                         ==
!     ==  class= 9:   central state 1; 3 neighbors with state 1                         ==
!     ==  class=10:   central state 1; 4 neighbors with state 1                         ==
!     ====================================================================================
      do i=1,nclasses
         isvar=int(real(i-1)/5.d0)    ! state of central atom
         nn=i-1-5*isvar               ! sum of spin up states on nearest neighbors
         deltae(i)=-2.d0*real(2*isvar-1)*(-hext-jint*real(2*nn-4))
         p1(i)=min(1,exp(-deltae(i)/kbt)) ! probability for a spin flip
      enddo
      p1(:)=p1(:)/sum(p1) ! relative jump rates for the different process classes
!
!     ====================================================================================
!     ==  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
!
!     ====================================================================================
!     ==  set up process table                                                          ==
!     ==  class assigns each site a process class                                       ==
!     ==  nmem counts the number of sites in each class in a vertical line              ==
!     ====================================================================================
      nmem(:,:)=0
      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)
          iclass=1+nn+5*state(i1,i2)
          nmem(iclass,i2)=nmem(iclass,i2)+1
          class(i1,i2)=iclass
        enddo
      enddo
!
!     ====================================================================================
!     ==  monte carlo loop                                                              ==
!     ====================================================================================
      open(10,file='nfold.dat')
      open(11,file='nfold.out')
      time=0.d0
      magav=0.d0
      etotav=0.d0
      nslice=0
      isteplast=0
      do istep=1,nstep
        if(modulo(istep,nstep/100_8).eq.0)print*,'percent finished ',100*istep/nstep
!
!       ==================================================================================
!       == select a class                                                              ==
!       ==================================================================================
        svar=0.d0
        do i=1,nclasses
          svar=svar+p1(i)*real(sum(nmem(i,:)))
          x(i)=svar
        enddo
        x=x/svar   
!       == the first class corresponds to the interval [0,x(1)],
!       == the second class corresponds to the interval [x(1),x(2)],
!       == the last class corresponds to the interval [x(n-1),x(n)], where x(n)=1.
        call random_number(ran)
        do i=1,nclasses
          if(x(i).gt.ran) then
            thisclass=i  ! process class "thisclass" is selected
            exit
          end if
        enddo
!
!       ==================================================================================
!       == select a random site from this class                                         ==
!       ==================================================================================
        nmem1=sum(nmem(thisclass,:))
        call random_number(ran)
        isvar=nint(0.5d0+ran*real(nmem1-1))  ! process isvar in this class is selected
        isvar=min(isvar,nmem1)   ! atom number in this class
        isvar=max(isvar,1)       ! atom number in this class
!
        icount=0
        do i2=1,n2x
          if(icount+nmem(thisclass,i2).lt.isvar) then
            icount=icount+nmem(thisclass,i2)
          else
            do i1=1,n1x
              if(class(i1,i2).eq.thisclass) then
                icount=icount+1
                if(icount.eq.isvar) then
                  goto 1000
                end if
              end if
            enddo
          end if
        enddo
1000    continue
!
!       ==================================================================================
!       == increment the time                                                           ==
!       ==================================================================================
        call random_number(ran)
        deltat=-log(ran)/p1(i)
        time=time+deltat
!
!       ==================================================================================
!       == update average values                                                        ==        
!       ==================================================================================
        magav=magav+mag*deltat
        etotav=etotav+etot*deltat
!
!       ==================================================================================
!       == flip spin                                                                    ==        
!       ==================================================================================
        state(i1,i2)=1-state(i1,i2)        
        mag=mag+2*(2*state(i1,i2)-1)*weight
        etot=etot+deltae(thisclass)
!
!       ==================================================================================
!       == update process list                                                          ==        
!       ==================================================================================
        i1p=1+modulo(i1,n1x)
        i1m=1+modulo(i1-2,n1x)
        i2p=1+modulo(i2,n2x)
        i2m=1+modulo(i2-2,n2x)
        nmem(class(i1,i2),i2)=nmem(class(i1,i2),i2)-1
        nmem(class(i1m,i2),i2)=nmem(class(i1m,i2),i2)-1
        nmem(class(i1p,i2),i2)=nmem(class(i1p,i2),i2)-1
        nmem(class(i1,i2m),i2m)=nmem(class(i1,i2m),i2m)-1
        nmem(class(i1,i2p),i2p)=nmem(class(i1,i2p),i2p)-1
        if(state(i1,i2).eq.1) then
          class(i1,i2)=class(i1,i2)+5
          class(i1p,i2)=class(i1p,i2)+1
          class(i1m,i2)=class(i1m,i2)+1
          class(i1,i2p)=class(i1,i2p)+1
          class(i1,i2m)=class(i1,i2m)+1
        else
          class(i1,i2)=class(i1,i2)-5
          class(i1p,i2)=class(i1p,i2)-1
          class(i1m,i2)=class(i1m,i2)-1
          class(i1,i2p)=class(i1,i2p)-1
          class(i1,i2m)=class(i1,i2m)-1
        end if
        nmem(class(i1,i2),i2)=nmem(class(i1,i2),i2)+1
        nmem(class(i1m,i2),i2)=nmem(class(i1m,i2),i2)+1
        nmem(class(i1p,i2),i2)=nmem(class(i1p,i2),i2)+1
        nmem(class(i1,i2m),i2m)=nmem(class(i1,i2m),i2m)+1
        nmem(class(i1,i2p),i2p)=nmem(class(i1,i2p),i2p)+1
!
!       ==================================================================================
!       == write result                                                                 ==
!       ==================================================================================
        if(time.gt.tslice) then
           svar=1.d0/tslice
           etotav=etotav*svar
           magav=magav*svar
           acceptanceratio=time/real(istep-isteplast)
           write(10,*)nslice*tslice,etot*weight,etotav*weight,mag,magav,acceptanceratio
           write(11,*)nslice*tslice,etot*weight,mag
           call plotstate(11,n1x,n2x,state)
           isteplast=istep
           nslice=nslice+1
           time=time-tslice
           magav=0.d0
           etotav=0.d0
        end if
      enddo
      close(10)
!
!     ====================================================================================
!     ==  analyze result                                                                ==
!     ====================================================================================
      close(11) 
      print*,'magnetization',2.d0*sum(state)/real(n1x*n2x)-1.d0
      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
