d_simpson_dh_sample Function

function d_simpson_dh_sample(f,dh) result(sum)
  real(8) :: f(:)
  integer :: n
  real(8) :: dh,sum,sum1,sum2,int1,int2
  integer :: i,p,m,mm,mmm
  N=size(f)
  sum=0d0
  if(N==1)return
  sum1=0d0
  sum2=0d0
  int1=0d0
  int2=0d0
  if(mod(n-1,2)==0)then                !if n-1 is even:
     do i=2,N-1,2
        sum1 = sum1 + f(i)
     enddo
     do i=3,N-2,2
        sum2 = sum2 + f(i)
     enddo
     sum = (f(1) + 4d0*sum1 + 2d0*sum2 + f(n))*dh/3d0
  else                        !if n-1 is odd, use Simpson's for N-3 slices + 3/8rule for the last
     if (N>=6) then
        do i=2,N-4,2
           sum1 = sum1 + f(i)
        enddo
        do i=3,N-5,2
           sum2 = sum2 + f(i)
        enddo
        int1 = (f(1) + 4d0*sum1 + 2d0*sum2 + f(n-3))*dh/3d0
     endif
     int2 = (f(n-3)+3d0*f(n-2)+3d0*f(n-1)+f(n))*dh*3d0/8d0
     sum  = int1 + int2
  end if
end function d_simpson_dh_sample