Fortran program to calculate an integral using "crude" Monte Carlo

The following program calculates the integral I

using "crude" Monte Carlo. In this case, n=3, x_ia=0, x_ib=1 and f(x)=exp(a·x).

        Program intmc
        
        double precision fsum,fsum2,a,xa,xb,V,x,f,I,IA,varf
        integer M,n,j,k
        dimension a(3),xa(3),xb(3),x(3)
c n: dimension of the array x
c M: number of random numbers we use
c V: volum of the hipervolume
        n=3
        M=1000000
        V=1 
        do j=1,n
            a(j)=rand()
            xa(j)=0
            xb(j)=1
            V=V*(xb(j)-xa(j))
        enddo
c let's find the value of the integral that we get using monte-carlo
        fsum=0
c fsum2 will help us to find the error
        fsum2=0
        do k=1,M
            do j=1,n
                x(j)=rand()
            enddo
            call func(a,x,f,n)
            fsum=fsum+f
            fsum2=fsum2+f**2
        enddo
c I is the value of the integral we wanted to find
c IA is the value of the integral calculated analiticaly
        I=V*fsum/M
        varf=(fsum2/M-((fsum/M)**2))
        sigma=sqrt(varf/M)
        IA=1
        do j=1,n
            IA=IA*((exp(a(j))-1)/a(j))
        enddo
        write(*,*) I,sigma*3,IA
        stop
        end


        SUBROUTINE func(a,x,f,n)
            double precision a,x,f,ax
            integer n
            dimension a(n),x(n)
            ax=0
            do j=1,n
                ax=ax+a(j)*x(j)
            enddo
            f=exp(ax)
        end