!Quasi-random Monte Carlo. Multidimensional integration routine.

subroutine MC_quasi(a, b, N, d, alfa, f_av)
use MY_STUFF
use MC_FUNCTION
Implicit none

     integer, intent(in) :: d, N
     integer :: i, j
     real(KIND=dbl) :: sumf, N_real
     real(KIND=dbl), intent(out) :: f_av
     real(KIND=dbl), dimension(d) :: p 
     real(KIND=dbl), dimension(d), intent(in) :: a, b, alfa

sumf = 0.0

do i = 1, N
      call QUASI_RANDOM_POINT(a, b, d, i, alfa, p)
      sumf = sumf + f(p)
enddo

N_real = REAL(N,dbl)

f_av = sumf/N_real

end subroutine MC_quasi