calculate I = \int_{-1}^{1} x^2 dx by monte carlo method in fortran

mc.f90

! ---------------------------
! calculate
!   I = \int_{-1}^{1} x^2 dx
! by monte carlo method.
! ---
! reference:
! http://aidiary.hatenablog.com/entry/20140728/1406555863
! ---------------------------

program main

  implicit none
  integer, parameter :: R_ = selected_real_kind(13)
  integer :: N = 100000
  integer :: a, b, i
  real(R_) :: rn, res
  real(R_), allocatable :: x(:)

  allocate(x(N))

  a = -1
  b = 1

  do i = 1, N

    x(i) = real((b - a)) / real(N) * real(i) - real(b)

    call random_number(rn)
    if (rn <= x(i)**2) then
      x(i) = 1
    else
      x(i) = 0
    end if

  end do

  res = real((b - a)) * sum(x(:)) / real(N)

  print *, res

end program main