subroutine truerr_r1(comm,f,ier)
!
! 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) :: comm
integer, intent(inout) :: ier
!
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) :: hmin, hsec !indep!
real(kind=wp) :: diff, errmax, mxerlc, tsec, ge_err, ge_test1, ge_test2
integer :: istep, level
!
integer, parameter :: just_fine=1
real(kind=wp), parameter :: pt1=0.1_wp, ten=10.0_wp
real(kind=wp), dimension(:,:), pointer :: ge_stages !dep!
real(kind=wp), dimension(:), pointer :: ge_y, ge_yp, ge_y_new !dep!
real(kind=wp), dimension(:), pointer :: ge_err_estimates, y !dep!
real(kind=wp), dimension(:), pointer :: ge_assess, weights !shp-dep!
!
ge_stages => comm%ge_stages
ge_y => comm%ge_y
ge_yp => comm%ge_yp
ge_y_new => comm%ge_y_new
ge_err_estimates => comm%ge_err_estimates
ge_assess => comm%ge_assess
y => comm%y
weights => comm%weights
!
tsec = comm%t - comm%h_old
hsec = comm%h_old/real(comm%no_of_ge_steps,kind=wp)
hmin = max(comm%sqtiny,comm%toosml*max(abs(tsec),abs(comm%t)))
body: do
if (abs(hsec)<hmin) then
ier = 6; exit body
end if
ge_test1 = comm%tol/real(comm%no_of_ge_steps,kind=wp)
ge_test2 = comm%tol/ten; level = 0
!
! The subroutine STEP is used to take a step.
!
! Perform secondary integration.
!
do istep = 1, comm%no_of_ge_steps
!
! Take a step.
call step_r1(comm,f,tsec,ge_y,ge_yp,ge_stages,ge_test1,hsec,ge_y_new, &
ge_err_estimates,ge_err)
!
! The primary integration is using a step size of H_OLD and the
! secondary integration is using the smaller step size
! HSEC = H_OLD/(NO_OF_GE_STEPS).
! If steps of this size were taken from the same starting point and the
! asymptotic behavior were evident, the smaller step size would result in
! a local error that is considerably smaller, namely by a factor of
! 1/(NO_OF_GE_STEPSSEC**(ORDER+1)). If the two approximate solutions are
! close and TOL is neither too large nor too small, this should be
! approximately true. The step size is chosen in the primary integration
! so that the local error ERR is no larger than TOL. The local error,
! GE_ERR, of the secondary integration is compared to TOL in an attempt to
! diagnose a secondary integration that is not rather more accurate than
! the primary integration.
!
if (ge_err>=ge_test1) then
level = 2
else if (ge_err>ge_test2) then
level = level + 1
end if
if (level>=2) then
ier = 6; exit body
end if
!
! Advance TSEC and the dependent variables GE_Y and GE_YP.
!
tsec = comm%t - real(comm%no_of_ge_steps-istep,kind=wp)*hsec
ge_y = ge_y_new
!
if (comm%fsal) then
!
! When FSAL = .TRUE., the derivative GE_YP is the last stage of the step.
!
ge_yp = ge_stages(:,comm%last_stage)
else
!
! Call F to evaluate GE_YP.
!
ge_yp = f(tsec,ge_y); comm%ge_f_count = comm%ge_f_count + 1
end if
!
end do
!
! Update the maximum error seen, GE_MAX_CONTRIB, and its location,
! T_GE_MAX_CONTRIB. Use local variables ERRMAX and MXERLC.
!
errmax = comm%ge_max_contrib; mxerlc = comm%t_ge_max_contrib;
diff = maxval(abs(ge_y-y)/weights) !spec-ar!
if (diff>errmax) then
errmax = diff; mxerlc = comm%t
end if
!
! If the global error is greater than 0.1, the solutions have diverged so
! far that comparing them may not provide a reliable estimate of the global
! error. The test is made before GE_ASSESS and GE_MAX_CONTRIB,
! T_GE_MAX_CONTRIB are updated so that on a failure, they refer to the
! last reliable results.
!
if (errmax>pt1) then
ier = 6
else
comm%ge_max_contrib = errmax; comm%t_ge_max_contrib = mxerlc;
ge_assess = ge_assess + (abs(ge_y-y)/weights)**2
ier = just_fine
end if
exit body
!
end do body
!
end subroutine truerr_r1