MODULE Shared
  USE, INTRINSIC :: ISO_FORTRAN_ENV, dp=>REAL64 !modern DOUBLE PRECISION
  INTEGER :: i, j, k, L, m, n_lattice, up, down, left, right, steps, MCS, MCS_read, equilibrium_read, equilibrium
  INTEGER, PARAMETER :: NEAREST_NEIGHBORS = 4, lattice = 100, spins_dat = 10 !manage UNIT numbers as a resource
  REAL(dp) :: spin_max, spin_old, beta, E_i, E_f, dE, dE_max, random, theta2, avg_theta2, coupling_strength, coupling_strength_read
  REAL(dp), DIMENSION(:, :), ALLOCATABLE :: spin
  REAL(dp), DIMENSION(:), ALLOCATABLE :: boltzmann
  REAL(dp), PARAMETER :: PI = ATAN2(0.0_dp, -1.0_dp), delta = 1.0_dp, T = 0.1_dp
  CHARACTER(*), PARAMETER :: my_format = '(2f22.15)' !.15 is dp
END MODULE Shared

PROGRAM kt1
  USE Shared
  IMPLICIT NONE
  CALL INIT_RANDOM_SEED()
  WRITE(6, '(a29)', advance='no') "Number of Monte Carlo steps: "
  READ *, MCS_read
  MCS = MCS_read
  WRITE(6, '(a29)', advance='no') "Number of equilibrium steps: "
  READ *, equilibrium_read
  equilibrium = equilibrium_read
  WRITE(6, '(a32)', advance='no') "Value of the coupling strength: "
  READ *, coupling_strength_read
  coupling_strength = coupling_strength_read
  ALLOCATE(spin(lattice, lattice), boltzmann(MCS * lattice))
  spin_max = 2.0_dp * PI
  n_lattice = NEAREST_NEIGHBORS !initial (minimum) lattice size
  beta = 1 / T !Boltzmann's constant = 1 in this simulation
  OPEN(spins_dat, file='spins.dat', access='sequential', status='unknown', action='write')
  theta2 = 0.0_dp
  spin = 0.0_dp !zero entire array
  boltzmann = 0.0_dp !zero entire array
  compute: DO
    IF (n_lattice > lattice) EXIT
    dE_max = 8.0_dp * coupling_strength !compute energy resolution
    dE = 0.0_dp
    resolution: DO m = 1, INT(1000.0_dp) !energy resolution of 0.001
      boltzmann(m) = EXP(-dE * beta)
      dE = dE + dE_max / 1000.0_dp
    END DO resolution
    random_x: DO i = 1, n_lattice
      random_y: DO j = 1, n_lattice
        CALL RANDOM_NUMBER(random)
        spin(i, j) = 2.0_dp * PI * random
      END DO random_y
    END DO random_x
    steps = MCS + equilibrium
    spin_lattice: DO L = 1, steps
      CALL RANDOM_NUMBER(random)
      i = INT(FLOAT(n_lattice) * random) + 1 !select a spin at random
      CALL x(i, n_lattice, left, right)
      CALL RANDOM_NUMBER(random)
      j = INT(FLOAT(n_lattice) * random) + 1 !select a spin at random
      CALL y(j, n_lattice, up, down)
      E_i = -coupling_strength * (COS(spin(right, j) - spin(i, j)) + COS(spin(left, j) - spin(i, j)) + COS(spin(i, up)             &
        - spin(i, j)) + COS(spin(i, down) - spin(i, j))) !calculate initial energy
      spin_old = spin(i, j) !preserve old spin
      CALL RANDOM_NUMBER(random)
      spin(i, j) = spin(i, j) + (1.0_dp - 2.0_dp * random) * delta !change the preserved spin by a random amount
      IF (spin(i, j) > spin_max) THEN !apply periodic boundary conditions
        spin(i, j) = spin(i, j) - spin_max
      ELSE IF (spin(i, j) < 0.0_dp) THEN
        spin(i, j) = spin(i, j) + spin_max
      END IF
      E_f = -coupling_strength * (COS(spin(right, j) - spin(i, j)) + COS(spin(left, j) - spin(i, j)) + COS(spin(i, up)             &
        - spin(i, j)) + COS(spin(i, down) - spin(i, j))) !calculate final energy
      dE = E_f - E_i
      k = INT(1000.0_dp * (dE / dE_max)) + 1
      IF (dE >= 0.0_dp) THEN !accept microstate
        CALL RANDOM_NUMBER(random)
        IF (boltzmann(k) <= random) THEN !reject microstate; reset energy, spin
          E_f = E_i
          spin(i, j) = spin_old
        END IF
      END IF
      IF (L >= equilibrium) theta2 = theta2 + spin(i, j) * spin(i, j) !running averages
    avg_theta2 = theta2 / FLOAT(MCS)
    END DO spin_lattice
    WRITE(spins_dat, my_format) LOG(FLOAT(n_lattice)), avg_theta2
    n_lattice = n_lattice + 1 !expand lattice
  END DO compute
  DEALLOCATE(spin, boltzmann)
  CLOSE(spins_dat)
