d_simpson_nonlin_sample Function

function d_simpson_nonlin_sample(f,x) result(sum)
  real(8) :: f(:)
  real(8) :: x(size(f))
  real(8) :: dx(size(f)+1)
  real(8) :: sum,sum1,sum2,sum3
  real(8) :: a,b,dh
  integer :: i,n,m
  n=size(f)
  m=n-1
  dx=0d0
  forall(i=1:m)dx(i)=x(i+1)-x(i)
  sum1=0d0
  sum2=0d0
  sum3=0d0
  i=0
  do while(i<n)
     !Simpson's 3/8 rule
     if((dx(i)==dx(i+1)).AND.(dx(i)==dx(i+2)))then
        sum1=sum1+(3.d0*dx(i)*(f(i)+&
             3.d0*(f(i+1)+f(i+2))+f(i+3)))/8.d0
        i=i+3
        !Simpson's 1/3 rule
     elseif(dx(i)==dx(i+1))then
        sum2=sum2+(2.d0*dx(i)*(f(i)+&
             4.d0*f(i+1)+f(i+2)))/6.d0
        i=i+2
        !trapezoidal rule
     elseif(dx(i)/=dx(i+1)) then
        sum3=sum3+dx(i)*(f(i)+f(i+1))/2.d0
        i = i + 1
     endif
  enddo
  sum = sum1+sum2+sum3
end function d_simpson_nonlin_sample