program main
  use, intrinsic :: iso_c_binding, only: sp=>c_float, dp=>c_double
  use iso_fortran_env
  implicit none
  integer*4 :: L
  real*8 :: magne, energ, TEMP, E, DIFE, DELTA, M, SUMA, W(-8:8)
  integer*4 :: I, J, IMC, MCTOT, MCINI, MCD, IPAS, N, SEED, NSEED, SEED0, pbc
  integer :: genrand_int31
  character(100) :: NOM

  integer :: SUMI
  real(dp) :: SUME, SUME2, SUMM, SUMAM, SUMM2, VARE, VARM, CV, CHI

  ! This should be replaced by L but it doesn't work dynamically, so we're
  ! setting 64 as the upper limit which should be higher than L.
  integer*2 :: S(1:64, 1:64)

  ! Inicialitzem algunes variables sobre el problema
  namelist /DADES/ L, NOM, TEMP, NSEED, SEED0, MCTOT, MCINI, MCD
  open(10, file="mc2.dat")
  read(10, DADES)
  close(10)

  N = L*L

  ! Obrim el fitxer on escriurem els resultats a cada iteració
  open(unit = 12, file = "data_out/" // trim(NOM) // ".out")

  ! Cache dels valors de l'exponencial
  do I = 0, 8
    W(I) = (2**30 - 1 + 2**30)*exp(-float(I)/TEMP)
  enddo

  SUMI = 0
  SUME = 0
  SUME2 = 0

  SUMM = 0
  SUMAM = 0
  SUMM2 = 0

  do SEED = SEED0, SEED0 + NSEED - 1
    print *, "Seed: ", SEED

    ! Inicialitzem la matriu d'spins aleatòriament
    call generateSpinMatrix(S, L, SEED)

    E = energ(S, L)
    M = magne(S, L)
    !write (12, *) "0", E, M

    ! Iterem amb el mètode de Montecarlo
    do IMC = 1, MCTOT
      do IPAS = 1, N
        ! LOS LOOPS HACEN COSAS
        I = mod(genrand_int31(), L) + 1
        J = mod(genrand_int31(), L) + 1
        SUMA = S(PBC(I - 1, L), J) + &
               S(PBC(I + 1, L), J) + &
               S(I, PBC(J - 1, L)) + &
               S(I, PBC(J + 1, L))
        DIFE = 2*SUMA*S(I, J)

        if (DIFE > 0) then
          DELTA = genrand_int31()
          if (DELTA >= W(int(DIFE))) then
            ! NO ho acceptem
            cycle
          endif
        endif

        ! Sí ho acceptem:
        S(I, J) = -S(I, J)
        E = E + DIFE
        M = M + 2*S(I, J)
      enddo

      if (mod(IMC, 10000) == 0) then
        print *, "Iter:", IMC, "Energia:", E, "Magn:", INT(M)
      endif
      !write (12, *) IMC, E, M

      if (IMC > MCINI .and. mod(IMC, MCD) == 0) then
        SUMI = SUMI + 1

        SUME = SUME + E
        SUME2 = SUME2 + E*E

        SUMM = SUMM + M
        SUMAM = SUMAM + abs(M)
        SUMM2 = SUMM2 + M*M
      endif
    enddo
  enddo

  ! Normalitzem els promitjos
  SUME = SUME/SUMI
  SUME2 = SUME2/SUMI
  SUMM = SUMM/SUMI
  SUMAM = SUMAM/SUMI
  SUMM2 = SUMM2/SUMI
  VARE = SUME2 - SUME**2
  VARM = SUMM2 - SUMM**2
  CV = (SUME2 - SUME**2)/real(N*TEMP**2)
  CHI = (SUMM2 - SUMAM**2)/real(N*TEMP)

  ! Guardem a un fiter els promitjos
  open(unit = 13, file = "data_out/" // trim(NOM) // ".res")
  !write(13, *) "L,     T,     <E>,     <E**2>,     Var(E),     <M>,     <M**2>,     Var(M),     C_V,     CHI"
  write(13, *) L, TEMP, SUME, SUME2, VARE, SUMM, SUMAM, SUMM2, VARM, CV, CHI

  ! Tanquem els fitxers de sortida
  close(12)
  close(13)
endprogram main