dop853 Subroutine

      subroutine dop853(me,x,y,xend,rtol,atol,iout,idid)

      implicit none

      class(dop853_class),intent(inout)       :: me
      real(wp),intent(inout)                  :: x      !! *input:* initial value of independent variable.
                                                        !! *output:* `x` for which the solution has been computed
                                                        !! (after successful return `x=xend`).
      real(wp),dimension(:),intent(inout)     :: y      !! *input:* initial values for `y`. [size n]
                                                        !!
                                                        !! *output:* numerical solution at `x`.
      real(wp),intent(in)                     :: xend   !! final x-value (`xend`-`x` may be positive or negative)
      real(wp),dimension(:),intent(in)        :: rtol   !! relative error tolerance. `rtol` and `atol`
                                                        !! can be both scalars or else both vectors of length `n`.
      real(wp),dimension(:),intent(in)        :: atol   !! absolute error tolerance. `rtol` and `atol`
                                                        !! can be both scalars or else both vectors of length `n`.
                                                        !! `atol` should be strictly positive (possibly very small)
      integer,intent(in)                      :: iout   !! switch for calling the subroutine `solout`:
                                                        !!
                                                        !! * `iout=0`: subroutine is never called
                                                        !! * `iout=1`: subroutine is called after every successful step
                                                        !! * `iout=2`: dense output is performed after every successful step
                                                        !! * `iout=3`: dense output is performed in steps defined by the user
                                                        !!          (see `xout` above)
      integer,intent(out)                     :: idid   !! reports on successfulness upon return:
                                                        !!
                                                        !! * `idid=1`  computation successful,
                                                        !! * `idid=2`  comput. successful (interrupted by [[solout]]),
                                                        !! * `idid=-1` input is not consistent,
                                                        !! * `idid=-2` larger `nmax` is needed,
                                                        !! * `idid=-3` step size becomes too small.
                                                        !! * `idid=-4` problem is probably stiff (interrupted).

      real(wp) :: beta,fac1,fac2,hmax,safe
      integer  :: i,ieco,iprint,istore,nrdens,nstiff,nmax
      logical  :: arret
      integer  :: itol    !! switch for `rtol` and `atol`:
                          !!
                          !! * `itol=0`: both `rtol` and `atol` are scalars.
                          !!    the code keeps, roughly, the local error of
                          !!    `y(i)` below `rtol*abs(y(i))+atol`.
                          !!
                          !! * `itol=1`: both `rtol` and `atol` are vectors.
                          !!    the code keeps the local error of `y(i)` below
                          !!    `rtol(i)*abs(y(i))+atol(i)`.

      iprint = me%iprint
      me%iout = iout
      arret = .false.

      !check procedures:
      if (.not. associated(me%fcn)) then
          if ( iprint/=0 ) &
                write (iprint,*) &
                'Error in dop853: procedure FCN is not associated.'
          idid = -1
          return
      end if

      if (iout/=0 .and. .not. associated(me%solout)) then
          if ( iprint/=0 ) &
                write (iprint,*) &
                'Error in dop853: procedure SOLOUT must be associated if IOUT/=0.'
          idid = -1
          return
      end if

      !scalar or vector tolerances:
      if (size(rtol)==1 .and. size(atol)==1) then
          itol = 0
      elseif (size(rtol)==me%n .and. size(atol)==me%n) then
          itol = 1
      else
          if ( iprint/=0 ) &
                write (iprint,*) &
                'Error in dop853: improper dimensions for rtol and/or atol.'
          idid = -1
          return
      end if

      ! setting the parameters
      me%nfcn   = 0
      me%nstep  = 0
      me%naccpt = 0
      me%nrejct = 0

      nmax = me%nmax
      nrdens = me%nrdens  !number of dense output components

      ! nstiff parameter for stiffness detection
      if ( me%nstiff<=0 ) then
          nstiff = nmax + 10  !no stiffness check
      else
          nstiff = me%nstiff
      end if

      if ( nrdens<0 .or. me%nrdens>me%n ) then
         if ( iprint/=0 ) write (iprint,*) ' curious input nrdens=' , nrdens
         arret = .true.
      else
         if ( nrdens>0 .and. iout<2 .and. iprint/=0 ) &
                write (iprint,*) ' warning: set iout=2 or iout=3 for dense output '
      end if

      if (size(y)/=me%n) then
          if ( iprint/=0 ) &
            write (iprint,*) ' error: y must have n elements: size(y)= ',size(y)
          arret = .true.
      end if

      safe = me%safe
      fac1 = me%fac1
      fac2 = me%fac2
      beta = me%beta

      if ( me%hmax==0.0_wp ) then
          hmax = xend - x
      else
          hmax = me%hmax
      end if

      me%h = me%hinitial     ! initial step size

      ! when a fail has occurred, we return with idid=-1
      if ( arret ) then

         idid = -1

      else

        ! call to core integrator
        call me%dp86co(x,y,xend,hmax,me%h,rtol,atol,itol,iprint, &
                       iout,idid,nmax,nstiff,safe,beta,fac1,fac2, &
                       me%nfcn,me%nstep,me%naccpt,me%nrejct)

      end if

      end subroutine dop853