recursive subroutine range_integrate_r1(comm,f,t_want,t_got,y_got,yderiv_got, &
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(in) :: t_want !indep!
real(kind=wp), intent(out) :: t_got !indep!
real(kind=wp), dimension(:), intent(out) :: y_got, yderiv_got !dep!
integer, intent(out), optional :: flag
type(rk_comm_real_1d), intent(inout) :: comm
!
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="RANGE_INTEGRATE"
!
real(kind=wp) :: hmin, t_now !indep!
integer :: step_flag, ier, nrec, state
logical :: goback, baderr
!
integer, parameter :: not_ready=-1, usable=-2, fatal=911, catastrophe=912, &
just_fine=1
logical, parameter :: tell=.false., ask=.true.
real(kind=wp), parameter :: zero=0.0_wp
!
ier = just_fine; nrec = 0
goback = .false.; baderr = .false.
body: do
!
! Is it permissible to call RANGE_INTEGRATE?
!
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 RANGE_INTEGRATE."
exit body
end if
if (.not.comm%use_range) then
ier = fatal; nrec = 2; write (comm%rec,"(a/a)") &
" ** You have called RANGE_INTEGRATE after you specified in SETUP that you",&
" ** were going to use STEP_INTEGRATE. This is not permitted."
exit body
end if
state = get_saved_state_r1(srname,comm%save_states)
if (state==5 .or. state==6) then
ier = fatal; nrec = 1; write (comm%rec,"(a/a)") &
" ** This routine has already returned with a hard failure. You must call",&
" ** SETUP to start another problem."
exit body
end if
state = usable
call set_saved_state_r1(srname,state,comm)
!
if (comm%at_t_start) then
!
! First call.
!
! A value of T_END is specified in SETUP. When INTRP_ABLE = .FALSE., as with
! METHOD = 'H', output is obtained at the specified T_WANT by resetting T_END
! to T_WANT. At this point, before the integration gets started, this can
! be done with a simple assignment. Later it is done with a call to
! RESET_T_END. The original T_END is saved in RANGE_T_END.
!
comm%range_t_end = comm%t_end
if (.not.comm%intrp_able) comm%t_end = t_want
!
! The last T_GOT returned is in the variable TLAST. T records how far the
! integration has advanced towards the specified T_END. When output is
! obtained by interpolation, the integration goes past the T_GOT returned
! (T is closer to the specified T_END than T_GOT).
comm%tlast = comm%t_start; t_got = comm%t_start
!
! If the code is to find an on-scale initial step size H, a bound was placed
! on H in SETUP. Here the first output point is used to refine this bound.
if (comm%h_start==zero) then
comm%h = min(abs(comm%h),abs(t_want-comm%t_start))
hmin = max(comm%sqtiny,comm%toosml* &
max(abs(comm%t_start),abs(comm%t_end)))
comm%h = max(comm%h,hmin)
end if
!
else
!
! Subsequent call.
!
if (comm%tlast==comm%range_t_end) then
ier = fatal; nrec = 3; write (comm%rec,"(a/a/a)") &
" ** You have called RANGE_INTEGRATE after reaching T_END. (Your last call",&
" ** to RANGE_INTEGRATE resulted in T_GOT = T_END.) To start a new",&
" ** problem, you will need to call SETUP."
exit body
end if
!
end if
!
! Check for valid T_WANT.
!
if (comm%dir*(t_want-comm%tlast)<=zero) then
ier = fatal; nrec = 3; write (comm%rec,"(a/a/a)") &
" ** You have made a call to RANGE_INTEGRATE with a T_WANT that does not lie",&
" ** between the previous value of T_GOT (T_START on the first call) and",&
" ** T_END. This is not permitted. Check your program carefully."
exit body
end if
if (comm%dir*(t_want-comm%range_t_end)>zero) then
hmin = max(comm%sqtiny,comm%toosml* &
max(abs(t_want),abs(comm%range_t_end)))
if (abs(t_want-comm%range_t_end)<hmin) then
ier = fatal; nrec = 4; write (comm%rec,"(a/a/a/a)") &
" ** You have made a call to RANGE_INTEGRATE with a T_WANT that does not lie",&
" ** between the previous value of T_GOT (T_START on the first call) and",&
" ** T_END. This is not permitted. T_WANT is very close to T_END, so you may",&
" ** have meant to set it to be T_END exactly. Check your program carefully."
else
ier = fatal; nrec = 3; write (comm%rec,"(a/a/a/a)") &
" ** You have made a call to RANGE_INTEGRATE with a T_WANT that does not lie",&
" ** between the previous value of T_GOT (T_START on the first call) and",&
" ** T_END. This is not permitted. Check your program carefully."
end if
exit body
end if
if (.not.comm%intrp_able) then
hmin = max(comm%sqtiny,comm%toosml*max(abs(comm%tlast),abs(t_want)))
if (abs(t_want-comm%tlast)<hmin) then
ier = fatal; nrec = 4; write (comm%rec,"(a/a/a/a,e13.5,a)") &
" ** You have made a call to RANGE_INTEGRATE with a T_WANT that is not",&
" ** sufficiently different from the last value of T_GOT (T_START on the",&
" ** first call). When using METHOD = 'H', it must differ by at least ",&
" ** ",HMIN,"."
exit body
end if
!
! We have a valid T_WANT. There is no interpolation with this METHOD and
! therefore we step to T_WANT exactly by resetting T_END with a call to
! RESET_T_END. On the first step this matter is handled differently as
! explained above.
!
if (.not.comm%at_t_start) then
call reset_t_end(comm,t_want)
baderr = get_saved_fatal_r1(comm)
if (baderr) exit body
end if
end if
!
! Process output, decide whether to take another step.
!
proceed: do
!
if (comm%intrp_able) then
!
! Interpolation is possible with this METHOD. The integration has
! already reached T. If this is past T_WANT, GOBACK is set .TRUE. and
! the answers are obtained by interpolation.
!
goback = comm%dir*(comm%t-t_want) >= zero
if (goback) then
call interpolate(comm,f,t_want,y_got,yderiv_got)
baderr = get_saved_fatal_r1(comm)
if (baderr) exit body
t_got = t_want
end if
else
!
! Interpolation is not possible with this METHOD, so output is obtained
! by integrating to T_WANT = T_END. Both Y_GOT and YDERIV_GOT are then
! already loaded with the solution at T_WANT by STEP_INTEGRATE.
!
goback = comm%t == t_want
if (goback) t_got = t_want
end if
!
! If done, go to the exit point.
if (goback) exit body
!
! Take a step with STEP_INTEGRATE in the direction of T_END. On exit, the
! solution is advanced to T_NOW. The approximate solution at T_NOW is
! available in Y_GOT. If output is obtained by stepping to the end (T_NOW
! = T_WANT = T_END), Y_GOT can be returned directly. If output is
! obtained by interpolation, the subroutine INTERPOLATE that does this uses
! the values in COMM for its computations and places the approximate solution
! at T_WANT in the arrays Y_GOT,YDERIV_GOT for return to the calling
! program. T_NOW is output from STEP_INTEGRATE and is actually a copy of T
! from inside COMM.
call step_integrate(comm,f,t_now,y_got,yderiv_got,step_flag)
ier = step_flag
!
! A successful step by STEP_INTEGRATE is indicated by step_flag= 1.
!
select case(step_flag)
case(1); cycle proceed
case(2); nrec = 4; write (comm%rec,"(a/a/a/a)") &
" ** The last message was produced on a call to STEP_INTEGRATE from",&
" ** RANGE_INTEGRATE. In RANGE_INTAGRATE the appropriate action is to",&
" ** change to METHOD = 'M', or, if insufficient memory is available,",&
" ** to METHOD = 'L'. "
case(3:6); nrec = 2; write (comm%rec,"(a)") &
" ** The last message was produced on a call to STEP_INTEGRATE from",&
" ** RANGE_INTEGRATE."
case default; baderr = .true.
end select
t_got = comm%t; exit body
end do proceed
!
end do body
!
if (baderr) then
ier = fatal; nrec = 3; write (comm%rec,"(a/a/a)") &
" ** An internal call by RANGE_INTEGRATE to a subroutine resulted in an",&
" ** error that should not happen. Check your program carefully for array",&
" ** sizes, correct number of arguments, type mismatches ... ."
end if
!
comm%tlast = t_got
!
! All exits are done here after a call to RKMSG_R1 to report
! what happened
!
call rkmsg_r1(ier,srname,nrec,comm,flag)
!
end subroutine range_integrate_r1