d_trapz2d_sample Function

function d_trapz2d_sample(func,dhx,dhy,xrange,yrange) result(int)
  real(8),dimension(:,:)        :: func
  real(8),optional              :: dhx,dhy
  real(8),dimension(2),optional :: xrange,yrange
  real(8)                       :: hx,hy
  real(8),dimension(2)          :: xrange_,yrange_
  integer                       :: Nx,Ny,i,j
  real(8)                       :: int
  Nx=size(func,1)
  Ny=size(func,2)
  xrange_=[0d0,1d0];if(present(xrange))xrange_=xrange
  yrange_=[0d0,1d0];if(present(yrange))yrange_=yrange
  hx=xrange_(2)-xrange_(1);hx=hx/Nx
  if(present(dhx))hx=dhx
  hy=yrange_(2)-yrange_(1);hy=hy/Ny
  if(present(dhy))hy=dhy
  !
  int = func(1,1)+func(1,Ny)+func(Nx,1)+func(Nx,Ny)
  do i=2,Nx
     do j=2,Ny
        int=int+4d0*func(i,j)
     enddo
  enddo
  do j=2,Ny
     int=int+2d0*( func(1,j) + func(Nx,j) )
  enddo
  do i=2,Nx
     int=int+2d0*( func(i,1) + func(i,Ny) )
  enddo
  int=int*hx*hy/4d0
end function d_trapz2d_sample