ddasid Subroutine

subroutine ddasid( &
   t, y, ydot, neq, icopt, idy, res, jac, dum1, h, tscale, &
   wt, dum2, rpar, ipar, dum3, delta, r, y0, ydot0, dum4, rwm, iwm, cj, uround, &
   dum5, dum6, dum7, epscon, ratemax, stptol, dum8, icnflg, icnstr, iernls)
!! This routine solves a nonlinear system of algebraic equations of the form:
!!
!!  $$ G(t, y, \dot{y}) = 0 $$
!!
!! for the unknown parts of \(y\) and \(\dot{y}\) in the initial conditions. The method
!! used is a modified Newton scheme.
!!
!! All arguments with "dum" in their names are dummy arguments which are not used in
!! this routine.

   use daskr_kinds, only: rk
   use daskr
   implicit none

   real(rk), intent(in) :: t
      !! Independent variable.
   real(rk), intent(inout) :: y(neq)
      !! Solution vector.
   real(rk), intent(inout) :: ydot(neq)
      !! Derivative of solution vector.
   integer, intent(in) :: neq
      !! Problem size.
   integer, intent(in) :: icopt
      !! Initial condition flag.
   integer, intent(in) :: idy(neq)
      !! Array indicating which variables are differential and which are algebraic.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   procedure(jac_t) :: jac
      !! User-defined Jacobian routine.
   procedure(psol_t) :: dum1
      !! Dummy argument.
   real(rk), intent(in) :: h
      !! Scaling factor for this initial condition calculation.
   real(rk), intent(in) :: tscale ! @todo: what is "t"?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(inout) :: wt(neq)
      !! Weights for error criterion.
   integer, intent(inout) :: dum2
      !! Dummy argument.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(out) :: dum3(neq)
      !! Dummy argument.
   real(rk), intent(inout) :: delta(neq)
      !! Real work array.
   real(rk), intent(out) :: r(neq)
      !! Real work array.
   real(rk), intent(out) :: y0(neq)
      !! Real work array.
   real(rk), intent(out) :: ydot0(neq)
      !! Real work array.
   real(rk), intent(inout) :: dum4(neq)
      !! Dummy argument.
   real(rk), intent(inout) :: rwm(*)
      !! Real workspace for the linear system solver.
   integer, intent(inout) :: iwm(*)
      !! Integer workspace for the linear system solver.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: uround ! @todo: used?
      !! Unit roundoff.
   real(rk), intent(in) :: dum5
      !! Dummy argument.
   real(rk), intent(in) :: dum6
      !! Dummy argument.
   real(rk), intent(in) :: dum7
      !! Dummy argument.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   real(rk), intent(in) :: ratemax
      !! Maximum convergence rate for which Newton iteration is considered converging.
   real(rk), intent(in) :: stptol
      !! Tolerance used in calculating the minimum lambda (`rl`) value allowed.
   integer, intent(in) :: dum8
      !! Dummy argument.
   integer, intent(in) :: icnflg
      !! Constraint flag. If nonzero, then constraint violations in the proposed new
      !! approximate solution will be checked for, and the maximum step length will be
      !! adjusted accordingly.
   integer, intent(out) :: icnstr(neq)
      !! Flags for checking constraints.
   integer, intent(out) :: iernls
      !! Error flag for nonlinear solver.
      !! `0`: nonlinear solver converged.
      !! `1`, `2`: recoverable error inside nonlinear solver.
      !! `1`: retry with current (y,y').
      !! `2`: retry with original (y,y').
      !! `-1`: unrecoverable error in nonlinear solver.

   integer :: ierj, iernew, ires, mxnit, mxnj, nj
   logical :: failed

   ! Perform initializations.
   mxnit = iwm(iwloc_mxnit)
   mxnj = iwm(iwloc_mxnj)
   iernls = 0
   nj = 0
   failed = .false.

   ! Call RES to initialize DELTA.
   ires = 0
   call res(t, y, ydot, cj, delta, ires, rpar, ipar)
   iwm(iwloc_nre) = iwm(iwloc_nre) + 1

   if (ires .lt. 0) then
      failed = .true.
   end if

   ! Looping point for updating the Jacobian.
   if (.not. failed) then
      do

         ! Set all error flags to zero.
         ierj = 0
         ires = 0
         iernew = 0

         ! Reevaluate the iteration matrix, J = dG/dY + CJ*dG/dYDOT.
         call dmatd(neq, t, y, ydot, delta, cj, h, ierj, wt, r, rwm, iwm, res, ires, &
                    uround, jac, rpar, ipar)
         nj = nj + 1
         iwm(iwloc_nje) = iwm(iwloc_nje) + 1

         if ((ires .lt. 0) .or. (ierj .ne. 0)) then
            failed = .true.
            exit
         end if

         ! Call the nonlinear Newton solver.
         call dnsid(t, y, ydot, neq, icopt, idy, res, wt, rpar, ipar, delta, r, &
                    y0, ydot0, rwm, iwm, cj, tscale, epscon, ratemax, mxnit, stptol, &
                    icnflg, icnstr, iernew)

         if ((iernew .eq. 1) .and. (nj .lt. mxnj)) then
            ! MXNIT iterations were done, the convergence rate is < 1,
            ! and the number of Jacobian evaluations is less than MXNJ.
            ! Call RES, reevaluate the Jacobian, and try again.
            ires = 0
            call res(t, y, ydot, cj, delta, ires, rpar, ipar)
            iwm(iwloc_nre) = iwm(iwloc_nre) + 1
            if (ires .lt. 0) then
               failed = .true.
               exit
            end if
            cycle
         else
            exit ! Either success or convergence failure
         end if

      end do
   end if

   ! Unsuccessful exits from nonlinear solver.
   ! Compute IERNLS accordingly.
   if (failed) then
      iernls = 2
      if (ires .le. -2) iernls = -1
      return
   end if

   if (iernew .ne. 0) then
      iernls = min(iernew, 2)
   end if

end subroutine ddasid