sumup3 Subroutine

  subroutine sumup3(dx,y,integral,n)

    !! Routine to integrate a 1-D array of y values using the
    !! Extended Simpson's Rule, assuming equally-spaced x values
    !! author: P J Knight, CCFE, Culham Science Centre
    !! dx : input real : (constant) spacing between adjacent x values
    !! y(1:n) : input real array : y values to be integrated
    !! integral : output real : calculated integral
    !! n : input integer : length of array y
    !! This routine uses Simpson's Rule to integrate an array y.
    !! If n is even, routine <CODE>sumup2</CODE> is called to
    !! perform the calculation.
    !! <P>Note: unlike sumup1 and sumup2, this routine returns only
    !! the complete integral, not the intermediate values as well.
    !! None
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    !  Arguments

    integer, intent(in) :: n
    real(dp), intent(in) :: dx
    real(dp), intent(in), dimension(n) :: y
    real(dp), intent(out) :: integral

    !  Local variables

    integer :: ix
    real(dp) :: sum1
    real(dp), allocatable, dimension(:) :: inty

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if (mod(n,2) == 0) then

       !  Use sumup2 if the number of tabulated points is even

       allocate (inty(n))

       inty(1) = 0.0D0
       call sumup2(dx,y,inty,n)
       integral = inty(n)

       deallocate (inty)

    else

       sum1 = y(1)
       do ix = 2,n-3,2
          sum1 = sum1 + 4.0D0*y(ix) + 2.0D0*y(ix+1)
       end do
       integral = dx/3.0D0*(sum1 + 4.0D0*y(n-1) + y(n))

    end if

  end subroutine sumup3