step_r1 Subroutine

subroutine step_r1(comm,f,tnow,y,yp,stages,tol,htry,y_new,    &
                     errest,err,hmin,phase_2)
!
! Part of rksuite_90 v1.0 (Aug 1994)
!         software for initial value problems in ODEs
!
! Authors: R.W. Brankin (NAG Ltd., Oxford, England)
!          I. Gladwell  (Math Dept., SMU, Dallas, TX, USA)
!          see main doc for contact details
!
type(rk_comm_real_1d), intent(inout), target :: comm
real(kind=wp), intent(out) :: err
real(kind=wp), intent(inout) :: htry                                 !indep!
real(kind=wp), intent(in) :: tnow                                    !indep!
real(kind=wp), intent(in) :: tol
real(kind=wp), intent(in), optional :: hmin                          !indep!
logical, intent(inout), optional :: phase_2
!
real(kind=wp), dimension(:), intent(in) :: y, yp                     !dep!
real(kind=wp), dimension(:), intent(out) ::  errest, y_new           !dep!
real(kind=wp), dimension(:,:), intent(out) :: stages                 !dep!
!
interface
   function f(t,y)
      use rksuite_90_prec, only:wp
      real(kind=wp), intent(in) :: t                                 !indep!
      real(kind=wp), dimension(:), intent(in) :: y                   !dep!
      real(kind=wp), dimension(size(y,1)) :: f                       !dep!
   end function f
end interface
!
real(kind=wp) :: tstg                                                !indep!
integer :: i, j
logical :: cutbak, main
!
intrinsic         abs, max, sign
!
real(kind=wp), dimension(:), pointer :: weights, thresh              !shp-dep!
real(kind=wp), dimension(:,:), pointer :: a                          !real!
real(kind=wp), dimension(:), pointer :: b, bhat, c                   !real!
integer, dimension(:), pointer :: ptr                                !integer!
!
real(kind=wp), parameter :: zero=0.0_wp, half=0.5_wp, one=1.0_wp
!
!  ERREST is used for working storage in this computation.
!
weights => comm%weights
thresh => comm%thresh
a => comm%a
b => comm%b
bhat => comm%bhat
c => comm%c
ptr => comm%ptr
!
main = present(hmin) .and. present(phase_2)
attempt_step: do
!
   if (main) then
      if (comm%phase2) weights = max(thresh,abs(y))
   end if
!
   do i = 2, comm%no_of_stages
      errest = a(i,1)*yp
      do j = 2, i - 1
         if (a(i,j)/=zero) errest = errest + a(i,j)*stages(:,ptr(j))
      end do
      y_new = y + htry*errest
!
!  METHOD = 'M' is special in that an estimate of the local error can be
!  formed before the step is completed.  If the step is a failure,
!  return immediately.  Otherwise, complete the step and compute a more
!  accurate error estimate.
!
      if (comm%rk_method==2 .and. i==7) then
         call stepb
         if (err>tol) return
      end if
!
      tstg = tnow + c(i)*htry
      if (main .and. comm%at_t_end .and. c(i)==one) tstg = comm%t_end
      stages(:,ptr(i)) = f(tstg,y_new)
!
!  Increment the counter for the number of function evaluations
!  depending on whether the primary or secondary integration is taking
!  place.
!
      if (main) then
         comm%f_count = comm%f_count + 1
      else
         comm%ge_f_count = comm%ge_f_count + 1
      end if
!
!  When PHASE2 is .TRUE. we are in the second phase of the automatic
!  selection of the initial step size.  The results of the first three
!  stages are monitored in the subroutine STEPA for evidence that H is
!  too large -- instability and/or an unreliable estimate of the error
!  of the step is then possible.  When the subroutine believes H to be
!  too large, it returns CUTBAK = .TRUE. and a suitably reduced H for
!  another try.
!
      if (main) then
         if (phase_2) then
            if (i<=3 .and. abs(htry)>hmin) then
               call stepa(stages(:,ptr(i)),htry,cutbak)
               if (cutbak) then
                  comm%at_t_end = .false.
!
!  Make sure that STEPA does not reduce the step size below the
!  minimum. If it does, reset H to HMIN and deactivate PHASE2.
!
                  if (abs(htry)<=hmin) then
                     htry = sign(hmin,htry); comm%phase2 = .false.
                  end if
                  cycle attempt_step
               end if
            end if
         end if
      end if
!
   end do
!
!  Some formulas are constructed so that the last stage represents
!  the result of the step (FSAL=.TRUE.), hence if the step is acceptable,
!  it will be the first stage for the next step. When FSAL=.FALSE., we
!  have to complete the computation of the step.
!
   if (.not.comm%fsal) then
      errest = bhat(1)*yp
      do i = 2, comm%no_of_stages
         if (bhat(i)/=zero) errest = errest + bhat(i)*stages(:,ptr(i))
      end do
      y_new = y + htry*errest
   end if
