dpchcs Subroutine

    subroutine dpchcs (switch, n, h, slope, d, incfd, ierr)

    real(wp),intent(in)   :: switch      !! indicates the amount of control desired over
                                         !! local excursions from data.
    integer,intent(in)    :: n           !! number of data points.  (assumes N>2).
    real(wp),intent(in)   :: h(*)        !! array of interval lengths.
    real(wp),intent(in)   :: slope(*)    !! array of data slopes.
                                         !! If the data are (X(I),Y(I)), I=1(1)N, then these inputs are:
                                         !!
                                         !! * H(I) =  X(I+1)-X(I),
                                         !! * SLOPE(I) = (Y(I+1)-Y(I))/H(I),  I=1(1)N-1.
    real(wp),intent(inout) :: d(incfd,*) !! (input) array of derivative values at the data points,
                                         !! as determined by [[dpchci]].
                                         !!
                                         !! (output) derivatives in the vicinity of switches in direction
                                         !! of monotonicity may be adjusted to produce a more "visually
                                         !! pleasing" curve.
                                         !! The value corresponding to X(I) is stored in
                                         !! `D(1+(I-1)*INCFD),  I=1(1)N`
                                         !! No other entries in D are changed.
    integer,intent(in)    :: incfd       !! increment between successive values in D.
                                         !! This argument is provided primarily for 2-D applications.
    integer,intent(out)   :: ierr        !! error flag.  should be zero.
                                         !! If negative, trouble in [[dpchsw]].  (should never happen.)

    integer :: i, indx, k, nless1
    real(wp) :: del(3), dext, dfloc, dfmx, fact, slmax, wtave(2)

    real(wp),parameter :: fudge = four

    ierr = 0
    nless1 = n - 1

    ! loop over segments.
    do i = 2 , nless1
        if ( dpchst(slope(i-1),slope(i))<0 ) then
            ! slope switches monotonicity at i-th point
            ! do not change d if 'up-down-up'.
            if ( i>2 ) then
                if ( dpchst(slope(i-2),slope(i))>zero ) cycle
            endif
            if ( i<nless1 ) then
                if ( dpchst(slope(i+1),slope(i-1))>zero ) cycle
            endif
            ! compute provisional value for d(1,i).
            dext = dpchsd(slope(i-1),slope(i),h(i-1),h(i))
            ! determine which interval contains the extremum.
            if ( dpchst(dext,slope(i-1))<0 ) then
                ! dext and slope(i-1) have opposite signs -- extremum is in (x(i-1),x(i)).
                k = i - 1
                ! set up to compute new values for d(1,i-1) and d(1,i).
                wtave(2) = dext
                if ( k>1 ) wtave(1) = dpchsd(slope(k-1),slope(k),h(k-1),h(k))
            elseif ( dpchst(dext,slope(i-1))==0 ) then
                cycle
            else
                ! dext and slope(i) have opposite signs -- extremum is in (x(i),x(i+1)).
                k = i
                ! set up to compute new values for d(1,i) and d(1,i+1).
                wtave(1) = dext
                if ( k<nless1 ) wtave(2) = dpchsd(slope(k),slope(k+1),h(k),h(k+1))
            endif
        elseif ( dpchst(slope(i-1),slope(i))==0 ) then
            ! at least one of slope(i-1) and slope(i) is zero --
            ! check for flat-topped peak
            if ( i==nless1 ) cycle
            if ( dpchst(slope(i-1),slope(i+1))>=zero ) cycle
            ! we have flat-topped peak on (x(i),x(i+1)).
            k = i
            ! set up to compute new values for d(1,i) and d(1,i+1).
            wtave(1) = dpchsd(slope(k-1),slope(k),h(k-1),h(k))
            wtave(2) = dpchsd(slope(k),slope(k+1),h(k),h(k+1))
        else
            cycle
        endif

        ! at this point we have determined that there will be an extremum
        ! on (x(k),x(k+1)), where k=i or i-1, and have set array wtave--
        !    wtave(1) is a weighted average of slope(k-1) and slope(k),
        !             if k>1
        !    wtave(2) is a weighted average of slope(k) and slope(k+1),
        !             if k<n-1

        slmax = abs(slope(k))
        if ( k>1 ) slmax = max(slmax,abs(slope(k-1)))
        if ( k<nless1 ) slmax = max(slmax,abs(slope(k+1)))

        if ( k>1 ) del(1) = slope(k-1)/slmax
        del(2) = slope(k)/slmax
        if ( k<nless1 ) del(3) = slope(k+1)/slmax

        if ( (k>1) .and. (k<nless1) ) then
            ! normal case -- extremum is not in a boundary interval.
            fact = fudge*abs(del(3)*(del(1)-del(2))*(wtave(2)/slmax))
            d(1,k) = d(1,k) + min(fact,one)*(wtave(1)-d(1,k))
            fact = fudge*abs(del(1)*(del(3)-del(2))*(wtave(1)/slmax))
            d(1,k+1) = d(1,k+1) + min(fact,one)*(wtave(2)-d(1,k+1))
        else
            ! special case k=1 (which can occur only if i=2) or
            ! k=nless1 (which can occur only if i=nless1).
            fact = fudge*abs(del(2))
            d(1,i) = min(fact,one)*wtave(i-k+1)
            ! note that i-k+1 = 1 if k=i  (=nless1),
            ! i-k+1 = 2 if k=i-1(=1).
        endif
        ! adjust if necessary to limit excursions from data.
        if ( switch>zero ) then
            dfloc = h(k)*abs(slope(k))
            if ( k>1 ) dfloc = max(dfloc,h(k-1)*abs(slope(k-1)))
            if ( k<nless1 ) dfloc = max(dfloc,h(k+1)*abs(slope(k+1)))
            dfmx = switch*dfloc
            indx = i - k + 1
            ! indx = 1 if k=i, 2 if k=i-1.
            call dpchsw(dfmx,indx,d(1,k),d(1,k+1),h(k),slope(k),ierr)
            if ( ierr/=0 ) return
        endif
        ! end of segment loop.
    enddo

    end subroutine dpchcs