program ising
  implicit none
  integer, parameter :: lmax=100
  real(kind(1.d0)), parameter :: JJ = 1.d0
  real(kind(1.d0)), external :: rndm
  real(kind(1.d0)) :: Etot, Hlocal, deltaE, beta
  real(kind(1.d0)) :: avgE, avgE2, avgM, avgM2, varE, varM
  integer :: L, M, icycle, Ncycle, k, ibeta, Nbeta
  integer :: i, j, ileft, iright, jup, jdown
  integer :: sigma(0:lmax,0:lmax)
  logical :: hot_start
  
  call set_rndm_seed(9879879)
  write (*,'(a,$)') " lattice size (L) >> "
  read (*,*) L
  write (*,*) L
  if (L>lmax) stop " L > lmax"
  if (L.le.0) stop " L .le. 0"
  write (*,'(a,$)') " hot start ? (T/F)>> "
  read (*,*) hot_start
  write (*,*) hot_start
! cold start
  sigma(0:L-1,0:L-1) = 1
! hot start
  if (hot_start) then
     do i=0,L-1
        do j=0,L-1
           if (rndm() > 0.5d0) sigma(i,j) = -1
        end do
     end do
  end if
! magnetization
  call magnetization(sigma,L,lmax,M)
  write (*,*) ' Starting magnetization = ', M
! Energy
  call energy(sigma,L,lmax,JJ,Etot)
  write (*,*) ' Starting energy = ', Etot

  write (*,'(a,$)') " Ncycle, Nbeta >> "
  read(*,*) Ncycle, Nbeta
  write(*,*) Ncycle, Nbeta

  write (*,*) "       beta,    Temp,     avgM/L2,    avgE/L2,    varM/L2,    varE/L2"
  do ibeta = 1, Nbeta

  read (*,*) beta
  avgM = 0.d0
  avgM2= 0.d0
  avgE = 0.d0
  avgE2= 0.d0
  do icycle = 1, Ncycle
     do k=1,L*L
        i = L*rndm()
        j = L*rndm()
        ileft  = mod(i+1,L)
        iright = i-1; if (iright<0) iright = iright + L
        jup    = mod(j+1,L)
        jdown  = j-1; if (jdown <0) jdown  = jdown  + L
        Hlocal = JJ * ( sigma(i,jup)    + sigma(i,jdown) + &
                        sigma(iright,j) + sigma(ileft,j) )
        deltaE = + 2 * Hlocal * sigma(i,j)
        if ( rndm() < exp( -beta * deltaE) ) then
           M = M - 2 * sigma(i,j)
           Etot = Etot + deltaE
           sigma(i,j) = - sigma(i,j)
!           call check_energy(sigma,L,lmax,JJ,Etot)
        end if
     end do 
     avgM = avgM + abs(M)
     avgE = avgE + Etot
     avgM2= avgM2+ M * M
     avgE2= avgE2+ Etot * Etot
  end do 
  avgM = avgM / Ncycle
  varM = avgM2/ Ncycle - avgM * avgM
  avgE = avgE / Ncycle
  varE = avgE2/ Ncycle - avgE * avgE

  write (*,'(6f12.4)') beta, 1.d0/beta, avgM/L/L,avgE/L/L,varM/L/L,varE/L/L

  end do

  stop
end program ising

subroutine magnetization(sigma,L,lmax,M)
   implicit none
   integer, intent(in) :: L, lmax, sigma(0:lmax,0:lmax)
   integer, intent(out):: M
   integer :: i,j
   M = 0
   do i=0,L-1
      do j=0,L-1
         M = M + sigma(i,j)
      end do
   end do
   return
end subroutine magnetization

subroutine energy(sigma,L,lmax,JJ,Etot)
   implicit none
   integer, intent(in) :: L, lmax, sigma(0:lmax,0:lmax)
   real(kind(1.d0)), intent(in) :: JJ
   real(kind(1.d0)), intent(out) :: Etot
   integer :: i,j, ileft, iright,jup,jdown
   real(kind(1.d0)):: Hlocal
   Etot = 0.d0
   do i=0,L-1
      do j=0,L-1
         ileft  = mod(i+1,L)
         iright = i-1; if (iright<0) iright = iright + L
         jup    = mod(j+1,L)
         jdown  = j-1; if (jdown <0) jdown  = jdown  + L
         Hlocal = JJ * ( sigma(i,jup)    + sigma(i,jdown) + &
                         sigma(iright,j) + sigma(ileft,j) )
         Etot = Etot - 0.5 * Hlocal * sigma(i,j)
! the factor 1/2 is because pair of interacting sites are counted only once
      end do
   end do
   return
end subroutine energy 

subroutine check_energy(sigma,L,lmax,JJ,Etot)
   implicit none
   integer, intent(in) :: L, lmax, sigma(0:lmax,0:lmax)
   real(kind(1.d0)), intent(in) :: JJ
   real(kind(1.d0)), intent(in) :: Etot
   real(kind(1.d0)), parameter :: small = 1.d-6
   real(kind(1.d0)) :: Enew
  
   call energy(sigma,L,lmax,JJ,Enew)
   if (abs(Enew-Etot) < small ) return
   write (*,*) Enew, Etot
   stop ' Energy calculation is not coinsistent'

end subroutine check_energy