program main
  use, intrinsic :: iso_c_binding, only: sp=>c_float, dp=>c_double
  implicit none
  integer*4, PARAMETER :: L = 32
  integer*2 :: S(1:L, 1:L)
  real*8 :: MAGNE, ENERG, TEMP, E, DIFE, DELTA, M, SUMA, W(-8:8), genrand_real2
  integer*4 :: PBC(0:L+1), I, J, IMC, MCTOT, MCINI, MCD, IPAS, N, SEED, NSEED, SEED0
  character :: TEMPSTRING*100, str*20
  character(:), allocatable :: NOM

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

  if (command_argument_count() < 1) then
    print *, "Error: please supply a command-line argument with the"
    print *, "temperature, via './mc1 {TEMPERATURE}', where {TEMPERATURE}"
    print *, "is the temperature."
    stop
  endif

  ! Inicialitzem algunes variables sobre el problema
  call get_command_argument(1, TEMPSTRING) ! Ex: 2.0d0
  TEMPSTRING = trim(adjustl(TEMPSTRING))
  read (TEMPSTRING, *) TEMP
  NSEED = 10
  SEED0 = 117654
  MCTOT = 10000
  MCINI = 1000
  MCD = 10
  N = L*L
  NOM = "out_data/" // "SIM-L-" // trim(str(L)) // "-TEMP-" // &
      trim(str(int(1000*TEMP))) // "-MCTOT-" // trim(str(MCTOT))

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

  ! Inicialitzem variable que ens fa latent la periodicitat
  PBC(0) = L
  PBC(L + 1) = 1
  do I = 1, L
    PBC(I) = I
  enddo

  ! Cache dels valors de l'exponencial
  do I = -8, 8
    W(I) = exp(-float(I)/TEMP)
  enddo

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

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

    E = ENERG(S, L, PBC)
    M = MAGNE(S, L)
    print *, "Energia inicial:", E, "Magnet. inicial:", M
    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 = INT(genrand_real2()*L) + 1
        J = INT(genrand_real2()*L) + 1
        SUMA = S(PBC(I - 1), J) + S(PBC(I + 1), J) + S(I, PBC(J - 1)) + S(I, PBC(J + 1))
        DIFE = 2*SUMA*S(I, J)

        if (DIFE > 0) then
          DELTA = genrand_real2()
          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, 1000) == 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*SUME
  VARM = SUMM2 - SUMM*SUMM

  ! Guardem a un fiter els promitjos
  open(unit = 13, file = NOM // ".res")
  write(13, *) L, TEMP, SUMI, SUME, SUME2, VARE, SUMM, SUMAM, SUMM2, VARM

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