sumup2 Subroutine

  subroutine sumup2(dx,y,inty,n)

    !! Routine to integrate a 1-D array of y values, using a process
    !! similar to Simpson's Rule, and assuming equally-spaced x values.
    !! It returns the integral at all tabulated points.
    !! 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
    !! inty(1:n) : input/output real array : calculated integral
    !! (see below)
    !! n : input integer : length of arrays y and inty
    !! This routine uses a process similar to (but not quite the same
    !! as) Simpson's Rule to integrate an array y,
    !! returning the integral up to point i in array element inty(i).
    !! Note that the first element of inty is not calculated, and must
    !! be set to the required value on entry. Usually, but not always,
    !! this value will be zero.
    !! The original source for this algorithm is not known...
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    !  Arguments

    integer, intent(in) :: n
    real(dp), intent(in) :: dx
    real(dp), intent(in), dimension(n) :: y
    real(dp), intent(inout), dimension(n) :: inty

    !  Local variables

    integer :: ix
    real(dp), parameter :: third = 1.0D0/3.0D0
    real(dp) :: thirddx
    real(dp), allocatable, dimension(:) :: yhalf

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

    allocate (yhalf(2:n))

    thirddx = third*dx

    do ix = 2,n-1
       yhalf(ix) = y(ix)-0.25D0*(y(ix+1)-y(ix-1))
    end do
    yhalf(n) = y(n-1)+0.25D0*(y(n)-y(n-2))

    do ix = 2,n
       inty(ix) = inty(ix-1) + thirddx*(y(ix)+yhalf(ix)+y(ix-1))
    end do

    deallocate (yhalf)

  end subroutine sumup2