muller Subroutine

    subroutine muller(me, Ioptn, Nparm, Numgr, Dvec, Fun, Ifun, Pttbl, Iptb, Indm, &
                      Zwork, Tolcon, Iphse, Iwork, Liwrk, Work, Lwrk, Parwrk, &
                      Err1, p1, f1, Procor, Emin)

        implicit none

        class(conmax_solver), intent(inout) :: me
        integer, intent(in) :: Ifun
        integer, intent(in) :: Indm
        integer, intent(in) :: Iptb
        integer, intent(in) :: Lwrk
        integer, intent(in) :: Nparm
        integer, intent(in) :: Numgr
        integer :: Ioptn
        integer :: Iphse
        integer :: Liwrk
        real(wp) :: Emin
        real(wp) :: f1
        real(wp) :: p1
        real(wp) :: Procor
        real(wp) :: Tolcon
        integer :: Iwork(Liwrk)
        real(wp) :: Dvec(Nparm)
        real(wp) :: Err1(Numgr + 3)
        real(wp) :: Fun(Ifun)
        real(wp) :: Parwrk(Nparm)
        real(wp) :: Pttbl(Iptb, Indm)
        real(wp) :: Work(Lwrk)
        real(wp) :: Zwork(Nparm)

        real(wp) :: acof, bcof, ccof, den, discr, f2, f3, &
                    f4, fval, p2, p3, p4, pval, temp
        integer ::  ilc08, ilc21, imain, ipmax, ismax, j, limmul, &
                   lll, nsrch

        real(wp), parameter :: tol1 = ten*ten*spcmn
        real(wp), parameter :: tol4 = tol1/four
        real(wp), parameter :: tolden = ten*spcmn

        ! SET MACHINE AND PRECISION DEPENDENT CONSTANTS.
        ilc08 = iloc(8, Nparm, Numgr)
        ilc21 = iloc(21, Nparm, Numgr)
        limmul = 5
        nsrch = 0
        imain = 0
        p3 = Procor
        f3 = Emin

        ! WE DO NOT ALLOW THE LENGTH OF THE INTERVAL (P1,P3) TO FALL BELOW
        ! TOL1.
100     if (p3 - p1 < tol1) then
            return
        else
            ! COMPUTE P2 = (P1+P3)/2.0 AND F(P2).
            p2 = (p1 + p3)/two
            pval = p2
            ! SET LLL AS THE THREAD THROUGH THE MINOTAURS CAVERN AND JUMP DOWN TO
            ! COMPUTE F(PVAL)=F(P2).  WE WILL JUMP BACK AFTER ALL SUCH JUMPS.
            lll = 1
            goto 1000
        end if

        ! HERE -TOLCON <= F2 <= TOLCON AND WE RETURN WITH PROCOR=P2 AND
        ! EMIN=F2.
200     Procor = p2
        Emin = f2
        return

        ! HERE WE HAVE NOT ACHIEVED SUCCESS YET AND WE SEE IF THE ITERATION
        ! LIMIT HAS BEEN REACHED.
