dpchsw Subroutine

    subroutine dpchsw (dfmax, iextrm, d1, d2, h, slope, ierr)

    integer,intent(in)  :: iextrm       !! index of the extreme data value.  (assumes
                                        !! IEXTRM = 1 or 2 .  Any value /=1 is treated as 2.)
    integer,intent(out) :: ierr         !! error flag.  should be zero.
                                        !!
                                        !! * If IERR=-1, assumption on D1 and D2 is not satisfied.
                                        !! * If IERR=-2, quadratic equation locating extremum has
                                        !!   negative discriminant (should never occur).
    real(wp),intent(in) :: dfmax        !! maximum allowed difference between F(IEXTRM) and
                                        !! the cubic determined by derivative values D1,D2.  (assumes
                                        !! DFMAX>0.)
    real(wp),intent(inout) :: d1        !! derivative values at the ends of the interval.
                                        !! (Assumes D1*D2 <= 0.)
                                        !! (output) may be modified if necessary to meet the restriction
                                        !! imposed by DFMAX.
    real(wp),intent(inout) :: d2        !! derivative values at the ends of the interval.
                                        !! (Assumes D1*D2 <= 0.)
                                        !! (output) may be modified if necessary to meet the restriction
                                        !! imposed by DFMAX.
    real(wp),intent(in) :: h            !! interval length.  (Assumes  H>0.)
    real(wp),intent(in) :: slope        !! data slope on the interval.

    real(wp) :: cp, hphi, nu, radcal, sigma
    real(wp) :: lambda  !! the ratio of d2 to d1.
    real(wp) :: phi     !! phi is the normalized value of p(x)-f1 at x = xhat = x-hat(rho),
                        !! where  that = (xhat - x1)/h .
                        !! that is, p(xhat)-f1 = d*h*phi,  where d=d1 or d2.
                        !! similarly,  p(xhat)-f2 = d*h*(phi-rho) .
    real(wp) :: rho     !! the ratio of the data slope to the derivative being tested.
    real(wp) :: that    !! that = t-hat(rho) is the normalized location of the extremum.

    real(wp),parameter :: fact = 100.0_wp
    real(wp),parameter :: third = one/three - d1mach4  !! third should be slightly less than 1/3 (original code had 0.33333)
    real(wp),parameter :: small = fact*d1mach4 !! small should be a few orders of magnitude greater than macheps.

    !  do main calculation.

    if (d1 == zero) then

        ! special case -- d1==zero .
        ! if d2 is also zero, this routine should not have been called.
        if (d2 == zero) then
            ! d1 and d2 both zero, or both nonzero and same sign.
            ierr = -1
            call xermsg ('PCHIP', 'dpchsw', 'd1 and/or d2 invalid', ierr, 1)
            return
        end if

        rho = slope/d2
        ! extremum is outside interval when rho >= 1/3 .
        if (rho >= third) then
            ierr = 0
            return
        end if
        that = (two*(three*rho-one)) / (three*(two*rho-one))
        phi = that**2 * ((three*rho-one)/three)

        ! convert to distance from f2 if iextrm/=1 .
        if (iextrm /= 1) phi = phi - rho

        ! test for exceeding limit, and adjust accordingly.
        hphi = h * abs(phi)
        if (hphi*abs(d2) > dfmax) then
            ! at this point, hphi>0, so divide is ok.
            d2 = sign (dfmax/hphi, d2)
        endif

    else

        rho = slope/d1
        lambda = -d2/d1
        if (d2 == zero) then
            ! special case -- d2==zero .
            ! extremum is outside interval when rho >= 1/3 .
            if (rho >= third) then
                ierr = 0
                return
            end if
            cp = two - three*rho
            nu = one - two*rho
            that = one / (three*nu)

        else
            if (lambda <= zero) then
                ! d1 and d2 both zero, or both nonzero and same sign.
                ierr = -1
                call xermsg ('PCHIP', 'dpchsw', 'd1 and/or d2 invalid', ierr, 1)
                return
            end if

            ! normal case -- d1 and d2 both nonzero, opposite signs.

            nu = one - lambda - two*rho
            sigma = one - rho
            cp = nu + sigma
            if (abs(nu) > small) then
                radcal = (nu - (two*rho+one))*nu + sigma**2
                if (radcal < zero) then
                    ! negative value of radical (should never occur).
                    ierr = -2
                    call xermsg ('PCHIP', 'dpchsw', 'negative radical', ierr, 1)
                end if
                that = (cp - sqrt(radcal)) / (three*nu)
            else
                that = one/(two*sigma)
            endif
        endif
        phi = that*((nu*that - cp)*that + one)

        ! convert to distance from f2 if iextrm/=1 .
        if (iextrm /= 1) phi = phi - rho

        ! test for exceeding limit, and adjust accordingly.
        hphi = h * abs(phi)
        if (hphi*abs(d1) > dfmax) then
            ! at this point, hphi>0, so divide is ok.
            d1 = sign (dfmax/hphi, d1)
            d2 = -lambda*d1
        endif

    endif

    ierr = 0

    end subroutine dpchsw