program ising_monte_carlo

implicit none
integer :: n, nstep 
real :: coupling=1., temp, rand
integer, dimension (:,:), allocatable :: grid
real :: energy, energyave, sqenergyave, sigmaenergy1, sigmaenergy2
real :: magnet, magnetave, sqmagnetave, sigmamagnet1, sigmamagnet2
integer :: istep, ix, iy, de, iaccept
logical :: cont
character (len=9) :: control

open(unit=2,file='enmag.d',status='replace')
open(unit=3,file='risp.d',status='replace')
open(unit=4,file="restart.d", status="unknown", form="unformatted")
open(unit=10,file='config.d',status='replace')

call get_command_argument(1,control)
read(control,'(i4)') n
call get_command_argument(2,control)
read(control,'(f8.4)') temp
call get_command_argument(3,control)
read(control,'(i9)') nstep
call get_command_argument(4,control)
read(control,'(l)') cont

if (cont) write(unit=*,fmt=*) "Continuation run"
write(unit=*,fmt=*) "Grid size:", n
write(unit=*,fmt=*) "N. step:", nstep
write(unit=*,fmt=*) "Temperature:", temp

allocate (grid(n,n))
if (cont) then
  read(4) grid
else
!    call random_seed()
  do ix = 1,n
    do iy = 1,n
      call random_number(rand)
      if (rand.gt.0.5) then
        grid(ix,iy) = 1
      else
        grid(ix,iy) = -1
      end if
    end do
  end do
end if
energyave=0.
sqenergyave=0.
magnetave=0.
sqmagnetave=0.
iaccept=0

do istep = 1, nstep
  call random_number(rand)
  ix = int(rand*n)+1
  call random_number(rand)
  iy = int(rand*n)+1
  de = 2*coupling*grid(ix,iy)* &
                 (grid(mod(ix,n)    +1,iy) + &
                  grid(mod(ix+n-2,n)+1,iy) + &
                  grid(ix,             mod(iy,n)    +1)+ &
                  grid(ix,             mod(iy+n-2,n)+1))
  call random_number(rand)
  if ((de.lt.0).or.(rand.lt.exp(-real(de)/temp))) then
    grid(ix,iy) = -grid(ix,iy) 
    iaccept=iaccept+1
  end if
  energy=0.
  do ix=1,n
    do iy=1,n
      energy=energy-coupling*grid(ix,iy)* &
                            (grid(mod(ix,n)    +1,iy) + &
                             grid(mod(ix+n-2,n)+1,iy) + &
                             grid(ix,             mod(iy,n)    +1)+ &
                             grid(ix,             mod(iy+n-2,n)+1))
    end do
  end do
  energy=energy/(2.*n**2)
  energyave=energyave+energy
  sqenergyave=sqenergyave+energy**2
  sigmaenergy1=energyave/real(istep)
  sigmaenergy2=sqenergyave/real(istep)
  magnet=real(sum(grid))/n**2
  magnetave=magnetave+magnet
  sqmagnetave=sqmagnetave+magnet**2
  sigmamagnet1=magnetave/real(istep)
  sigmamagnet2=sqmagnetave/real(istep)
  write(2,*) istep, energy, magnet
  write(3,*) istep, (sigmaenergy2-sigmaenergy1**2)/temp**2, (sigmamagnet2-sigmamagnet1**2)/temp
end do

write(unit=*,fmt=*) "Acceptance ratio: ", iaccept/real(nstep)
write(unit=*,fmt=*) "Magnetization: ", magnetave/real(nstep)
write(unit=*,fmt=*) "Average energy: ", energyave/real(nstep)
write(unit=*,fmt=*) "Specific heat: ",  (sigmaenergy2-sigmaenergy1**2)/temp**2
write(unit=*,fmt=*) "Magnetizability: ", (sigmamagnet2-sigmamagnet1**2)/temp

do ix = 1,n
    do iy = 1,n
        write(10, '(3I5)') ix, iy, grid(ix,iy)
    end do
end do

write(4) grid

close(2)
close(3)
close(4)
close(10)
    
end program ising_monte_carlo