300     if (nsrch < limmul) then

            ! HERE WE HAVE NOT REACHED THE ITERATION LIMIT SO WE TRY AGAIN.
            ! IF IMAIN=0 HERE WE WILL HAVE NO P4 TO SHUFFLE IN, AND WE WILL HAVE
            ! ALREADY CHECKED P3-P1 >= TOL1, SO WE RESET IMAIN TO 1 AND DO A FIT.
            if (imain <= 0) goto 900

            ! HERE WE HAVE POINTS P1, P2, P3, P4 WITH P1+TOL1/4.0 <= P2 <=
            ! P3-TOL1/4.0, P1+TOL1/4.0 <= P4 <= P3-TOL1/4.0, ABS(P4-P2) >=
            ! TOL1/4.0, F(P1) > TOLCON, F(P3) < -TOLCON, ABS(F(P2)) >
            ! TOLCON, AND ABS(F(P4)) > TOLCON.  WE WILL NOW DISCARD EITHER
            ! P1 OR P3 AND RELABEL TO GET NEW POINTS P1, P2, P3, EXCEPT IN ONE
            ! CASE WHERE TWO POINTS WILL BE DISCARDED AND WE WILL RELABEL TO GET
            ! NEW POINTS P1, P3.
            ! IF P2 > P4 HERE WE WILL, IN THE INTEREST OF A MORE READABLE
            ! PROGRAM, INTERCHANGE P2 AND P4 (AND F2 AND F4) SO WE WILL BE ABLE
            ! TO ASSUME P2 <= P4.
            if (p2 > p4) then
                temp = p2
                p2 = p4
                p4 = temp
                temp = f2
                f2 = f4
                f4 = temp
            end if
            if (f2 <= 0) then
                ! HERE F2 < 0.0.
                if (f4 <= 0) goto 700
                ! HERE F2 < 0.0 AND F4 > 0.0, AND IN THIS SAWTOOTH PATTERN WE
                ! DISCARD BOTH P4 AND P3, SET IMAIN=0, AND GO BACK TO THE BEGINNING
                ! (EXCEPT NSRCH CONTINUES TO INCREASE, INSURING EVENTUAL TERMINATION).
                imain = 0
                p3 = p2
                f3 = f2
                Procor = p3
                Emin = f3
                goto 100
            else
                ! HERE F2 > 0.0.
                if (f4 <= 0) then
                    ! HERE F2 > 0.0 AND F4 < 0.0.
                    if (p2 - p1 <= (p3 - p4)) goto 700
                end if
                ! HERE EITHER F2 > 0.0 AND F4 > 0.0, OR ELSE F2 > 0.0,
                ! F4 < 0.0, AND P2-P1 > P3-P4.  WE DISCARD P1, SINCE IN THE
                ! FORMER CASE THE FIRST THREE F VALUES ARE ALL POSITIVE, AND IN THE
                ! LATTER CASE ONLY THE FIRST TWO F VALUES ARE POSITIVE, BUT BY DROPPING
                ! P1 WE CAN GET MAXIMUM SHRINKAGE OF P3-P1.
                p1 = p2
                f1 = f2
                p2 = p4
                f2 = f4
                goto 800
            end if
        else

            ! HERE WE HAVE REACHED THE ITERATION LIMIT WITHOUT SUCCESS.  WE RETURN
            ! WITH PROCOR = THE LEFTMOST OF THE THREE POINTS P2, P4, AND P3 WHICH
            ! HAS NEGATIVE F VALUE (UNLESS IMAIN=0, IN WHICH CASE WE IGNORE P4).
            ! WRITE(NWRIT,'(A)') '***WARNING  TOO MANY ITERATIONS IN MULLER***'
            if (imain <= 0) goto 600
            if (p2 <= p4) then
                ! HERE P2 < P4.
                if (f2 < 0) goto 200
                if (f4 >= 0) goto 500
                ! HERE P4 < P2.
            else if (f4 >= 0) then
                goto 600
            end if
        end if
400     Procor = p4
        Emin = f4
        return

500     Procor = p3
        Emin = f3
        return

600     if (f2 >= 0) goto 500
        goto 200

        ! HERE EITHER F2 < 0.0 AND F4 < 0.0, OR ELSE F2 > 0.0,
        ! F4 < 0.0, AND P2-P1 <= P3-P4.  WE DISCARD P3, SINCE IN THE
        ! FORMER CASE THE LAST THREE F VALUES ARE NEGATIVE, AND IN THE LATTER
        ! CASE ONLY THE LAST TWO F VALUES ARE NEGATIVE, BUT BY DROPPING P3 WE
        ! GET MAXIMUM SHRINKAGE OF P3-P1.
700     p3 = p4
        f3 = f4

        ! HERE WE HAVE THREE POINTS.  IF P3-P1 < TOL1 WE WILL RETURN AFTER
        ! SETTING PROCOR AND EMIN.
