corrct Subroutine

    subroutine corrct(me, Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, &
                      Icntyp, Unit, Tolcon, Rchin, Error, Mact, Iact, Projct, &
                      Iphse, Iwork, Liwrk, Work, Lwrk, Parwrk, Err1, Dvec, &
                      Pmat, Confun, Zwork, Jcntyp, Parprj, Icorct)

        implicit none

        integer, intent(in) :: Ifun
        integer, intent(in) :: Iptb
        integer, intent(in) :: Liwrk
        integer, intent(in) :: Lwrk
        integer, intent(in) :: Nparm
        integer, intent(in) :: Numgr
        integer :: Icorct
        integer :: Indm
        integer :: Ioptn
        integer :: Iphse
        integer :: Mact
        real(wp) :: Projct
        real(wp) :: Rchin
        real(wp) :: Tolcon
        real(wp) :: Unit
        integer :: Iact(Numgr)
        integer :: Icntyp(Numgr)
        integer :: Iwork(Liwrk)
        integer :: Jcntyp(Numgr)
        real(wp) :: Confun(Numgr, Nparm + 1)
        real(wp) :: Dvec(Nparm)
        real(wp) :: Err1(Numgr + 3)
        real(wp) :: Error(Numgr + 3)
        real(wp) :: Fun(Ifun)
        real(wp) :: Parprj(Nparm)
        real(wp) :: Parwrk(Nparm)
        real(wp) :: Pmat(Nparm + 1, Numgr)
        real(wp) :: Pttbl(Iptb, Indm)
        real(wp) :: Work(Lwrk)
        real(wp) :: Zwork(Nparm)

        class(conmax_solver), intent(inout) :: me
        real(wp) :: emin, eold, f1, p1, procor, rchdwn, s, wdist
        integer :: i, ilc06, ilc16, ilc22, ilc24, ilc30, ilc31, ilc33, ilc35, &
                   ilc41, ioptth, ipmax, ipt, ismax, isrcr, j, jflag, &
                   k, l, ncor, newtit, nmaj, nmin, npar1

        real(wp), parameter :: gain = one/(ten*ten)
        integer, parameter :: newtlm = 3 !! Set the limit newtlm on the number of quasi-newton steps (i.e. calls
                                         !! to searcr), and if newtlm > 1 set the parameter gain such that no
                                         !! further newton steps will be tried unless the last step reduced the
                                         !! maximum standard error by a factor of gain or better.

        ! SET MACHINE AND PRECISION DEPENDENT CONSTANTS.
        ilc06 = iloc(6, Nparm, Numgr)
        ilc16 = iloc(16, Nparm, Numgr)
        ilc22 = iloc(22, Nparm, Numgr)
        ilc24 = iloc(24, Nparm, Numgr)
        ilc30 = iloc(30, Nparm, Numgr)
        ilc31 = iloc(31, Nparm, Numgr)
        ilc33 = iloc(33, Nparm, Numgr)
        ilc35 = iloc(35, Nparm, Numgr)
        ilc41 = iloc(41, Nparm, Numgr)
        ioptth = (Ioptn - (Ioptn/100000)*100000)/10000
        npar1 = Nparm + 1
        newtit = 0

        ! FOR NOW, SET JCNTYP(I)=0 IF ICNTYP(I) > 0 AND SET JCNTYP(I)
        ! =ICNTYP(I) OTHERWISE TO DIRECT ERCMP1 TO COMPUTE THE ERRORS FOR THE
        ! STANDARD CONSTRAINTS ONLY.
        do i = 1, Numgr
            if (Icntyp(i) <= 0) then
                Jcntyp(i) = Icntyp(i)
            else
                Jcntyp(i) = 0
            end if
        end do
        ! PUT PARPRJ IN PARWRK FOR USE IN ERCMP1 AND FNSET.
        do j = 1, Nparm
            Parwrk(j) = Parprj(j)
        end do
        ! CALL ERCMP1 WITH ICNUSE=1 TO COMPUTE THE STANDARD ERRORS.
        ! 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, Confun, Jcntyp, ipmax, ismax, Err1)

        ! IF THE TYPE -2 AND TYPE -1 ERROR NORMS ARE BOTH <= TOLCON
        ! WE RETURN WITH ICORCT=-1.
        ! NOTE THAT IN THEORY THE TYPE -1 CONSTRAINTS SHOULD BE NO PROBLEM,
        ! BUT OCCASIONALLY THEY ARE VIOLATED DUE TO ROUNDOFF ERROR OR
        ! PROBLEMS IN WOLFE, SO WE CHECK THEM TO BE SAFE.
        if (Err1(Numgr + 3) > Tolcon) then
            ! HERE THE TYPE -2 ERROR NORM IS > TOLCON AND WE CALL RCHMOD WITH
            ! IRCH=-1 TO SEE IF RCHIN SHOULD BE INCREASED.
            call rchmod(Numgr, Error, Err1, Icntyp, Mact, Iact, ipmax, ismax, Unit, &
                        -1, rchdwn, Rchin)
        else if (Err1(Numgr + 2) <= Tolcon) then
            Icorct = -1
            return
        end if

        ! PUT PARPRJ INTO THE WORK VECTOR ZWORK SO PARPRJ ITSELF WILL REMAIN
        ! UNCHANGED UNLESS CORRCT IS SUCCESSFUL IN CORRECTING BACK INTO THE
        ! FEASIBLE REGION.
        do j = 1, Nparm
            Zwork(j) = Parprj(j)
        end do
        ! COMPUTE EOLD = MAX(ERR1(NUMGR+2),ERR1(NUMGR+3)).  NOTE THAT EOLD IS
        ! POSITIVE SINCE OTHERWISE WE WOULD HAVE RETURNED ABOVE (ASSUMING
        ! TOLCON >= 0.0).  THUS IF ONLY ONE TYPE OF STANDARD CONSTRAINT IS
        ! PRESENT, THE FACT THAT ERR1(NUMGR+2) OR ERR1(NUMGR+3) IS ZERO WILL
        ! DO NO HARM.
        eold = Err1(Numgr + 3)
        if (Err1(Numgr + 2) > eold) eold = Err1(Numgr + 2)

        ! STATEMENTS ABOVE THIS POINT WILL NOT BE EXECUTED AGAIN IN THIS CALL
        ! TO CORRCT.

        main_loop: do

            ! NOW WE SET UP PMAT FOR USE IN WOLFE TO TRY TO COMPUTE A VECTOR DVEC
            ! POINTING BACK INTO THE FEASIBLE REGION.
            ! IF IOPTTH=1 WE CALL DERST ONCE TO PUT THE STANDARD
            ! GRADIENTS IN CONFUN.
            if (ioptth > 0) then
                ! WE SET IPT=-1 TO TELL DERST TO COMPUTE STANDARD CONSTRAINTS ONLY.
                ipt = -1
                call me%derst(Ioptn, Nparm, Numgr, Pttbl, Iptb, Indm, Parwrk, ipt, &
                              Work(ilc24), Work(ilc35), Iwork(ilc22), Confun)
            end if

            l = 0
            do i = 1, Numgr
                if (Icntyp(i) + 1 < 0) then
                    ! HERE ICNTYP(I)=-2 AND WE WILL INCLUDE CONSTRAINT I IF AND ONLY IF
                    ! ERR1(I) >= -RCHIN*PROJCT.  WHEN ICNTYP(I)=-1 WE HAVE A LINEAR
                    ! STANDARD CONSTRAINT AND IT WILL ALWAYS BE INCLUDED.
                    if (Err1(i) + Rchin*Projct < 0) cycle
                else if (Icntyp(i) + 1 /= 0) then
                    cycle
                end if

                if (ioptth <= 0) then
                    ! HERE IOPTTH=0 AND WE HAVE NOT YET PLACED THE GRADIENT OF THE LEFT
                    ! SIDE OF CONSTRAINT I IN CONFUN(I,.) SO WE DO IT NOW.
                    ipt = i
                    call me%derst(Ioptn, Nparm, Numgr, Pttbl, Iptb, Indm, Parwrk, ipt, &
                                  Work(ilc24), Work(ilc35), Iwork(ilc22), Confun)
                end if

                l = l + 1
                ! PUT THE GRADIENT OF THE LEFT SIDE OF CONSTRAINT I IN PMAT(1,L),...,
                ! PMAT(NPARM,L).
                do k = 1, Nparm
                    Pmat(k, l) = Confun(i, k + 1)
                end do

                ! SET ROW NPARM+1 OF PMAT.  WE WILL USUALLY SET PMAT(NPARM+1,L)=
                ! ERR1(I), SO THE WOLFE CONSTRAINT WILL BE GRADIENT(I).DVEC + ERR1(I)
                ! <= 0.0, I.E. (-GRADIENT(I)).DVEC >= ERR1(I).  THE EXCEPTION
                ! OCCURS WHEN ICNTYP(I)=-1 AND ERR1(I) < 0.0, IN WHICH CASE WE
                ! REPLACE ERR1(I) BY ERR1(I)/2.0, IN ORDER TO INSURE THAT EVEN IF PROCOR
                ! TAKES ON ITS MAXIMUM ALLOWED VALUE OF 2.0, NO LINEAR STANDARD
                ! CONSTRAINT WITH NEGATIVE VALUE WILL BECOME POSITIVE VALUED (IGNORING
                ! ROUNDOFF ERROR).  NOTE THAT IF WE DENOTE CONSTRAINT I BY G(I)  <=
                ! 0.0, THEN OUR INEQUALITIES BECOME (GRAD G)(I).DVEC <= -G(I) (OR
                ! -G(I)/2.0), SO A SOLUTION DVEC IS A SOLUTION OF (GRAD G)(I).DVEC =
                ! -G(I) - EPS(I) WHERE EPS(I) = -(GRAD G)(I).DVEC - G(I) = -(GRAD G)(I).
                ! DVEC - G(I)/2.0 - G(I)/2.0 >= 0.0.  NOW WITH H(I) = G(I) + EPS(I)
                ! WE HAVE (GRAD H)(I).DVEC = -H(I), SO IF THIS SYSTEM IS SQUARE THEN
                ! PROCOR=1.0 GIVES A NEWTON STEP FOR SOLVING H(I)=0.0, I.E. G(I) =
                ! -EPS(I) <= 0.0.  THUS WE HAVE A KIND OF GENERALIZED NEWTON METHOD.
                if (Icntyp(i) + 1 == 0) then
                    if (Err1(i) < 0) then
                        Pmat(npar1, l) = Err1(i)/two
                        cycle
                    end if
                end if
                Pmat(npar1, l) = Err1(i)
            end do

            ! CALL WOLFE WITH ISTRT=0 TO COMPUTE DVEC FROM SCRATCH.
            call wolfe(Nparm, l, Pmat, 0, s, ncor, Iwork(ilc16), Iwork, Liwrk, Work, &
                       Lwrk, Work(ilc33), Work(ilc06), Work(ilc31), Work(ilc30), &
                       Nparm, Numgr, Work(ilc41), Dvec, wdist, nmaj, nmin, jflag)
            if (jflag <= 0) then

                ! IN SEARCR AND MULLER WE WILL COMPUTE THE ERROR NORM FOR TYPE -2 AND
                ! TYPE -1 CONSTRAINTS, SO WE LUMP THESE TOGETHER BY SETTING
                ! JCNTYP(I)=-2 IF IT WAS -1.
                do i = 1, Numgr
                    if (Jcntyp(i) + 1 == 0) Jcntyp(i) = -2
                end do
                ! CALL SEARCR TO TRY TO FIND PROCOR SO THAT WITH PARAMETER VECTOR
                ! (OLD) ZWORK + PROCOR*DVEC WE WILL HAVE EMIN = MAX(ERR1(NUMGR+2),
                ! ERR1(NUMGR+3)) <= TOLCON.  IF SEARCR SUCCEEDS IT WILL RETURN WITH
                ! ISRCR=0, WHILE IF IT FAILS IT WILL RETURN WITH ISRCR=1.  IN BOTH
                ! CASES ZWORK WILL BE THE SAME AS BEFORE THE CALL TO SEARCR.
                call me%searcr(Ioptn, Nparm, Numgr, Dvec, Fun, Ifun, Pttbl, Iptb, Indm, &
                               Zwork, Tolcon, Iphse, Iwork, Liwrk, Work, Lwrk, Parwrk, &
                               Err1, p1, f1, procor, emin, isrcr)
                if (isrcr <= 0) then
                    ! HERE THE MAXIMUM STANDARD CONSTRAINT ERROR IS SMALLER
                    ! THAN -TOLCON.  SINCE OVERCORRECTION MAY ADVERSELY AFFECT CONVERGENCE,
                    ! WE CALL MULLER TO TRY TO GET THE MAXIMUM STANDARD CONSTRAINT
                    ! ERROR INTO THE CLOSED INTERVAL (-TOLCON, TOLCON) BY FURTHER
                    ! MODIFYING PROCOR.
                    if (emin + Tolcon < 0) &
                        call me%muller(Ioptn, Nparm, Numgr, Dvec, Fun, &
                                       Ifun, Pttbl, Iptb, Indm, Zwork, Tolcon, Iphse, Iwork, Liwrk, &
                                       Work, Lwrk, Parwrk, Err1, p1, f1, procor, emin)

                    ! NOW COMPUTE PARPRJ = ZWORK + PROCOR*DVEC, SET ICORCT=0, AND RETURN.
                    do j = 1, Nparm
                        Parprj(j) = Zwork(j) + procor*Dvec(j)
                    end do
                    Icorct = 0
                    return
                else
                    newtit = newtit + 1
                    if (newtit < newtlm) then
                        if (emin <= gain*eold) then
                            ! HERE WE UPDATE ZWORK, EOLD, PARWRK, AND ERR1, AND TRY ANOTHER NEWTON
                            ! STEP WITH SEARCR.
                            eold = emin
                            do j = 1, Nparm
                                Zwork(j) = Zwork(j) + procor*Dvec(j)
                                Parwrk(j) = Zwork(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, Confun, &
                                           Jcntyp, ipmax, ismax, Err1)
                            cycle main_loop
                        end if
                    end if
                end if
            end if

            exit ! done
        end do main_loop

        ! HERE WE WERE UNABLE TO OBTAIN A FEASIBLE PARPRJ AND WE RETURN WITH
        ! THE WARNING ICORCT=1.
        Icorct = 1

    end subroutine corrct