subroutine generateSpinMatrix(S, L, SEED)
  use iso_fortran_env
  implicit none
  integer*4, intent(in) :: L
  integer*4, intent(in) :: SEED
  integer*2, intent(out) :: S(1:64,1:64)
  integer*4 :: i, j
  integer(int32) :: genrand_int31

  call init_genrand(SEED)

  do i = 1, L
    do j = 1, L
      if (genrand_int31() < 2**30) then
        S(i, j) = 1
      else
        S(i, j) = -1
      endif
    enddo
  enddo
end subroutine generateSpinMatrix