800     if (p3 - p1 < tol1) goto 600

        ! HERE WE RESET IMAIN TO 1 AND COMPUTE P4, THE UNIQUE ZERO IN THE
        ! INTERVAL (P1,P3) OF THE QUADRATIC POLYNOMIAL WHICH PASSES THROUGH
        ! (P1,F1), (P2,F2), AND (P3,F3).  RECALL THAT F1 > 0.0,
        ! F3 < 0.0, AND P1+TOL1/4.0 <= P2 <= P3-TOL1/4.0.
900     imain = 1

        ! COMPUTE THE COEFFICIENTS ACOF, BCOF, AND CCOF OF OUR POLYNOMIAL
        ! ACOF*X**2 + BCOF*X + CCOF.
        acof = ((f3 - f2)*(p2 - p1) - (f2 - f1)*(p3 - p2))/((p2 - p1)*(p3 - p2)*(p3 - p1))
        bcof = (f3 - f1)/(p3 - p1) - acof*(p1 + p3)
        ccof = f2 - p2*(acof*p2 + bcof)
        discr = bcof**2 - four*acof*ccof
        ! IN THEORY THE DISCRIMINANT SHOULD BE POSITIVE HERE, BUT TO BE SAFE WE
        ! CHECK IT IN CASE ROUNDOFF ERROR HAS MADE IT NEGATIVE.
        if (discr < 0) goto 600
        if (bcof <= 0) then
            ! HERE BCOF <= 0.0 AND WE USE THE ALTERNATE FORM OF THE QUADRATIC
            ! FORMULA.  NOTE THAT THE DENOMINATOR CANNOT BE ZERO SINCE THAT
            ! WOULD IMPLY BOTH BCOF=0.0 AND SQRT(...)=0.0, SO ALSO EITHER
            ! ACOF=0.0 OR CCOF=0.0, BUT THIS CONTRADICTS THE FACT THAT F1  >
            ! 0.0 AND F3 < 0.0.
            ! STILL, TO BE SAFE, WE CHECK THE SIZE OF THE DENOMINATOR.
            den = -bcof + sqrt(discr)
            if (den < tolden) goto 600
            p4 = two*ccof/den
        else
            ! HERE BCOF > 0.0 AND WE USE THE USUAL FORM OF THE QUADRATIC
            ! FORMULA TO TRY TO REDUCE PROBLEMS WITH SUBTRACTION AND SMALL
            ! DENOMINATORS.  THE MINUS SIGN IS USED IN FRONT OF THE SQUARE ROOT
            ! BECAUSE IF ACOF > 0.0 THEN THE POLYNOMIAL IS CONCAVE UP, WHICH
            ! IMPLIES P1 MUST BE ON THE LEFT BRANCH (SINCE F1 > F3), WHICH
            ! IMPLIES WE WANT THE LEFT (I.E. SMALLER) ZERO, AGREEING WITH
            ! -SQRT(...)/ACOF <= 0.0.  IF ON THE OTHER HAND ACOF < 0.0 THEN
            ! THE POLYNOMIAL IS CONCAVE DOWN, WHICH IMPLIES P3 MUST BE ON THE
            ! RIGHT BRANCH (SINCE F1 > F3), WHICH IMPLIES WE WANT THE RIGHT
            ! (I.E. LARGER) ZERO, AGREEING WITH -SQRT(...)/ACOF >= 0.0.
            ! NOTE THAT ACOF=0.0 CANNOT OCCUR HERE SINCE IF IT DID THE POLYNOMIAL
            ! WOULD BE LINEAR, AND BCOF > 0.0 WOULD THEN CONTRADICT F1 > F3.
            ! STILL, TO BE SAFE, WE CHECK THE SIZE OF THE DENOMINATOR.
            den = two*acof
            if (abs(den) < tolden) goto 600
            p4 = (-bcof - sqrt(discr))/den
        end if

        ! THE NEXT SECTION (FROM HERE TO STATEMENT 3200) MODIFIES P4, IF
        ! NECESSARY, TO GET P1+TOL4 <= P2,P4 <= P3-TOL4 AND ABS(P4-P2) >=
        ! TOL4.

        ! IF ABS(P4-P2) < TOL1/4.0 WE REDEFINE P4 BY MOVING IT A DISTANCE
        ! TOL1/4.0 FROM P2 INTO THE LONGER SUBINTERVAL.  NOTE THAT THE LENGTH
        ! OF THIS SUBINTERVAL MUST BE AT LEAST TOL1/2.0 SINCE P3-P1 >= TOL1.
        if (abs(p4 - p2) < tol4) then
            if (p3 - p2 <= (p2 - p1)) then
                p4 = p2 - tol4
            else
                p4 = p2 + tol4
            end if
            ! HERE WE HAD ABS(P4-P2) >= TOL4 AND WE MAKE SURE THAT P1+TOL4
            ! <= P4 <= P3-TOL4.
        else if (p4 <= (p3 - tol4)) then
            if (p4 < (p1 + tol4)) then
                ! HERE P4 < P1+TOL4 AND WE SET P4=P1+TOL4 IF P2-P1 >= TOL1/2.0
                ! AND OTHERWISE WE SET P4=P2+TOL4.
                if (p2 - p1 < tol1/two) then
                    p4 = p2 + tol4
                else
                    p4 = p1 + tol4
                end if
            end if
            ! HERE P4 > P3-TOL4 AND WE SET P4=P3-TOL4 IF P3-P2 >= TOL1/2.0,
            ! AND OTHERWISE WE SET P4=P2-TOL4.
        else if (p3 - p2 < tol1/two) then
            p4 = p2 - tol4
        else
            p4 = p3 - tol4
        end if

        ! COMPUTE F4=F(P4).
        pval = p4
        lll = 2

        ! NOW INCREMENT NSRCH SINCE WE ARE ABOUT TO COMPUTE F.
