dsteps Subroutine

   subroutine dsteps(me, neqn, y, x, h, eps, wt, start, hold, k, &
                     kold, crash, phi, p, yp, psi, alpha, beta, sig, v, w, g, &
                     phase1, ns, nornd, ksteps, twou, fouru, xold, kprev, ivc, iv, &
                     kgi, gi)

      implicit none

      class(ddeabm_class), intent(inout) :: me
      real(wp), intent(inout) :: x      !! independent variable
      real(wp), intent(inout) :: h      !! appropriate step size for next step.  normally determined by code
      real(wp), intent(inout) :: eps    !! local error tolerance
      real(wp), intent(inout) :: hold   !! step size used for last successful step
      integer, intent(inout)  :: k      !! appropriate order for next step (determined by code)
      integer, intent(inout)  :: kold   !! order used for last successful step
      integer, intent(inout)  :: ns     !! number of dsteps taken with size h
      integer, intent(inout)  :: ksteps !! counter on attempted steps
      real(wp), intent(inout) :: twou   !! 2.*u where u is machine unit roundoff quantity
      real(wp), intent(inout) :: fouru  !! 4.*u where u is machine unit roundoff quantity
      real(wp), intent(inout) :: xold   !! required for the interpolation subroutine [[dintp]]
      integer, intent(inout)  :: kprev
      integer, intent(inout)  :: ivc    !! required for the interpolation subroutine [[dintp]]
      integer, intent(inout)  :: kgi    !! required for the interpolation subroutine [[dintp]]
      integer, intent(in)     :: neqn   !! number of equations to be integrated
      logical, intent(inout)  :: start  !! logical variable set .true. for first step, .false. otherwise
      logical, intent(inout)  :: crash  !! logical variable set .true. when no step can be taken, .false. otherwise.
      logical, intent(inout)  :: phase1
      logical, intent(inout)  :: nornd
      real(wp), dimension(neqn), intent(inout) :: y      !! solution vector at x
      real(wp), dimension(neqn), intent(inout) :: wt     !! vector of weights for error criterion
      real(wp), dimension(neqn, 16), intent(inout) :: phi !! required for the interpolation subroutine [[dintp]]
      real(wp), dimension(neqn), intent(inout) :: p      !! required for the interpolation subroutine [[dintp]]
      real(wp), dimension(neqn), intent(inout) :: yp     !! derivative of solution vector at x after successful step
      real(wp), dimension(12), intent(inout) :: psi
      real(wp), dimension(12), intent(inout) :: alpha    !! required for the interpolation subroutine [[dintp]]
      real(wp), dimension(12), intent(inout) :: beta
      real(wp), dimension(13), intent(inout) :: sig
      real(wp), dimension(12), intent(inout) :: v
      real(wp), dimension(12), intent(inout) :: w        !! required for the interpolation subroutine [[dintp]]
      real(wp), dimension(13), intent(inout) :: g        !! required for the interpolation subroutine [[dintp]]
      real(wp), dimension(11), intent(inout) :: gi       !! required for the interpolation subroutine [[dintp]]
      integer, dimension(10), intent(inout)  :: iv      !! required for the interpolation subroutine [[dintp]]

      integer  :: i, ifail, im1, ip1, iq, j, km1, km2, knew, &
                  kp1, kp2, l, limit1, limit2, nsm2, &
                  nsp1, nsp2, jv
      real(wp) :: absh, big, &
                  erk, erkm1, erkm2, erkp1, err, &
                  hnew, p5eps, r, &
                  reali, realns, rho, round, sum, tau, temp1, &
                  temp2, temp3, temp4, temp5, temp6, u

      if (me%error) return !user-triggered error

      !       ***     begin block 0     ***
      !   check if step size or error tolerance is too small for machine
      !   precision.  if first step, initialize phi array and estimate a
      !   starting step size.
      !                   ***

      ! if step size is too small, determine an acceptable one
      crash = .true.
      if (abs(h) < fouru*abs(x)) then
         h = sign(fouru*abs(x), h)
         return
      end if

      p5eps = 0.5_wp*eps

      ! if error tolerance is too small, increase it to an acceptable value
      round = 0.0_wp
      do l = 1, neqn
         round = round + (y(l)/wt(l))**2
      end do
      round = twou*sqrt(round)
      if (p5eps < round) then
         eps = 2.0_wp*round*(1.0_wp + fouru)
         return
      end if

      crash = .false.
      g(1) = 1.0_wp
      g(2) = 0.5_wp
      sig(1) = 1.0_wp

      if (start) then

         ! initialize.  compute appropriate step size for first step
         ! [There are three modes]

         if (me%initial_step_mode == 2) then
            call me%df(x, y, yp)
            if (me%error) return ! user-triggered error
         end if
         u = d1mach4
         big = sqrt(d1mach2)
         if (me%initial_step_mode == 2) sum = 0.0_wp
         do l = 1, neqn
            phi(l, 1) = yp(l)
            phi(l, 2) = 0.0_wp
            if (me%initial_step_mode == 2) sum = sum + (yp(l)/wt(l))**2
         end do
         select case (me%initial_step_mode)
         case (1)
            call me%dhstrt(neqn, x, x + h, y, yp, wt, 1, u, big, &
                           phi(1, 3), phi(1, 4), phi(1, 5), phi(1, 6), h)
            if (me%error) return !user-triggered error
            me%initial_step_size = h ! save it
         case (2)
            sum = sqrt(sum)
            absh = abs(h)
            if (eps < 16.0_wp*sum*h*h) absh = 0.25_wp*sqrt(eps/sum)
            h = sign(max(absh, fouru*abs(x)), h)
            me%initial_step_size = h ! save it
         case (3)
            h = sign(abs(me%initial_step_size), h)   ! user specified
         case default
            error stop 'invalid value for initial_step_mode'
         end select

         hold = 0.0_wp
         k = 1
         kold = 0
         kprev = 0
         start = .false.
         phase1 = .true.
         nornd = .true.
         if (p5eps <= 100.0_wp*round) then
            nornd = .false.
            do l = 1, neqn
               phi(l, 15) = 0.0_wp
            end do
         end if

      end if

      ifail = 0

      !       ***     end block 0     ***

      !       ***     begin block 1     ***
      !   compute coefficients of formulas for this step.  avoid computing
      !   those quantities not changed when step size is not changed.
      !                   ***
      main: do

         kp1 = k + 1
         kp2 = k + 2
         km1 = k - 1
         km2 = k - 2

         ! ns is the number of dsteps taken with size h, including the current
         ! one.  when k<ns, no coefficients change

         if (h /= hold) ns = 0
         if (ns <= kold) ns = ns + 1
         nsp1 = ns + 1

         if (k >= ns) then

            ! compute those components of alpha(*),beta(*),psi(*),sig(*) which
            ! are changed

            beta(ns) = 1.0_wp
            realns = ns
            alpha(ns) = 1.0_wp/realns
            temp1 = h*realns
            sig(nsp1) = 1.0_wp
            if (k >= nsp1) then
               do i = nsp1, k
                  im1 = i - 1
                  temp2 = psi(im1)
                  psi(im1) = temp1
                  beta(i) = beta(im1)*psi(im1)/temp2
                  temp1 = temp2 + h
                  alpha(i) = h/temp1
                  reali = i
                  sig(i + 1) = reali*alpha(i)*sig(i)
               end do
            end if
            psi(k) = temp1

            ! compute coefficients g(*)
            !
            ! initialize v(*) and set w(*).

            if (ns > 1) then

               ! if order was raised, update diagonal part of v(*)

               if (k > kprev) then
                  if (ivc /= 0) then
                     jv = kp1 - iv(ivc)
                     ivc = ivc - 1
                  else
                     jv = 1
                     temp4 = k*kp1
                     v(k) = 1.0_wp/temp4
                     w(k) = v(k)
                     if (k == 2) then
                        kgi = 1
                        gi(1) = w(2)
                     end if
                  end if
                  nsm2 = ns - 2
                  if (nsm2 >= jv) then
                     do j = jv, nsm2
                        i = k - j
                        v(i) = v(i) - alpha(j + 1)*v(i + 1)
                        w(i) = v(i)
                     end do
                     if (i == 2) then
                        kgi = ns - 1
                        gi(kgi) = w(2)
                     end if
                  end if
               end if

               !update v(*) and set w(*)

               limit1 = kp1 - ns
               temp5 = alpha(ns)
               do iq = 1, limit1
                  v(iq) = v(iq) - temp5*v(iq + 1)
                  w(iq) = v(iq)
               end do
               g(nsp1) = w(1)
               if (limit1 /= 1) then
                  kgi = ns
                  gi(kgi) = w(2)
               end if
               w(limit1 + 1) = v(limit1 + 1)
               if (k < kold) then
                  ivc = ivc + 1
                  iv(ivc) = limit1 + 2
               end if

            else

               do iq = 1, k
                  temp3 = iq*(iq + 1)
                  v(iq) = 1.0_wp/temp3
                  w(iq) = v(iq)
               end do
               ivc = 0
               kgi = 0
               if (k /= 1) then
                  kgi = 1
                  gi(1) = w(2)
               end if

            end if

            ! compute the g(*) in the work vector w(*)

            nsp2 = ns + 2
            kprev = k
            if (kp1 >= nsp2) then
               do i = nsp2, kp1
                  limit2 = kp2 - i
                  temp6 = alpha(i - 1)
                  do iq = 1, limit2
                     w(iq) = w(iq) - temp6*w(iq + 1)
                  end do
                  g(i) = w(1)
               end do
            end if

         end if

         !       ***     end block 1     ***

         !       ***     begin block 2     ***
         !   predict a solution p(*), evaluate derivatives using predicted
         !   solution, estimate local error at order k and errors at orders k,
         !   k-1, k-2 as if constant step size were used.
         !                   ***
         !
         !   increment counter on attempted dsteps
         !
         ksteps = ksteps + 1
         !
         !   change phi to phi star
         !
         if (k >= nsp1) then
            do i = nsp1, k
               temp1 = beta(i)
               do l = 1, neqn
                  phi(l, i) = temp1*phi(l, i)
               end do
            end do
         end if
         !
         !   predict solution and differences
         !
         do l = 1, neqn
            phi(l, kp2) = phi(l, kp1)
            phi(l, kp1) = 0.0_wp
            p(l) = 0.0_wp
         end do
         do j = 1, k
            i = kp1 - j
            ip1 = i + 1
            temp2 = g(i)
            do l = 1, neqn
               p(l) = p(l) + temp2*phi(l, i)
               phi(l, i) = phi(l, i) + phi(l, ip1)
            end do
         end do
         if (nornd) then
            p = y + h*p
         else
            do l = 1, neqn
               tau = h*p(l) - phi(l, 15)
               p(l) = y(l) + tau
               phi(l, 16) = (p(l) - y(l)) - tau
            end do
         end if
         xold = x
         x = x + h
         absh = abs(h)
         call me%df(x, p, yp)
         if (me%error) return ! user-triggered error
         !
         !   estimate errors at orders k,k-1,k-2
         !
         erkm2 = 0.0_wp
         erkm1 = 0.0_wp
         erk = 0.0_wp
         do l = 1, neqn
            temp3 = 1.0_wp/wt(l)
            temp4 = yp(l) - phi(l, 1)
            if (km2 >= 0.0_wp) then
               if (km2 > 0.0_wp) erkm2 = erkm2 + ((phi(l, km1) + temp4)*temp3)**2
               erkm1 = erkm1 + ((phi(l, k) + temp4)*temp3)**2
            end if
            erk = erk + (temp4*temp3)**2
         end do
         if (km2 >= 0.0_wp) then
            if (km2 > 0.0_wp) erkm2 = absh*sig(km1)*me%gstr(km2)*sqrt(erkm2)
            erkm1 = absh*sig(k)*me%gstr(km1)*sqrt(erkm1)
         end if
         temp5 = absh*sqrt(erk)
         err = temp5*(g(k) - g(kp1))
         erk = temp5*sig(kp1)*me%gstr(k)
         knew = k

         !   test if order should be lowered
         if (km2 > 0.0_wp) then
            if (max(erkm1, erkm2) <= erk) knew = km1
         else if (km2 == 0.0_wp) then
            if (erkm1 <= 0.5_wp*erk) knew = km1
         end if
         !       ***     end block 2     ***

         !test if step successful
         if (err > eps) then

            !       ***     begin block 3     ***
            !   the step is unsuccessful.  restore  x, phi(*,*), psi(*) .
            !   if third consecutive failure, set order to one.  if step fails more
            !   than three times, consider an optimal step size.  double error
            !   tolerance and return if estimated step size is too small for machine
            !   precision.

            !   restore x, phi(*,*) and psi(*)

            phase1 = .false.
            x = xold
            do i = 1, k
               temp1 = 1.0_wp/beta(i)
               ip1 = i + 1
               do l = 1, neqn
                  phi(l, i) = temp1*(phi(l, i) - phi(l, ip1))
               end do
            end do
            if (k >= 2) then
               do i = 2, k
                  psi(i - 1) = psi(i) - h
               end do
            end if

            !   on third failure, set order to one.  thereafter, use optimal step
            !   size

            ifail = ifail + 1
            temp2 = 0.5_wp
            if (ifail - 3 >= 0.0_wp) then
               if (ifail - 3 > 0.0_wp .and. p5eps < 0.25_wp*erk) &
                  temp2 = sqrt(p5eps/erk)
               knew = 1
            end if
            h = temp2*h
            k = knew
            ns = 0
            if (abs(h) >= fouru*abs(x)) cycle main
            crash = .true.
            h = sign(fouru*abs(x), h)
            eps = eps + eps
            return

         else
            exit main
         end if

      end do main

      !       ***     end block 3     ***

      block4: block
         !       ***     begin block 4     ***
         !   the step is successful.  correct the predicted solution, evaluate
         !   the derivatives using the corrected solution and update the
         !   differences.  determine best order and step size for next step.

         kold = k
         hold = h

         !   correct and evaluate

         temp1 = h*g(kp1)
         if (nornd) then
            do l = 1, neqn
               temp3 = y(l)
               y(l) = p(l) + temp1*(yp(l) - phi(l, 1))
               p(l) = temp3
            end do
         else
            do l = 1, neqn
               temp3 = y(l)
               rho = temp1*(yp(l) - phi(l, 1)) - phi(l, 16)
               y(l) = p(l) + rho
               phi(l, 15) = (y(l) - p(l)) - rho
               p(l) = temp3
            end do
         end if
         call me%df(x, y, yp)
         if (me%error) return ! user-triggered error

         !   update differences for next step

         do l = 1, neqn
            phi(l, kp1) = yp(l) - phi(l, 1)
            phi(l, kp2) = phi(l, kp1) - phi(l, kp2)
         end do
         do i = 1, k
            do l = 1, neqn
               phi(l, i) = phi(l, i) + phi(l, kp1)
            end do
         end do

         !   estimate error at order k+1 unless:
         !     in first phase when always raise order,
         !     already decided to lower order,
         !     step size not constant so estimate unreliable

         erkp1 = 0.0_wp
         if (knew == km1 .or. k == 12) phase1 = .false.
         if (phase1) then
            call raise_order()
            exit block4
         end if
         if (knew == km1) then
            call lower_order()
            exit block4
         end if
         if (kp1 > ns) exit block4
         do l = 1, neqn
            erkp1 = erkp1 + (phi(l, kp2)/wt(l))**2
         end do
         erkp1 = absh*me%gstr(kp1)*sqrt(erkp1)

         !   using estimated error at order k+1, determine appropriate order
         !   for next step

         if (k > 1) then
            if (erkm1 <= min(erk, erkp1)) then
               call lower_order()
               exit block4
            end if
            if (erkp1 >= erk .or. k == 12) exit block4
         else
            if (erkp1 >= 0.5_wp*erk) exit block4
         end if

         !   here erkp1 < erk < max(erkm1,erkm2) else order would have
         !   been lowered in block 2.  thus order is to be raised
         call raise_order()

      end block block4

      ! with new order determine appropriate step size for next step
      hnew = h + h
      if (.not. phase1) then
         if (p5eps < erk*me%two(k + 1)) then
            hnew = h
            if (p5eps < erk) then
               temp2 = k + 1
               r = (p5eps/erk)**(1.0_wp/temp2)
               hnew = absh*max(0.5_wp, min(0.9_wp, r))
               hnew = sign(max(hnew, fouru*abs(x)), h)
            end if
         end if
      end if
      h = hnew
      !       ***     end block 4     ***

   contains

      subroutine raise_order()
         !! raise order
         k = kp1
         erk = erkp1
      end subroutine raise_order

      subroutine lower_order()
         !! lower order
         k = km1
         erk = erkm1
      end subroutine lower_order

   end subroutine dsteps