program p2e1
  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, IPAS, N, SEED

  ! Inicialitzem algunes variables sobre el problema
  TEMP = 1.3d0
  SEED = 234567
  MCTOT = 3000
  N = L*L

  ! 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

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

  E = ENERG(S, L, PBC)
  M = MAGNE(S, L)
  print *, "Energia inicial:", E, "Magnet. inicial:", 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
    print *, "Iter:", IMC, "Energia:", E, "Recalc.:", ENERG(S, L, PBC), "Magn:", INT(M), "Recalc.:", INT(MAGNE(S, L))
  enddo
endprogram p2e1