rkpar Subroutine

    subroutine rkpar(me, Ioptn, Numgr, Nparm, Icntyp, Mactrk, Iact, Actdif, &
                     Projct, Param, Fun, Ifun, Pttbl, Iptb, Indm, Vder, Pmat, &
                     Ncor, s, Itypm1, Itypm2, Unit, Tolcon, Rchin, Nstep, &
                     Error, Iphse, Iwork, Liwrk, Work, Lwrk, Confun, Vdern, &
                     Vders, Wvec, Parprj, Ifrkpr)

        implicit none

        class(conmax_solver), intent(inout) :: me

        integer, intent(in)  :: Ifun
        integer, intent(in)  :: Indm
        integer, intent(in)  :: Iptb
        integer, intent(in)  :: Liwrk
        integer, intent(in) :: Lwrk
        integer, intent(in)  :: Nparm
        integer, intent(in)  :: Numgr
        integer  :: Ifrkpr
        integer  :: Ioptn
        integer  :: Iphse
        integer  :: Itypm1
        integer  :: Itypm2
        integer  :: Mactrk
        integer  :: Ncor
        integer  :: Nstep
        real(wp) :: Projct
        real(wp) :: Rchin
        real(wp) :: s
        real(wp) :: Tolcon
        real(wp) :: Unit
        real(wp) :: Actdif(Numgr)
        real(wp) :: Confun(Numgr, Nparm + 1)
        real(wp) :: Error(Numgr + 3)
        real(wp) :: Fun(Ifun)
        integer  :: Iact(Numgr)
        integer  :: Icntyp(Numgr)
        integer  :: Iwork(Liwrk)
        real(wp) :: Param(Nparm)
        real(wp) :: Parprj(Nparm)
        real(wp) :: Pmat(Nparm + 1, Numgr)
        real(wp) :: Pttbl(Iptb, Indm)
        real(wp) :: Vder(Nparm)
        real(wp) :: Vdern(Nparm)
        real(wp) :: Vders(Nparm)
        real(wp) :: Work(Lwrk)
        real(wp) :: Wvec(Nparm)

        real(wp) :: p6, wdist, proj1
        integer :: icorct, ilc06, ilc10, ilc11, ilc15, ilc21, ilc27, &
                   ilc30, ilc31, ilc33, ilc40, ilc48, j, jflag, nmaj, &
                   nmin, npar1, nstcnt

        ! SET MACHINE AND PRECISION DEPENDENT CONSTANTS FOR RKPAR.
        ilc06 = iloc(6, Nparm, Numgr)
        ilc10 = iloc(10, Nparm, Numgr)
        ilc11 = iloc(11, Nparm, Numgr)
        ilc15 = iloc(15, Nparm, Numgr)
        ilc21 = iloc(21, Nparm, Numgr)
        ilc27 = iloc(27, Nparm, Numgr)
        ilc30 = iloc(30, Nparm, Numgr)
        ilc31 = iloc(31, Nparm, Numgr)
        ilc33 = iloc(33, Nparm, Numgr)
        ilc40 = iloc(40, Nparm, Numgr)
        ilc48 = iloc(48, Nparm, Numgr)

        ! IFRKPR=0 IS A SIGNAL THAT THE SUBROUTINE OPERATED NORMALLY.
        Ifrkpr = 0
        proj1 = Projct/Nstep
        p6 = proj1/(two + two + two)
        npar1 = Nparm + 1
        nstcnt = 1
        ! PARPRJ WILL BE USED AS THE BASE POINT FOR THE NEXT RK STEP DURING THE
        ! OPERATION OF THIS SUBROUTINE.
        do j = 1, Nparm
            Parprj(j) = Param(j)
            Vdern(j) = Vder(j)
        end do

        main: block

            do
                ! NOTE THAT HERE H*VDERN IS THE K1 OF THE USUAL RUNGE-KUTTA FORMULAE.
                ! SET THE WORK VECTOR WVEC = PARPRJ-PROJ1*VDERN/2.0, THEN CALL PMTST
                ! AND WOLFE TO GET THE VECTOR (AGAIN CALLED VDERN) OF DERIVATIVE VALUES.
                ! THEN H*VDERN WILL BE THE K2 OF THE USUAL RUNGE-KUTTA FORMULAE.
                ! WE WILL ACCUMULATE K1/H + 2.0*K2/H + 2.0*K3/H IN VDERS, AND ADD IN
                ! K4/H AT THE END.
                do j = 1, Nparm
                    Vders(j) = Vdern(j)
                    Wvec(j) = Parprj(j) - proj1*Vdern(j)/two
                end do
                ! IF THERE ARE ANY STANDARD CONSTRAINTS, WE CORRECT BACK INTO THE
                ! FEASIBLE REGION IF POSSIBLE BEFORE CALLING PMTST.
                if (Itypm1 + Itypm2 > 0) then
                    call me%corrct(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Icntyp, &
                                   Unit, Tolcon, Rchin, Error, Mactrk, Iact, Projct, Iphse, &
                                   Iwork, Liwrk, Work, Lwrk, Work(ilc27), Work(ilc11), &
                                   Work(ilc10), Pmat, Confun, Work(ilc48), Iwork(ilc21), &
                                   Wvec, icorct)
                    if (icorct > 0) exit main ! failure
                end if
                call me%pmtst(Ioptn, Numgr, Nparm, Wvec, Icntyp, Mactrk, Iact, Pttbl, Iptb, &
                              Indm, Actdif, Iphse, Iwork, Liwrk, Work, Lwrk, Confun, Pmat)
                call wolfe(Nparm, Mactrk, Pmat, 1, s, Ncor, Iwork(ilc15), Iwork, Liwrk, &
                           Work, Lwrk, Work(ilc33), Work(ilc06), Work(ilc31), &
                           Work(ilc30), Nparm, Numgr, Work(ilc40), Vdern, wdist, nmaj, &
                           nmin, jflag)
                ! IF WOLFE FAILED, SO WILL THIS SUBROUTINE.
                if (jflag <= 0) then
                    ! NOW VDERN REPRESENTS K2/H.  SET WVEC = PARPRJ-PROJ1*VDERN/2.0 AND
                    ! COMPUTE THE NEW VDERN, WHICH WILL REPRESENT K3/H.
                    do j = 1, Nparm
                        Vders(j) = Vders(j) + two*Vdern(j)
                        Wvec(j) = Parprj(j) - proj1*Vdern(j)/two
                    end do
                    ! IF THERE ARE ANY STANDARD CONSTRAINTS, WE CORRECT BACK INTO THE
                    ! FEASIBLE REGION IF POSSIBLE BEFORE CALLING PMTST.
                    if (Itypm1 + Itypm2 > 0) exit main ! failure
                    call me%corrct(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Icntyp, &
                                   Unit, Tolcon, Rchin, Error, Mactrk, Iact, Projct, Iphse, &
                                   Iwork, Liwrk, Work, Lwrk, Work(ilc27), Work(ilc11), &
                                   Work(ilc10), Pmat, Confun, Work(ilc48), Iwork(ilc21), &
                                   Wvec, icorct)
                    if (icorct > 0) exit main ! failure
                end if

                call me%pmtst(Ioptn, Numgr, Nparm, Wvec, Icntyp, Mactrk, Iact, Pttbl, Iptb, &
                              Indm, Actdif, Iphse, Iwork, Liwrk, Work, Lwrk, Confun, Pmat)
                call wolfe(Nparm, Mactrk, Pmat, 1, s, Ncor, Iwork(ilc15), Iwork, Liwrk, &
                           Work, Lwrk, Work(ilc33), Work(ilc06), Work(ilc31), &
                           Work(ilc30), Nparm, Numgr, Work(ilc40), Vdern, wdist, nmaj, &
                           nmin, jflag)
                if (jflag > 0) exit main ! failure

                ! NOW VDERN REPRESENTS K3/H.  SET WVEC = PARPRJ-PROJ1*VDERN AND
                ! COMPUTE THE NEW VDERN, WHICH WILL REPRESENT K4/H.
                do j = 1, Nparm
                    Vders(j) = Vders(j) + two*Vdern(j)
                    Wvec(j) = Parprj(j) - proj1*Vdern(j)
                end do
                ! IF THERE ARE ANY STANDARD CONSTRAINTS, WE CORRECT BACK INTO THE
                ! FEASIBLE REGION IF POSSIBLE BEFORE CALLING PMTST.
                if (Itypm1 + Itypm2 > 0) then
                    call me%corrct(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Icntyp, &
                                   Unit, Tolcon, Rchin, Error, Mactrk, Iact, Projct, Iphse, &
                                   Iwork, Liwrk, Work, Lwrk, Work(ilc27), Work(ilc11), &
                                   Work(ilc10), Pmat, Confun, Work(ilc48), Iwork(ilc21), &
                                   Wvec, icorct)
                    if (icorct > 0) exit main ! failure
                end if
                call me%pmtst(Ioptn, Numgr, Nparm, Wvec, Icntyp, Mactrk, Iact, Pttbl, Iptb, &
                              Indm, Actdif, Iphse, Iwork, Liwrk, Work, Lwrk, Confun, Pmat)
                call wolfe(Nparm, Mactrk, Pmat, 1, s, Ncor, Iwork(ilc15), Iwork, Liwrk, &
                           Work, Lwrk, Work(ilc33), Work(ilc06), Work(ilc31), &
                           Work(ilc30), Nparm, Numgr, Work(ilc40), Vdern, wdist, nmaj, &
                           nmin, jflag)
                if (jflag > 0) exit main ! failure

                ! NOW VDERN REPRESENTS K4/H, SO VDERS + VDERN WILL REPRESENT (K1 +
                ! 2.0*K2 + 2.0*K3 + K4)/H.  PUT THE NEW PARAMETER VECTOR IN PARPRJ.
                do j = 1, Nparm
                    Parprj(j) = Parprj(j) - p6*(Vders(j) + Vdern(j))
                end do
                if (nstcnt < Nstep) then
                    ! HERE NSTCNT < NSTEP AND WE SET UP FOR THE NEXT RK STEP.
                    ! AFTER WE HAVE DONE THIS STEP, VDERN WILL REPRESENT THE VDER1 FOR THE
                    ! NEXT STEP.  PARPRJ ALREADY IS THE BASE POINT FOR THE NEXT STEP.
                    nstcnt = nstcnt + 1
                    ! IF THERE ARE ANY STANDARD CONSTRAINTS, WE CORRECT BACK INTO THE
                    ! FEASIBLE REGION IF POSSIBLE BEFORE CALLING PMTST.
                    if (Itypm1 + Itypm2 > 0) then
                        call me%corrct(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, &
                                       Icntyp, Unit, Tolcon, Rchin, Error, Mactrk, Iact, &
                                       Projct, Iphse, Iwork, Liwrk, Work, Lwrk, Work(ilc27), &
                                       Work(ilc11), Work(ilc10), Pmat, Confun, Work(ilc48), &
                                       Iwork(ilc21), Parprj, icorct)
                        if (icorct > 0) exit main ! failure
                    end if
                    call me%pmtst(Ioptn, Numgr, Nparm, Parprj, Icntyp, Mactrk, Iact, Pttbl, &
                                  Iptb, Indm, Actdif, Iphse, Iwork, Liwrk, Work, Lwrk, Confun, &
                                  Pmat)
                    call wolfe(Nparm, Mactrk, Pmat, 1, s, Ncor, Iwork(ilc15), Iwork, Liwrk, &
                               Work, Lwrk, Work(ilc33), Work(ilc06), Work(ilc31), &
                               Work(ilc30), Nparm, Numgr, Work(ilc40), Vdern, wdist, &
                               nmaj, nmin, jflag)
                    if (jflag > 0) exit main ! failure
                else
                    exit ! done
                end if

            end do

            return ! success
        end block main

        ! failure
        Ifrkpr = 1
        ! WRITE(NWRIT,'(A)') '*****RKPAR HAS FAILED'

    end subroutine rkpar