!
!  Form an estimate of the error in the lower order formula by comparing
!  it to the higher order formula of the pair. ERREST has been used
!  as working storage above.  The higher order approximation has been
!  formed as Y_NEW = Y + HTRY*ERREST where ERREST is a linear
!  combination of the stages of the formula. The lower order result also
!  has the form Y plus HTRY times a different linear combination of
!  the stages. Hence, this different linear combination of stages for
!  the lower order formula can just be subtracted from the combination
!  stored in ERREST to produce the errors. The result is then
!  multiplied by HTRY to obtain the error estimate.
!
   if (b(1)/=zero) errest = errest - b(1)*yp
   do i = 2, comm%no_of_stages
      if (b(i)/=zero) errest = errest - b(i)*stages(:,ptr(i))
   end do
   errest = htry*errest
!
!  The error in a solution component is measured relative to a weight
!  that is the larger of a threshold and the size of the solution over
!  the step.  Using the magnitude of a solution component at both ends
!  of the step in the definition of "size" increases the robustness of
!  the test. When global error estimation is specified, the weight
!  vector WEIGHTS is defined by the primary integration and is then
!  used in the secondary integration.
!
   if (main) weights = max(half*(abs(y)+abs(y_new)),thresh)
!
   err = maxval(abs(errest/weights))                                 !spec-ar!
!
   exit attempt_step
!
end do attempt_step
!
contains
!
   subroutine stepa(ypnew,htry,cutbak)
!
   real(kind=wp), intent(inout) :: htry                              !indep!
   real(kind=wp), dimension(:), intent(in) :: ypnew                  !dep!
   logical, intent(out) :: cutbak
!
   real(kind=wp) :: argdif, fdiff, scl, tdiff, twt, ynrm, ystgnm
!
!  Update the weights to account for the current intermediate solution
!  approximation Y_NEW.  Compute the sizes of Y and Y_NEW in the
!  new norm.  The size of the Lipschitz constant is assessed by a difference
!  in the arguments Y, Y_NEW and a difference in the function evaluated
!  at these arguments.
!
   weights = max(weights,abs(y_new))
   ynrm = maxval(abs(y)/weights)                                     !spec-ar!
   ystgnm = maxval(abs(y_new)/weights)                               !spec-ar!
   argdif = maxval(abs(y_new-y)/weights)                             !spec-ar!
   fdiff = maxval(abs(ypnew-yp)/weights)                             !spec-ar!
!
!  The transformation of the equation to autonomous form is done
!  implicitly.  The difference of the arguments must take into account
!  the difference between the values of the independent variable T and
!  TSTG. The difference of the corresponding component of the function
!  is zero because of the way the standard transformation is done.
!
   tdiff = tstg - tnow
   twt = abs(comm%t_end-tnow)
   ynrm = max(ynrm,abs(tnow)/twt)
   ystgnm = max(ystgnm,abs(tstg)/twt)
   argdif = max(argdif,abs(tdiff)/twt)
!
!  The ratio FDIFF/ARGDIF is a lower bound for, and an approximation to,
!  a Lipschitz constant L for the differential equation written in 
!  autonomous form.  First we must ask if the difference ARGDIF is 
!  significant in the precision available.  If it appears to be, we insist
!  that abs(HTRY)*L be less than an approximate radius, STABILITY_RADIUS,
!  of the stability region of the method.  This is more stringent than 
!  necessary for stability, possibly a lot more stringent, but the aim is
!  to get an HTRY small enough that the error estimate for the step is
!  credible.  The reduction is required to be at least as much as the step
!  control parameter RS1. It is necessary to limit the reduction of HTRY
!  at any one time because we may be misled in the size of the reduction
!  that is appropriate due to nonlinearity of the differential equation
!  and to inaccurate weights caused by HTRY much too large.  The reduction
!  is not permitted to be more than the step control parameter RS4.
!
   cutbak = .false.
   if (argdif > comm%round_off*max(ynrm,ystgnm)) then
      if ((abs(htry)*fdiff) > (comm%stability_radius*argdif)) then
         scl = max(comm%rs4,min(comm%rs1, &
                   (comm%stability_radius*argdif)/(abs(htry)*fdiff)))
         htry = scl*htry; cutbak = .true.
      end if
   end if
!
   end subroutine stepa
!
   subroutine stepb 
!
   real(kind=wp), dimension(:), pointer :: e                         !real!
!
   e => comm%e
!
   if (main) then
      err = maxval( abs( e(1)*yp + e(3)*stages(:,ptr(3)) + &         !spec-ar!
                         e(4)*stages(:,ptr(4)) + e(5)*stages(:,ptr(5)) + &
                         e(6)*stages(:,ptr(6)) ) /  &
                    max ( half*(abs(y)+abs(y_new)), thresh ) )
   else
      err = maxval( abs( e(1)*yp + e(3)*stages(:,ptr(3)) + &         !spec-ar!
                         e(4)*stages(:,ptr(4)) + e(5)*stages(:,ptr(5)) + &
                         e(6)*stages(:,ptr(6)) ) / weights ) 
   end if
!
   err = abs(comm%h)*err
!
   end subroutine stepb
!
end subroutine step_r1