d_trapz2d_func Function

function d_trapz2d_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(Nx),yy(Ny)
  real(8)              :: hx,hy
  real(8)              :: int
  hx=xrange(2)-xrange(1)
  hx=hx/Nx
  hy=yrange(2)-yrange(1)
  hy=hy/Ny
  int=&
       func([xrange(1),yrange(1)])+&
       func([xrange(1),yrange(2)])+&
       func([xrange(2),yrange(1)])+&
       func([xrange(2),yrange(2)])
  xx=linspace(xrange(1),xrange(2),Nx)
  yy=linspace(yrange(1),yrange(2),Ny)
  do i=2,Nx
     do j=2,Ny
        int=int+4d0*func([xx(i),yy(j)])
     enddo
  enddo
  do j=2,Ny
     int=int+2d0*( func([xrange(1),yy(j)]) + func([xrange(2),yy(j)]) )
  enddo
  do i=2,Nx
     int=int+2d0*( func([xx(i),yrange(1)]) + func([xx(i),yrange(2)]) )
  enddo
  int=int*hx*hy/4d0
end function d_trapz2d_func