bndset Subroutine

    subroutine bndset(Nparm, x, Itersl, Numin, Prjslp, Cofbnd, Xkeep, Bndkp)

        implicit none

        integer, intent(in)                           :: Nparm
        real(wp), dimension(Nparm + 1), intent(in)    :: x
        integer, intent(in)                           :: Itersl
        integer, intent(in)                           :: Numin
        real(wp), intent(in)                          :: Prjslp
        real(wp), dimension(Nparm), intent(inout)     :: Cofbnd
        real(wp), dimension(Nparm + 1), intent(inout) :: Xkeep
        real(wp), dimension(Nparm), intent(inout)     :: Bndkp

        ! set initial parameters.  fact1, fact3a, fact3b, chlm1, and chlm2
        ! should be between 0.0 and 1.0, while fact2 should be > 1.0.
        real(wp), parameter :: fact1 = (one + two)/four
        real(wp), parameter :: fact2 = two
        real(wp), parameter :: fact3a = one/ten
        real(wp), parameter :: fact3b = one/(ten*ten)
        real(wp), parameter :: fact4 = two/ten
        real(wp), parameter :: chlm1 = one/ten
        real(wp), parameter :: chlm2 = (four + four)/ten
        real(wp), parameter :: tstprj = one/two - one/(ten*ten*ten)
        real(wp), parameter :: epsil = ten*ten*spcmn
        real(wp), parameter :: epsil1 = (one + one/(ten*ten*ten))*epsil
        real(wp), parameter :: bnd = two/(ten*ten) !! The initial bound on all coefficient changes.

        real(wp) :: bsave, fact3
        integer :: itight, j

        if (Numin < 1) then

            if (Itersl < 1) then

            else if (Itersl == 1) then
                ! HERE NUMIN=0 AND ITERSL=1, SO THE LAST BNDSET CALL RESULTED IN
                ! THE FIRST SUCCESSFUL PRINCIPAL ERROR NORM IMPROVEMENT,
                ! AND SO WE SAVE COFBND IN BNDKP AND X IN XKEEP.  WE WILL NOT
                ! CHANGE COFBND HERE.
                do j = 1, Nparm
                    Xkeep(j) = x(j)
                    Bndkp(j) = Cofbnd(j)
                end do
                return
            else
                ! HERE NUMIN=0 AND ITERSL >= 2, SO WE HAVE HAD AT LEAST 2 SUCCESSES,
                ! WITH THE COEFFICIENTS AND BOUNDS FOR THE LAST ONE IN X AND
                ! COFBND RESPECTIVELY, AND THE COEFFICIENTS AND BOUNDS FOR THE
                ! PREVIOUS ONE IN XKEEP AND BNDKP RESPECTIVELY.  WE WILL FORM A
                ! NEW COFBND, AND SHIFT THE OLD COFBND INTO BNDKP AND X INTO XKEEP.
                do j = 1, Nparm
                    Cofbnd_block: block
                        ! SAVE THE OLD COFBND(J) IN BSAVE.
                        bsave = Cofbnd(j)
                        ! IF AT BOTH THE LAST AND PREVIOUS SUCCESSFUL ITERATION THE CHANGES
                        ! IN A COEFFICIENT RELATIVE TO ITS BOUND WERE >= CHLM2 IN ABSOLUTE
                        ! VALUE AND IN THE SAME DIRECTION, WE LOOSEN THE BOUND BY A FACTOR
                        ! OF FACT2.  IF THE RELATIVE CHANGES WERE >= CHLM1 IN ABSOLUTE
                        ! VALUE AND IN OPPOSITE DIRECTIONS, WE TIGHTEN THE BOUND BY A FACTOR
                        ! OF FACT1 BECAUSE OF SUSPECTED OSCILLATION.  WE ALSO TIGHTEN THE
                        ! BOUND IF BOTH RELATIVE CHANGES WERE LESS THAN CHLM1 IN ABSOLUTE
                        ! VALUE IN ORDER TO PREVENT A LONG SEQUENCE OF OSCILLATIONS OF THE
                        ! SAME SMALL ORDER.  OTHERWISE WE LEAVE THE BOUND ALONE.
                        ! THE NEXT FOUR IF STATEMENTS CHECK TO SEE IF THE BOUND SHOULD BE
                        ! LOOSENED.
                        if (x(j) < chlm2*Cofbnd(j)) then
                            if (x(j) + chlm2*Cofbnd(j) <= 0) then
                                if (Xkeep(j) + chlm2*Bndkp(j) <= 0) then
                                    ! LOOSEN THE BOUND.
                                    Cofbnd(j) = fact2*Cofbnd(j)
                                    exit Cofbnd_block
                                end if
                            end if
                        else if (Xkeep(j) >= chlm2*Bndkp(j)) then
                            Cofbnd(j) = fact2*Cofbnd(j)
                            exit Cofbnd_block
                        end if
                        ! HERE THE BOUND SHOULD NOT BE LOOSENED.  THE NEXT FIVE IF
                        ! STATEMTENTS CHECK TO SEE IF IT SHOULD BE TIGHTENED.
                        if (x(j) < chlm1*Cofbnd(j)) then
                            if (x(j) + chlm1*Cofbnd(j) <= 0) then
                                if (Xkeep(j) < chlm1*Bndkp(j)) exit Cofbnd_block
                                ! HERE WE HAVE ABS(X(J)) < CHLM1*COFBND(J).
                            else if (abs(Xkeep(j)) >= chlm1*Bndkp(j)) then
                                exit Cofbnd_block
                            end if
                        else if (Xkeep(j) + chlm1*Bndkp(j) > 0) then
                            exit Cofbnd_block
                        end if
                        ! TIGHTEN THE BOUND.
                        Cofbnd(j) = fact1*Cofbnd(j)
                        ! DO NOT ALLOW THE BOUND TO DROP BELOW EPSIL.
                        if (Cofbnd(j) < epsil) Cofbnd(j) = epsil
                    end block Cofbnd_block
                    ! SAVE X(J) AND THE OLD COFBND(J).
                    Bndkp(j) = bsave
                    Xkeep(j) = x(j)
                end do
                ! IF THE LAST PROJECTION FACTOR IS SMALLER THAN .499, WE TIGHTEN THE
                ! BOUNDS BY A FACTOR OF 0.2, WITH THE RESTRICTION THAT WE DO NOT
                ! ALLOW THE BOUNDS TO DROP BELOW EPSIL.
                if (Prjslp < tstprj) then
                    do j = 1, Nparm
                        Cofbnd(j) = fact4*Cofbnd(j)
                        if (Cofbnd(j) < epsil) Cofbnd(j) = epsil
                    end do
                end if
                return
            end if

            ! HERE NUMIN=0 AND ITERSL=0, SO WE ARE IN THE FIRST BNDSET CALL SINCE THE
            ! LAST RK SUCCESS (IF ANY), SO WE SET INITIAL BOUNDS.
            Cofbnd = bnd
            return

        else if (Numin == 1) then
            ! HERE NUMIN=1 SO THE LAST BNDSET CALL RESULTED IN FAILURE TO
            ! IMPROVE THE PRINCIPAL ERROR NORM, AND WE SET FACT3=
            ! FACT3A AND TIGHTEN THE BOUNDS.
            fact3 = fact3a
        else
            ! HERE NUMIN > 1 SO THERE HAVE BEEN AT LEAST 2 SUCCESSIVE
            ! FAILURES, AND WE SET FACT3=FACT3B AND TIGHTEN THE BOUNDS.
            fact3 = fact3b
        end if

        ! TIGHTEN THE BOUNDS BY A FACTOR OF FACT3.
        itight = 1
        do j = 1, Nparm
            bsave = Cofbnd(j)
            Cofbnd(j) = fact3*bsave
            ! WE DO NOT ALLOW A BOUND TO DROP BELOW EPSIL.
            if (Cofbnd(j) < epsil) then
                ! IF THE BOUND WAS ALREADY (ESSENTIALLY) AT EPSIL, KEEP TRACK OF
                ! THIS BY NOT SETTING ITIGHT=0.
                if (bsave > epsil1) itight = 0
                Cofbnd(j) = epsil
            else
                itight = 0
            end if
        end do

        ! IF ALL THE BOUNDS WERE ALREADY (ESSENTIALLY) AT EPSIL, WE TRY
        ! RESETTING THE BOUNDS TO THEIR ORIGINAL VALUES.
        !WRITE(NWRIT,'(A)') '*****RESETTING BOUNDS TO THEIR ORIGINAL VALUES*****'
        if (itight > 0) Cofbnd = bnd

    end subroutine bndset