ddasik Subroutine

subroutine ddasik( &
   t, y, ydot, neq, icopt, idy, res, jack, psol, h, tscale, &
   wt, jskip, rpar, ipar, savr, delta, r, y0, ydot0, pwk, rwm, iwm, cj, uround, &
   epli, sqrtn, rsqrtn, epscon, ratemax, stptol, jflg, 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. An initial
!! value for \(y\) and initial guess for \(\dot{y}\) are input. The method used is a
!! Newton scheme with Krylov iteration and a linesearch algorithm.

   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(jack_t) :: jack
      !! User-supplied routine to update the preconditioner.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   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) :: jskip
      !! Flag to signal if initial `jack` call is to be skipped.
      !! `1`: skip the call,
      !! `0`: do not skip call.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(out) :: savr(neq)
      !! Real work array.
   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) :: pwk(neq)
      !! Real work array.
   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
      !! Unit roundoff (not used).
   real(rk), intent(in) :: epli ! @note: not the same as 'epslin'!
      !! Convergence test constant for the iterative solution of linear systems. See
      !! [[daskr]] for details.
   real(rk), intent(in) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   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) :: jflg
      !! Flag indicating whether a `jac` routine is supplied by the user.
   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 :: iernew, ierpj, ires, liwp, lrwp, maxnit, maxnj, nj
   real(rk) :: epslin
   logical :: failed

   ! Perform initializations.
   lrwp = iwm(iwloc_lrwp)
   liwp = iwm(iwloc_liwp)
   maxnit = iwm(iwloc_mxnit)
   maxnj = iwm(iwloc_mxnj)
   nj = 0
   epslin = epli*epscon
   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) failed = .true.

   ! Preconditioner update loop
   if (.not. failed) then
      do

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

         ! If a Jacobian routine was supplied, call it.
         if ((jflg .eq. 1) .and. (jskip .eq. 0)) then
            nj = nj + 1
            iwm(iwloc_nje) = iwm(iwloc_nje) + 1
            call jack(res, ires, neq, t, y, ydot, wt, delta, r, h, cj, &
                      rwm(lrwp), iwm(liwp), ierpj, rpar, ipar)

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

         jskip = 0

         ! Call the nonlinear Newton solver.
         call dnsik(t, y, ydot, neq, icopt, idy, res, psol, wt, rpar, ipar, &
                    savr, delta, r, y0, ydot0, pwk, rwm, iwm, cj, tscale, sqrtn, rsqrtn, &
                    epslin, epscon, ratemax, maxnit, stptol, icnflg, icnstr, iernew)

         if ((iernew .eq. 1) .and. (nj .lt. maxnj) .and. (jflg .eq. 1)) then
            ! Up to MXNIT iterations were done, the convergence rate is < 1,
            ! a Jacobian routine is supplied, and the number of JACK calls is less than MXNJ.
            ! Copy the residual SAVR to DELTA, call JACK, and try again.
            delta = savr
            cycle
         else
            exit ! Either success or convergence failure
         end if

      end do
   end if

   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 ddasik