recursive subroutine step_integrate_r1(comm,f,t_now,y_now,yderiv_now,flag)
!
! 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
!
real(kind=wp), intent(out) :: t_now !indep!
integer, intent(out), optional :: flag
type(rk_comm_real_1d), intent(inout) :: comm
real(kind=wp), dimension(:), intent(out) :: y_now, yderiv_now !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
!
character(len=*), parameter :: srname="STEP_INTEGRATE"
!
real(kind=wp) :: hmin, htry !indep!
real(kind=wp) :: alpha, beta, err, tau, t1, t2, ypnorm, extra_wk
integer :: ier, nrec, state
logical :: failed, phase1, phase3, toomch, sure_stiff
!
integer, parameter :: not_ready=-1, usable=-2, fatal=911, catastrophe=912, &
max_f_count=5000, just_fine=1
logical, parameter :: tell=.false., ask=.true.
real(kind=wp),parameter :: zero=0.0_wp, pt1=0.1_wp, pt9=0.9_wp, one=1.0_wp, &
two=2.0_wp, hundrd=100.0_wp
!
ier = just_fine; nrec = 0
!
! Is it permissible to call STEP_INTEGRATE?
!
body: do
state = get_saved_state_r1("SETUP",comm%save_states)
if (state==fatal) then
ier = catastrophe; nrec = 1; write (comm%rec,"(a)") &
" ** A catastrophic error has already been detected elsewhere."
exit body
end if
if (state==not_ready) then
ier = fatal; nrec = 1; write (comm%rec,"(a)") &
" ** You have not called SETUP, so you cannot use STEP_INTEGRATE."
exit body
end if
if (comm%use_range) then
if (get_saved_state_r1("RANGE_INTEGRATE",comm%save_states)/=usable) then
ier = fatal; nrec = 2; write (comm%rec,"(a/a)") &
" ** You have called STEP_INTEGRATE after you specified in SETUP that you", &
" ** were going to use RANGE_INTEGRATE. This is not permitted."
comm%use_range = .false.
exit body
end if
end if
state = get_saved_state_r1(srname,comm%save_states)
if (state==5 .or. state==6) then
ier = fatal; nrec = 3; write (comm%rec,"(a/a/a)") &
" ** STEP_INTEGRATE has already returned with a flag value of 5 or 6. You",&
" ** cannot continue integrating this problem. You must call SETUP to start ",&
" ** another problem."
exit body
end if
!
if (comm%at_t_start) then
!
comm%yp = f(comm%t,comm%y); comm%f_count = comm%f_count + 1
if (comm%erason) comm%ge_yp = comm%yp
!
! The weights for the control of the error depend on the size of the
! solution at the beginning and at the end of the step. On the first
! step we do not have all this information. Whilst determining the
! initial step size we initialize each component of WEIGHTS to the
! larger of the corresponding component of both abs(Y) and the threshold.
!
comm%weights = max(abs(comm%y),comm%thresh)
!
! If H_START is equal to zero, the code is to find an on-scale initial
! step size H. STEP_INTEGRATE has an elaborate scheme of three phases for
! finding such an H, and some preparations are made earlier. In SETUP an
! upper bound is placed on H that reflects the scale of the independent
! variable. RANGE_INTEGRATE, when used, refines this bound using the
! first output point. Here in STEP_INTEGRATE PHASE1 applies a rule of
! thumb based on the error control, the order of the the formula, and the
! size of the initial slope to get a crude approximation to an on-scale H.
! PHASE2 may reduce H in the course of taking the first step. PHASE3
! repeatedly adjusts H and retakes the first step until H is on scale.
!
! A guess for the magnitude of the first step size H can be provided to SETUP
! as H_START. If it is too big or too small, it is ignored and the automatic
! determination of an on-scale initial step size is activated. If it is
! acceptable, H is set to H_START in SETUP. Even when H is supplied to
! STEP_INTEGRATE, PHASE3 of the scheme for finding an on-scale initial step
! size is made active so that the code can deal with a bad guess.
!
phase1 = comm%h_start == zero; comm%phase2 = phase1; phase3 = .true.
if (phase1) then
comm%h = abs(comm%h)
ypnorm = max(zero, &
maxval(abs(comm%yp)/comm%weights,mask=comm%y/=zero)) !spec-ar1!
tau = comm%tol**comm%expon
if (comm%h*ypnorm > tau) comm%h = tau/ypnorm
hmin = max(comm%sqtiny,comm%toosml* &
max(abs(comm%t_start),abs(comm%t_end)))
comm%h = comm%dir*max(comm%h,hmin)
phase1 = .false.
end if
!
else
!
! Continuation call
!
if (comm%at_t_end) then
ier = fatal; nrec = 3; write (comm%rec,"(a,e13.5,a/a/a)") &
" ** You have already reached T_END ( = ",comm%t_end, "). To integrate",&
" ** furhter with the same problem you must call the routine RESET_T_END",&
" ** with a new value of T_END."
exit body
end if
end if
!
! Begin computation of a step here.
!
failed = .false.
!
take_step: do
!
comm%h = sign(abs(comm%h),comm%dir)
!
! Reduce the step size if necessary so that the code will not step
! past T_END. "Look ahead" to prevent unnecessarily small step sizes.
!
comm%at_t_end = comm%dir*((comm%t+comm%h)-comm%t_end) >= zero
if (comm%at_t_end) then
comm%h = comm%t_end - comm%t
else if (comm%dir*((comm%t+two*comm%h)-comm%t_end) >= zero) then
comm%h = (comm%t_end-comm%t)/two
end if
!
! When the integrator is at T and attempts a step of H, the function
! defining the differential equations will be evaluated at a number of
! arguments between T and T+H. If H is too small, these arguments cannot
! be clearly distinguished in the precision available.
!
hmin = max(comm%sqtiny,comm%toosml*max(abs(comm%t),abs(comm%t+comm%h)))
if (abs(comm%h)<hmin) then
ier = 5; nrec = 3; write (comm%rec,"(a/a,e13.5,a,e13.5,a/a)") &
" ** In order to satisfy your error requirements STEP_INTEGRATE would have",&
" ** to use a step size of ",comm%H," at T_NOW = ",comm%T," This is too",&
" ** small for the machine precision."
exit body
end if
!
! Monitor the impact of output on the efficiency of the integration.
!
if (comm%chkeff) then
comm%hit_t_end_count = comm%hit_t_end_count + 1
if (comm%hit_t_end_count>=100 .and. &
comm%hit_t_end_count>=comm%step_count/3) then
ier = 2; nrec = 5; write (comm%rec,"(a/a/a/a/a)") &
" ** More than 100 output points have been obtained by integrating to T_END.",&
" ** They have been sufficiently close to one another that the efficiency",&
" ** of the integration has been degraded. It would probably be (much) more",&
" ** efficient to obtain output by interpolating with INTERPOLATE (after",&
" ** changing to METHOD='M' if you are using METHOD = 'H')."
comm%hit_t_end_count = 0; exit body
end if
end if
!
! Check for stiffness and for too much work. Stiffness can be
! checked only after a successful step.
!
if (.not.failed) then
!
! Check for too much work.
toomch = comm%f_count > max_f_count
if (toomch) then
ier = 3; nrec = 3; write (comm%rec,"(a,i6,a/a/a)") &
" ** Approximately ",max_f_count," function evaluations have been used to",&
" ** compute the solution since the integration started or since this", &
" ** message was last printed."
!
! After this warning message, F_COUNT is reset to permit the integration
! to continue. The total number of function evaluations in the primary
! integration is FULL_F_COUNT + F_COUNT
!
comm%full_f_count = comm%full_f_count + comm%f_count
comm%f_count = 0
end if
!
! Check for stiffness. If stiffness is detected, the message about too
! much work is augmented inside STIFF to explain that it is due to
! stiffness.
!
call stiff_r1(comm,f,toomch,sure_stiff)
if (sure_stiff) then
!
! Predict how much extra work will be needed to reach TND.
extra_wk = (comm%cost*abs((comm%t_end-comm%t)/comm%h_average)) / &
real(comm%full_f_count+comm%f_count,kind=wp)
ier = 4; nrec = nrec + 4
write (comm%rec(nrec-3:nrec),"(a/a,e13.5,a/a/a)") &
" ** Your problem has been diagnosed as stiff. If the situation persists,",&
" ** it will cost roughly ",extra_wk," times as much to reach T_END as it", &
" ** has cost to reach T_NOW. You should probably change to a code intended",&
" ** for stiff problems."
end if
if (ier/=just_fine) exit body
end if
!
! Take a step. Whilst finding an on-scale H (PHASE2 = .TRUE.), the input
! value of H might be reduced (repeatedly), but it will not be reduced
! below HMIN. The local error is estimated, a weight vector is formed,
! and a weighted maximum norm, ERR, of the local error is returned.
! The presence of the optional argument PHASE2 in the call to STEP
! indicates that this is the primary integration.
!
! H is used by both STEP_INTEGRATE and STEP. Since it may be changed inside
! STEP, a local copy is made.
!
htry = comm%h
call step_r1(comm,f,comm%t,comm%y,comm%yp,comm%stages,comm%tol,htry, &
comm%y_new,comm%err_estimates,err,hmin,comm%phase2)
comm%h = htry
!
! Compare the norm of the local error to the tolerance.
!
if (err > comm%tol) then
!
! Failed step. Reduce the step size and try again.
!
! First step: Terminate PHASE3 of the search for an on-scale step size.
! The step size is not on scale, so ERR may not be accurate;
! reduce H by a fixed factor. Failed attempts to take the
! first step are not counted.
! Later step: Use ERR to compute an "optimal" reduction of H. More than
! one failure indicates a difficulty with the problem and an
! ERR that may not be accurate, so reduce H by a fixed factor.
!
if (comm%at_t_start) then
phase3 = .false.; alpha = comm%rs1
else
comm%bad_step_count = comm%bad_step_count + 1
comm%stiff_bad_step_count = comm%stiff_bad_step_count + 1
if (failed) then
alpha = comm%rs1
else
alpha = comm%safety*(comm%tol/err)**comm%expon
alpha = max(alpha,comm%rs1)
end if
end if
comm%h = alpha*comm%h; failed = .true.; cycle take_step
end if
!
! Successful step.
!
! Predict a step size appropriate for the next step. After the first
! step the prediction can be refined using an idea of H.A. Watts that
! takes account of how well the prediction worked on the previous step.
!
beta = (err/comm%tol)**comm%expon
if (.not.comm%at_t_start) then
t1 = (err**comm%expon)/comm%h
t2 = (comm%errold**comm%expon)/comm%h_old
if (t1<t2*hundrd .and. t2<t1*hundrd) beta = beta*(t1/t2)
end if
alpha = comm%rs3
if (comm%safety < beta*alpha) alpha = comm%safety/beta
!
! On the first step a search is made for an on-scale step size. PHASE2
! of the scheme comes to an end here because a step size has been found
! that is both successful and has a credible local error estimate. Except
! in the special case that the first step is also the last, the step is
! repeated in PHASE3 as long as an increase greater than RS2 appears
! possible. An increase as big as RS3 is permitted. A step failure
! terminates PHASE3.
!
if (comm%at_t_start) then
comm%phase2 = .false.
phase3 = phase3 .and. .not. comm%at_t_end .and. (alpha > comm%rs2)
if (phase3) then
comm%h = alpha*comm%h; cycle take_step
end if
end if
!
! After getting on scale, step size changes are more restricted.
!
alpha = min(alpha,comm%rs)
if (failed) alpha = min(alpha,one)
alpha = max(alpha,comm%rs1)
comm%h_old = comm%h; comm%h = alpha*comm%h
!
! For the diagnosis of stiffness, an average accepted step size, H_AVERAGE,
! must be computed.
!
if (comm%at_t_start) then
comm%h_average = comm%h_old
else
comm%h_average = pt9*comm%h_average + pt1*comm%h_old
end if
!
comm%at_t_start = .false.; comm%errold = err; comm%t_old = comm%t
!
! Take care that T is set to precisely T_END when the end of the
! integration is reached.
!
if (comm%at_t_end) then
comm%t = comm%t_end
else
comm%t = comm%t + comm%h_old
end if
!
! Increment counter on accepted steps. Note that successful steps
! that are repeated whilst getting on scale are not counted.
!
comm%step_count = comm%step_count + 1
!
! Advance the current solution and its derivative. Note that the previous
! derivative will overwrite Y_NEW (see pointer assignments in SETUP).
!
comm%y_old = comm%y; comm%y = comm%y_new
comm%yp_old = comm%yp
!
if (comm%fsal) then
!
! When FSAL = .TRUE., YP is the last stage of the step.
!
comm%yp = comm%stages(:,comm%last_stage)
else
!
! Call F to evaluate YP.
!
comm%yp = f(comm%t,comm%y); comm%f_count = comm%f_count + 1
end if
!
! If global error assessment is desired, advance the secondary
! integration from TOLD to T.
!
if (comm%erason) then
call truerr_r1(comm,f,ier)
if (ier==6) then
!
! The global error estimating procedure has broken down. Treat it as a
! failed step. The solution and derivative are reset to their values at
! the beginning of the step since the last valid error assessment refers
! to them.
!
comm%step_count = comm%step_count - 1; comm%erasfl = .true.
comm%at_t_end = .false.
comm%t = comm%t_old; comm%h = comm%h_old
comm%y = comm%y_old; comm%yp = comm%yp_old
if (comm%step_count > 0) then
nrec = 2; write (comm%rec,"(a/a,e13.5/a)") &
" ** The global error assessment may not be reliable for T past ",&
" ** T_NOW = ",comm%t,". The integration is being terminated."
exit body
else
nrec = 2; write (comm%rec,"(a/a)") &
" ** The global error assessment algorithm failed at the start of the ",&
" ** integration. The integration is being terminated."
exit body
end if
end if
end if
exit take_step
end do take_step
exit body
end do body
!
! Exit point for STEP_INTEGRATE
! Set the output variables and flag that interpolation is permitted
!
if (ier < fatal) then
t_now = comm%t; comm%at_t_end = t_now == comm%t_end
comm%chkeff = comm%at_t_end;
y_now = comm%y; yderiv_now = comm%yp
comm%ymax = max(abs(comm%y),comm%ymax)
if (ier==just_fine) then
state = usable; call set_saved_state_r1("INTERPOLATE",state,comm)
end if
end if
!
! Call RKMSG_R1 to report what happened and set FLAG.
!
call rkmsg_r1(ier,srname,nrec,comm,flag)
!
end subroutine step_integrate_r1