ddes Subroutine

   subroutine ddes(me, neq, t, y, tout, info, rtol, atol, idid, &
                   ypout, yp, yy, wt, p, phi, alpha, beta, psi, v, w, sig, g, gi, &
                   h, eps, x, xold, hold, told, delsgn, tstop, twou, fouru, start, &
                   phase1, nornd, stiff, intout, ns, kord, kold, init, ksteps, &
                   kle4, iquit, kprev, ivc, iv, kgi)

      implicit none

      class(ddeabm_class), intent(inout)        :: me
      integer, intent(in)                       :: neq
      real(wp), intent(inout)                   :: t
      real(wp), dimension(neq), intent(inout)    :: y
      real(wp), intent(in)                      :: tout
      integer, dimension(4), intent(inout)       :: info
      real(wp), dimension(:), intent(inout)      :: rtol
      real(wp), dimension(:), intent(inout)      :: atol
      integer, intent(inout)                    :: idid
      real(wp), dimension(neq), intent(inout)    :: ypout
      real(wp), dimension(neq), intent(inout)    :: yp
      real(wp), dimension(neq), intent(inout)    :: yy
      real(wp), dimension(neq), intent(inout)    :: wt
      real(wp), dimension(neq), intent(inout)    :: p
      real(wp), dimension(neq, 16), intent(inout) :: phi
      real(wp), dimension(12), intent(inout)     :: alpha
      real(wp), dimension(12), intent(inout)     :: beta
      real(wp), dimension(12), intent(inout)     :: psi
      real(wp), dimension(12), intent(inout)     :: v
      real(wp), dimension(12), intent(inout)     :: w
      real(wp), dimension(13), intent(inout)     :: sig
      real(wp), dimension(13), intent(inout)     :: g
      real(wp), dimension(11), intent(inout)     :: gi
      real(wp), intent(inout)                   :: h
      real(wp), intent(inout)                   :: eps
      real(wp), intent(inout)                   :: x
      real(wp), intent(inout)                   :: xold
      real(wp), intent(inout)                   :: hold
      real(wp), intent(inout)                   :: told
      real(wp), intent(inout)                   :: delsgn
      real(wp), intent(inout)                   :: tstop
      real(wp), intent(inout)                   :: twou
      real(wp), intent(inout)                   :: fouru
      logical, intent(inout)                    :: start
      logical, intent(inout)                    :: phase1
      logical, intent(inout)                    :: nornd
      logical, intent(inout)                    :: stiff
      logical, intent(inout)                    :: intout
      integer, intent(inout)                    :: ns
      integer, intent(inout)                    :: kord
      integer, intent(inout)                    :: kold
      integer, intent(inout)                    :: init
      integer, intent(inout)                    :: ksteps
      integer, intent(inout)                    :: kle4
      integer, intent(inout)                    :: iquit
      integer, intent(inout)                    :: kprev
      integer, intent(inout)                    :: ivc
      integer, dimension(10), intent(inout)      :: iv
      integer, intent(inout)                    :: kgi

      integer :: k, l, m, ltol, natolp, nrtolp
      real(wp) :: a, absdel, del, dt, ha, u
      logical :: crash
      character(len=8) :: xern1
      character(len=16) :: xern3, xern4

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

      ! on the first call, perform initialization

      if (info(1) == 0) then
         u = d1mach4    ! machine unit roundoff quantity
         twou = 2.0_wp*u   ! set associated machine dependent parameters
         fouru = 4.0_wp*u   !
         iquit = 0          ! set termination flag
         init = 0          ! set initialization indicator
         ksteps = 0          ! set counter for attempted steps
         intout = .false.    ! set indicator for intermediate-output
         stiff = .false.    ! set indicator for stiffness detection
         kle4 = 0          ! set step counter for stiffness detection
         start = .true.     ! set indicators for steps code
         phase1 = .true.     !
         nornd = .true.     !
         info(1) = 1          ! reset info(1) for subsequent calls
      end if

      ! check validity of input parameters on each entry

      if (info(1) /= 0 .and. info(1) /= 1) then
         write (xern1, '(i8)') info(1)
         call report_error('ddes', 'in ddeabm, info(1) must be '// &
                           'set to 0 for the start of a new problem, and must be '// &
                           'set to 1 following an interrupted task.  you are '// &
                           'attempting to continue the integration illegally by '// &
                           'calling the code with info(1) = '//xern1, 3, 1)
         idid = -33
      end if

      if (info(2) /= 0 .and. info(2) /= 1) then
         write (xern1, '(i8)') info(2)
         call report_error('ddes', 'in ddeabm, info(2) must be '// &
                           '0 or 1 indicating scalar and vector error tolerances, '// &
                           'respectively.  you have called the code with info(2) = '// &
                           xern1, 4, 1)
         idid = -33
      end if

      if (info(3) /= 0 .and. info(3) /= 1) then
         write (xern1, '(i8)') info(3)
         call report_error('ddes', 'in ddeabm, info(3) must be '// &
                           '0 or 1 indicating the interval or intermediate-output '// &
                           'mode of integration, respectively.  you have called '// &
                           'the code with  info(3) = '//xern1, 5, 1)
         idid = -33
      end if

      if (info(4) /= 0 .and. info(4) /= 1) then
         write (xern1, '(i8)') info(4)
         call report_error('ddes', 'in ddeabm, info(4) must be '// &
                           '0 or 1 indicating whether or not the integration '// &
                           'interval is to be restricted by a point tstop.  you '// &
                           'have called the code with info(4) = '//xern1, 14, 1)
         idid = -33
      end if

      if (neq < 1) then
         write (xern1, '(i8)') neq
         call report_error('ddes', 'in ddeabm,  the number of '// &
                           'equations neq must be a positive integer.  you have '// &
                           'called the code with  neq = '//xern1, 6, 1)
         idid = -33
      end if

      nrtolp = 0
      natolp = 0
      do k = 1, neq
         if (nrtolp == 0 .and. rtol(k) < 0.0_wp) then
            write (xern1, '(i8)') k
            write (xern3, '(1pe15.6)') rtol(k)
            call report_error('ddes', 'in ddeabm, the relative '// &
                              'error tolerances rtol must be non-negative.  you '// &
                              'have called the code with  rtol('//xern1//') = '// &
                              xern3//'.  in the case of vector error tolerances, '// &
                              'no further checking of rtol components is done.', 7, 1)
            idid = -33
            nrtolp = 1
         end if

         if (natolp == 0 .and. atol(k) < 0.0_wp) then
            write (xern1, '(i8)') k
            write (xern3, '(1pe15.6)') atol(k)
            call report_error('ddes', 'in ddeabm, the absolute '// &
                              'error tolerances atol must be non-negative.  you '// &
                              'have called the code with  atol('//xern1//') = '// &
                              xern3//'.  in the case of vector error tolerances, '// &
                              'no further checking of atol components is done.', 8, 1)
            idid = -33
            natolp = 1
         end if

         if (info(2) == 0) exit
         if (natolp > 0 .and. nrtolp > 0) exit
      end do

      if (info(4) == 1) then
         if (sign(1.0_wp, tout - t) /= sign(1.0_wp, tstop - t) &
             .or. abs(tout - t) > abs(tstop - t)) then
            write (xern3, '(1pe15.6)') tout
            write (xern4, '(1pe15.6)') tstop
            call report_error('ddes', 'in ddeabm, you have '// &
                              'called the code with tout = '//xern3//' but '// &
                              'you have also told the code (info(4) = 1) not to '// &
                              'integrate past the point tstop = '//xern4// &
                              ' these instructions conflict.', 14, 1)
            idid = -33
         end if
      end if

      ! check some continuation possibilities

      if (init /= 0) then
         if (t == tout) then
            write (xern3, '(1pe15.6)') t
            call report_error('ddes', 'in ddeabm, you have '// &
                              'called the code with  t = tout = '//xern3// &
                              '  this is not allowed on continuation calls.', 9, 1)
            idid = -33
         end if

         if (t /= told) then
            write (xern3, '(1pe15.6)') told
            write (xern4, '(1pe15.6)') t
            call report_error('ddes', 'in ddeabm, you have '// &
                              'changed the value of t from '//xern3//' to '// &
                              xern4//'  this is not allowed on continuation calls.', &
                              10, 1)
            idid = -33
         end if

         if (init /= 1) then
            if (delsgn*(tout - t) < 0.0_wp) then
               write (xern3, '(1pe15.6)') tout
               call report_error('ddes', 'in ddeabm, by '// &
                                 'calling the code with tout = '//xern3// &
                                 ' you are attempting to change the direction of '// &
                                 'integration.  this is not allowed without '// &
                                 'restarting.', 11, 1)
               idid = -33
            end if
         end if
      end if

      ! invalid input detected

      if (idid == -33) then
         if (iquit /= -33) then
            iquit = -33
            info(1) = -1
         else
            call report_error('ddes', 'in ddeabm, invalid '// &
                              'input was detected on successive entries.  it is '// &
                              'impossible to proceed because you have not '// &
                              'corrected the problem, so execution is being '// &
                              'terminated.', 12, 2)
         end if
         return
      end if

      ! rtol = atol = 0. is allowed as valid input and interpreted as
      ! asking for the most accurate solution possible. in this case,
      ! the relative error tolerance rtol is reset to the smallest value
      ! fouru which is likely to be reasonable for this method and machine

      do k = 1, neq
         if (rtol(k) + atol(k) <= 0.0_wp) then
            rtol(k) = fouru
            idid = -2
         end if
         if (info(2) == 0) exit
      end do

      if (idid == -2) then
         ! rtol=atol=0 on input, so rtol is changed to a small positive value
         info(1) = -1
         return
      end if

      ! branch on status of initialization indicator
      !  init=0 means initial derivatives and nominal step size and direction not yet set
      !  init=1 means nominal step size and direction not yet set
      !  init=2 means no further initialization required

      if (init == 0) then
         ! more initialization --
         ! -- evaluate initial derivatives
         init = 1
         a = t
         call me%df(a, y, yp)
         if (me%error) return ! user-triggered error
         if (t == tout) then
            idid = 2
            ypout = yp
            told = t
            return
         end if
      end if

      if (init == 1) then
         ! -- set independent and dependent variables
         !                      x and yy(*) for steps
         ! -- set sign of integration direction
         ! -- initialize the step size
         init = 2
         x = t
         yy = y
         delsgn = sign(1.0_wp, tout - t)
         h = sign(max(fouru*abs(x), abs(tout - x)), tout - x)
      end if

      ! on each call set information which determines the allowed interval
      ! of integration before returning with an answer at tout

      del = tout - t
      absdel = abs(del)

      do

         ! if already past output point, interpolate and return
         if (abs(x - t) >= absdel) then
            call dintp(x, yy, tout, y, ypout, neq, kold, phi, ivc, iv, kgi, gi, alpha, g, w, xold, p)
            idid = 3
            if (x == tout) then
               idid = 2
               intout = .false.
            end if
            t = tout
            told = t
            return
         end if

         ! if cannot go past tstop and sufficiently close, extrapolate and return
         if (info(4) == 1 .and. abs(tstop - x) < fouru*abs(x)) then
            dt = tout - x
            y = yy + dt*yp
            call me%df(tout, y, ypout)
            if (me%error) return ! user-triggered error
            idid = 3
            t = tout
            told = t
            return
         end if

         ! intermediate-output mode
         if (info(3) /= 0 .and. intout) then
            idid = 1
            y = yy
            ypout = yp
            t = x
            told = t
            intout = .false.
            return
         end if

         ! monitor number of steps attempted

         if (ksteps > me%maxnum) then
            ! a significant amount of work has been expended
            idid = -1
            ksteps = 0
            if (stiff) then
               ! problem appears to be stiff
               idid = -4
               stiff = .false.
               kle4 = 0
            end if
            y = yy
            ypout = yp
            t = x
            told = t
            info(1) = -1
            intout = .false.
            return
         end if

         ! limit step size, set weight vector and take a step

         ha = abs(h)
         if (info(4) == 1) ha = min(ha, abs(tstop - x))
         h = sign(ha, h)
         eps = 1.0_wp
         ltol = 1
         do l = 1, neq
            if (info(2) == 1) ltol = l
            wt(l) = rtol(ltol)*abs(yy(l)) + atol(ltol)
            if (wt(l) <= 0.0_wp) then
               ! relative error criterion inappropriate
               idid = -3
               y = yy
               ypout = yp
               t = x
               told = t
               info(1) = -1
               intout = .false.
               return
            end if
         end do

         call me%dsteps(neq, yy, x, h, eps, wt, start, hold, kord, kold, crash, phi, p, &
                        yp, psi, alpha, beta, sig, v, w, g, phase1, ns, nornd, ksteps, &
                        twou, fouru, xold, kprev, ivc, iv, kgi, gi)
         if (me%error) return !user-triggered error

         if (crash) then
            ! tolerances too small
            idid = -2
            rtol = eps*rtol
            atol = eps*atol
            y = yy
            ypout = yp
            t = x
            told = t
            info(1) = -1
            intout = .false.
            return
         end if

         ! (stiffness test) count number of consecutive steps taken with the
         ! order of the method being less or equal to four

         kle4 = kle4 + 1
         if (kold > 4) kle4 = 0
         stiff = (kle4 >= 50)
         intout = .true.

      end do

   end subroutine ddes