END PROGRAM kt1

SUBROUTINE x(is, n_latts, lefts, rights) !special case nearest neighbors in x-direction
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: is, n_latts
  INTEGER, INTENT(OUT) :: lefts, rights
  IF (is == n_latts) THEN
    lefts = n_latts - 1
    rights = 1
  ELSE IF (is == 1) THEN
    lefts = n_latts
    rights = 2
  ELSE
    lefts = is - 1
    rights = is + 1
  END IF
END SUBROUTINE x

SUBROUTINE y(js, n_latts, ups, downs) !special case nearest neighbors in y-direction
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: js, n_latts
  INTEGER, INTENT(OUT) :: ups, downs
  IF (js == n_latts) THEN
    ups = 1
    downs = n_latts - 1
  ELSE IF (js == 1) THEN
    ups = 2
    downs = n_latts
  ELSE
    ups = js + 1
    downs = js - 1
  END IF
END SUBROUTINE y

SUBROUTINE INIT_RANDOM_SEED()
  USE ISO_Fortran_env, ONLY: INT64
  IMPLICIT NONE
  INTEGER, ALLOCATABLE :: seed(:)
  INTEGER :: i, n, un, istat, dt(8), pid
  INTEGER(INT64) :: t
  CALL RANDOM_SEED(size = n)
  ALLOCATE(seed(n))
  OPEN(newunit=un, file='/dev/urandom', access='stream', status='old', action='read', form='unformatted', iostat=istat)
  IF (istat == 0) THEN
    READ(un) seed
    CLOSE(un)
  ELSE
    CALL SYSTEM_CLOCK(t)
    IF (t == 0) THEN
      CALL DATE_AND_TIME(values = dt)
      t = (dt(1) - 1970) * 365_INT64 * 24 * 60 * 60 * 1000 + dt(2) * 31_INT64 * 24 * 60 * 60 * 1000 + dt(3) * 24_INT64 * 60 * 60   &
        * 1000 + dt(5) * 60 * 60 * 1000 + dt(6) * 60 * 1000 + dt(7) * 1000 + dt(8)
    END IF
    pid = GETPID()
    t = IEOR(t, INT(pid, KIND(t)))
    DO i = 1, n
      seed(i) = lcg(t)
    END DO
  END IF
  CALL RANDOM_SEED(put = seed)
  DEALLOCATE(seed)
CONTAINS
  FUNCTION lcg(s)
    INTEGER :: lcg
    INTEGER(INT64) :: s
    IF (s == 0) THEN
      s = 104729
    ELSE
      s = MOD(s, 4294967296_INT64)
    END IF
    s = MOD(s * 279470273_INT64, 4294967291_INT64)
    lcg = INT(MOD(s, INT(HUGE(0), INT64)), KIND(0))
  END FUNCTION lcg
END SUBROUTINE INIT_RANDOM_SEED