d_simps2d_func Function

function d_simps2d_func(func,xrange,yrange,Nx,Ny) result(int)
  interface
     function func(x)
       real(8),dimension(:) :: x
       real(8)              :: func
     end function func
  end interface
  real(8),dimension(2) :: xrange,yrange
  integer :: Nx,Ny,i,j
  real(8) :: xx(2*Nx),yy(2*Ny)
  real(8) :: hx,hy
  real(8) :: int
  hx=xrange(2)-xrange(1)
  hx=hx/Nx/2
  hy=yrange(2)-yrange(1)
  hy=hy/Ny/2
  xx=linspace(xrange(1),xrange(2),2*Nx)
  yy=linspace(yrange(1),yrange(2),2*Ny)
  !
  int=&
       func([xrange(1),yrange(1)])+&
       func([xrange(1),yrange(2)])+&
       func([xrange(2),yrange(1)])+&
       func([xrange(2),yrange(2)])
  !
  do j=1,Ny
     int=int+4d0*(func([xrange(1),yy(2*j-1)])+func([xrange(2),yy(2*j-1)]))
  enddo
  do j=1,Ny-1
     int=int+2d0*(func([xrange(1),yy(2*j)])+func([xrange(2),yy(2*j)]))
  enddo
  !
  do i=1,Nx
     int=int+4d0*(func([xx(2*i-1),yrange(1)])+func([xx(2*i-1),yrange(2)]))
  enddo
  do i=1,Nx-1
     int=int+2d0*(func([xx(2*i),yrange(1)])+func([xx(2*i),yrange(2)]))
  enddo
  !
  do j=1,Ny
     do i=1,Nx
        int=int+16d0*func([xx(2*i-1),yy(2*j-1)])
     enddo
  enddo
  do j=1,Ny-1
     do i=1,Nx
        int=int+8d0*func([xx(2*i-1),yy(2*j)])
     enddo
  enddo
  !
  do j=1,Ny
     do i=1,Nx-1
        int=int+8d0*func([xx(2*i),yy(2*j-1)])
     enddo
  enddo
  do j=1,Ny-1
     do i=1,Nx-1
        int=int+4d0*func([xx(2*i),yy(2*j)])
     enddo
  enddo
  int=int*hx*hy/9d0
end function d_simps2d_func