1000    nsrch = nsrch + 1

        ! HERE IS WHERE WE MUST SUPPLY A ROUTINE TO COMPUTE FVAL = F(PVAL) =
        ! THE MAXIMUM OF THE LEFT SIDES OF THE TYPE -2 AND -1 CONSTRAINTS.

        ! PROJECT DVEC TO GET PARWRK FOR USE IN ERCMP1.
        do j = 1, Nparm
            Parwrk(j) = Zwork(j) + pval*Dvec(j)
        end do
        ! WE TAKE IPHSE=-3 AS A KLUDGE TO TELL ERCMP1 TO COMPUTE ONLY STANDARD
        ! ERRORS IF THE TEN THOUSANDS DIGIT OF IOPTN IS 1, THUS SAVING ERCMP1
        ! THE WORK OF SCANNING ICNTYP.
        call me%ercmp1(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Parwrk, 1, &
                       -3, Iwork, Liwrk, Work(ilc08), Iwork(ilc21), ipmax, ismax, &
                       Err1)
        fval = Err1(Numgr + 3)

        ! CARRY THE COMPUTED F VALUE BACK TO THE APPROPRIATE PART OF THE PROGRAM.
        if (lll == 1) then
            f2 = fval
            if (f2 > Tolcon) goto 300
            if (f2 + Tolcon >= 0) goto 200
            goto 300
        else if (lll == 2) then
            f4 = fval
            ! IF -TOLCON <= F4 <= TOLCON WE RETURN WITH PROCOR=P4 AND EMIN
            ! =F4, AND OTHERWISE WE GO BACK UP TO SEE IF WE HAVE REACHED THE LIMIT
            ! ON THE NUMBER OF STEPS.
            if (f4 > Tolcon) goto 300
            if (f4 + Tolcon >= 0) goto 400
            goto 300
        end if

    end subroutine muller