! This set of 10 Fortran 90 programs are also from NAG site
! The first one is basic Ronger Kuta solver
! Following 9 are examples

module rksuite_90_prec
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
!
integer, parameter :: wp = selected_real_kind(10,50)

end module rksuite_90_prec

module rksuite_90
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
!
use rksuite_90_prec, only:wp

implicit none

private

public :: wp, setup, range_integrate, step_integrate, interpolate, &
          global_error, statistics, reset_t_end, collect_garbage,  &
          set_stop_on_fatal, get_saved_fatal
!starta!
public :: rk_comm_real_1d
type rk_comm_real_1d

private

real(kind=wp) :: t, t_old, t_start, t_end, dir                       !indep!
real(kind=wp) :: h, h_old, h_start, h_average                        !indep!
real(kind=wp) :: tol
integer :: f_count, full_f_count, step_count, bad_step_count
logical :: at_t_start, at_t_end


real(kind=wp), dimension(:), pointer :: thresh, weights, ymax        !shp-dep!

real(kind=wp), dimension(:), pointer :: scratch, y, yp, y_new        !dep!
real(kind=wp), dimension(:), pointer :: y_old, yp_old, v0, v1        !dep!
real(kind=wp), dimension(:), pointer :: err_estimates, v2, v3        !dep!
real(kind=wp), dimension(:), pointer ::  vtemp                       !dep!
real(kind=wp), dimension(:,:), pointer :: stages                     !dep!

real(kind=wp) :: a(13,13), b(13), c(13), bhat(13), r(11,6), e(7)
integer :: ptr(13), no_of_stages, rk_method, intrp_degree
logical :: intrp_able, intrp_needs_stages

real(kind=wp) :: toosml, cost, safety, expon, stability_radius, tan_angle, &
    rs, rs1, rs2, rs3, rs4
integer :: order, last_stage, max_stiff_iters, no_of_ge_steps
logical :: fsal

real(kind=wp) :: ge_max_contrib
real(kind=wp) :: t_ge_max_contrib                                    !indep!
integer :: ge_f_count
real(kind=wp), dimension(:), pointer :: ge_assess                    !shp-dep!

real(kind=wp), dimension(:), pointer :: ge_y, ge_yp, ge_y_new        !dep!
real(kind=wp), dimension(:), pointer :: ge_err_estimates             !dep!
real(kind=wp), dimension(:,:), pointer :: ge_stages                  !dep!

logical :: erason, erasfl

real(kind=wp) :: mcheps, dwarf, round_off, sqrrmc, cubrmc, sqtiny
integer :: outch

logical :: print_message, use_range

character(len=80) :: rec(10)

real(kind=wp) :: tlast, range_t_end                                  !indep!

real(kind=wp), dimension(:), pointer :: xstage, ytemp                !dep!
real(kind=wp), dimension(:,:), pointer :: p                          !dep!

integer :: stiff_bad_step_count, hit_t_end_count
real(kind=wp) :: errold
logical :: chkeff, phase2

integer, dimension(7) :: save_states

logical :: stop_on_fatal, saved_fatal_err

end type rk_comm_real_1d
!enda!
interface setup
   module procedure setup_r1
end interface
interface range_integrate
   module procedure range_integrate_r1
end interface
interface step_integrate
   module procedure step_integrate_r1
end interface
interface statistics
   module procedure statistics_r1
end interface
interface global_error
   module procedure global_error_r1
end interface
interface reset_t_end
   module procedure reset_t_end_r1
end interface
interface interpolate
   module procedure interpolate_r1
end interface
interface set_stop_on_fatal
   module procedure set_stop_on_fatal_r1
end interface
interface get_saved_fatal
   module procedure get_saved_fatal_r1
end interface
interface collect_garbage
   module procedure collect_garbage_r1
end interface

contains
!startb!
subroutine machine_const(round_off,sqrrmc,cubrmc,sqtiny,outch)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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) :: round_off, sqrrmc, cubrmc, sqtiny
integer, intent(out) :: outch
!
real(kind=wp) :: dummy
real(kind=wp), parameter :: third=1.0_wp/3.0_wp, ten=10.0_wp
!
outch = 6
!
round_off = ten*epsilon(dummy)
sqrrmc = sqrt(epsilon(dummy))
cubrmc = epsilon(dummy)**third
sqtiny = sqrt(tiny(dummy))
!
end subroutine machine_const

subroutine method_const(rk_method, a, b, c, bhat, r, e, ptr, no_of_stages, &
                        intrp_degree, intrp_able, intrp_needs_stages, &
                        cost, safety, expon, stability_radius, &
                        tan_angle, rs, rs1, rs2, rs3, rs4, order, last_stage, &
                        max_stiff_iters, no_of_ge_steps, fsal, cdiff)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
!
integer, intent(in) :: rk_method
real(kind=wp), intent(out) :: a(13,13), b(13), c(13), bhat(13), r(11,6), e(7)
integer, intent(out) :: ptr(13), no_of_stages, intrp_degree
logical, intent(out) :: intrp_able, intrp_needs_stages

real(kind=wp), intent(out) :: cost, safety, expon, stability_radius, &
    tan_angle, rs, rs1, rs2, rs3, rs4, cdiff
integer, intent(out) :: order, last_stage, max_stiff_iters, no_of_ge_steps
logical, intent(out) :: fsal
!
integer :: i
real(kind=wp), parameter :: fivepc=0.05_wp,  one=1.0_wp, two=2.0_wp, &
   fifty=50.0_wp
!
select case (rk_method)
case(1)
!
!  METHD = 1.
!    This pair is from "A 3(2) Pair of Runge-Kutta Formulas" by P. Bogacki
!    and L.F. Shampine, Appl. Math. Lett., 2, pp. 321-325, 1989.  The authors
!    are grateful to P. Bogacki for his assistance in implementing the pair.
!
   no_of_stages = 4; fsal = .true.; order = 2
   tan_angle = 8.9_wp; stability_radius = 2.3_wp
   safety = 0.8_wp; intrp_able = .true.; intrp_degree = 3
   intrp_needs_stages = .false.; no_of_ge_steps = 3
!
   ptr(1:4) = (/ 0, 1, 2, 3 /)
!
   a(2,1) = 1.0_wp/2.0_wp
   a(3,1) = 0.0_wp
   a(3,2) = 3.0_wp/4.0_wp
   a(4,1) = 2.0_wp/9.0_wp
   a(4,2) = 1.0_wp/3.0_wp
   a(4,3) = 4.0_wp/9.0_wp
!
!  The coefficients BHAT refer to the formula used to advance the
!  integration, here the one of order 3.  The coefficients B refer
!  to the other formula, here the one of order 2. For this pair, BHAT
!  is not needed since FSAL = .TRUE.
!
   b(1) = 7.0_wp/24.0_wp
   b(2) = 1.0_wp/4.0_wp
   b(3) = 1.0_wp/3.0_wp
   b(4) = 1.0_wp/8.0_wp
!
   c(1) = 0.0_wp
   c(2) = 1.0_wp/2.0_wp
   c(3) = 3.0_wp/4.0_wp
   c(4) = 1.0_wp
!
case (2)
!
!  METHD = 2
!    This pair is from "An Efficient Runge-Kutta (4,5) Pair" by P. Bogacki
!    and L.F. Shampine, Rept. 89-20, Math. Dept., Southern Methodist
!    University, Dallas, Texas, USA, 1989.  The authors are grateful to
!    P. Bogacki for his assistance in implementing the pair.  Shampine and
!    Bogacki subsequently modified the formula to enhance the reliability of
!    the pair.  The original fourth order formula is used in an estimate of
!    the local error.  If the step fails, the computation is broken off.  If
!    the step is acceptable, the first evaluation of the next step is done,
!    i.e., the pair is implemented as FSAL and the local error of the step
!    is again estimated with a fourth order formula using the additional data.
!    The step must succeed with both estimators to be accepted.  When the
!    second estimate is formed, it is used for the subsequent adjustment of
!    the step size because it is of higher quality.  The two fourth order
!    formulas are well matched to leading order, and only exceptionally do
!    the estimators disagree -- problems with discontinuous coefficients are
!    handled more reliably by using two estimators as is global error
!    estimation.
!
   no_of_stages = 8; fsal = .true.; order = 4
   tan_angle = 5.2_wp; stability_radius = 3.9_wp
   safety = 0.8_wp; intrp_able = .true.
   intrp_needs_stages = .true.; intrp_degree = 6
   no_of_ge_steps = 2
!
   ptr(1:8) = (/ 0, 1, 2, 3, 4, 5, 6, 7 /)
!
   a(2,1) = 1.0_wp/6.0_wp
   a(3,1) = 2.0_wp/27.0_wp
   a(3,2) = 4.0_wp/27.0_wp
   a(4,1) = 183.0_wp/1372.0_wp
   a(4,2) = -162.0_wp/343.0_wp
   a(4,3) = 1053.0_wp/1372.0_wp
   a(5,1) = 68.0_wp/297.0_wp
   a(5,2) = -4.0_wp/11.0_wp
   a(5,3) = 42.0_wp/143.0_wp
   a(5,4) = 1960.0_wp/3861.0_wp
   a(6,1) = 597.0_wp/22528.0_wp
   a(6,2) = 81.0_wp/352.0_wp
   a(6,3) = 63099.0_wp/585728.0_wp
   a(6,4) = 58653.0_wp/366080.0_wp
   a(6,5) = 4617.0_wp/20480.0_wp
   a(7,1) = 174197.0_wp/959244.0_wp
   a(7,2) = -30942.0_wp/79937.0_wp
   a(7,3) = 8152137.0_wp/19744439.0_wp
   a(7,4) = 666106.0_wp/1039181.0_wp
   a(7,5) = -29421.0_wp/29068.0_wp
   a(7,6) = 482048.0_wp/414219.0_wp
   a(8,1) = 587.0_wp/8064.0_wp
   a(8,2) = 0.0_wp
   a(8,3) = 4440339.0_wp/15491840.0_wp
   a(8,4) = 24353.0_wp/124800.0_wp
   a(8,5) = 387.0_wp/44800.0_wp
   a(8,6) = 2152.0_wp/5985.0_wp
   a(8,7) = 7267.0_wp/94080.0_wp
!
!  The coefficients B refer to the formula of order 4.
!
   b(1) = 2479.0_wp/34992.0_wp
   b(2) = 0.0_wp
   b(3) = 123.0_wp/416.0_wp
   b(4) = 612941.0_wp/3411720.0_wp
   b(5) = 43.0_wp/1440.0_wp
   b(6) = 2272.0_wp/6561.0_wp
   b(7) = 79937.0_wp/1113912.0_wp
   b(8) = 3293.0_wp/556956.0_wp
!
!  The coefficients E refer to an estimate of the local error based on
!  the first formula of order 4.  It is the difference of the fifth order
!  result, here located in A(8,:), and the fourth order result.  By
!  construction both E(2) and E(7) are zero.
!
   e(1) = -3.0_wp/1280.0_wp
   e(2) = 0.0_wp
   e(3) = 6561.0_wp/632320.0_wp
   e(4) = -343.0_wp/20800.0_wp
   e(5) = 243.0_wp/12800.0_wp
   e(6) = -1.0_wp/95.0_wp
   e(7) = 0.0_wp
!
   c(1) = 0.0_wp
   c(2) = 1.0_wp/6.0_wp
   c(3) = 2.0_wp/9.0_wp
   c(4) = 3.0_wp/7.0_wp
   c(5) = 2.0_wp/3.0_wp
   c(6) = 3.0_wp/4.0_wp
   c(7) = 1.0_wp
   c(8) = 1.0_wp
!
!  To do interpolation with this pair, some extra stages have to be computed.
!  The following additional A and C coefficients are for this purpose.
!  In addition there is an array R that plays a role for interpolation
!  analogous to that of BHAT for the basic step.
!
   c(9) = 1.0_wp/2.0_wp
   c(10) = 5.0_wp/6.0_wp
   c(11) = 1.0_wp/9.0_wp
!
   a(9,1) = 455.0_wp/6144.0_wp
   a(10,1) = -837888343715.0_wp/13176988637184.0_wp
   a(11,1) = 98719073263.0_wp/1551965184000.0_wp
   a(9,2) = 0.0_wp
   a(10,2) = 30409415.0_wp/52955362.0_wp
   a(11,2) = 1307.0_wp/123552.0_wp
   a(9,3) = 10256301.0_wp/35409920.0_wp
   a(10,3) = -48321525963.0_wp/759168069632.0_wp
   a(11,3) = 4632066559387.0_wp/70181753241600.0_wp
   a(9,4) = 2307361.0_wp/17971200.0_wp
   a(10,4) = 8530738453321.0_wp/197654829557760.0_wp
   a(11,4) = 7828594302389.0_wp/382182512025600.0_wp
   a(9,5) = -387.0_wp/102400.0_wp
   a(10,5) = 1361640523001.0_wp/1626788720640.0_wp
   a(11,5) = 40763687.0_wp/11070259200.0_wp
   a(9,6) = 73.0_wp/5130.0_wp
   a(10,6) = -13143060689.0_wp/38604458898.0_wp
   a(11,6) = 34872732407.0_wp/224610586200.0_wp
   a(9,7) = -7267.0_wp/215040.0_wp
   a(10,7) = 18700221969.0_wp/379584034816.0_wp
   a(11,7) = -2561897.0_wp/30105600.0_wp
   a(9,8) = 1.0_wp/32.0_wp
   a(10,8) = -5831595.0_wp/847285792.0_wp
   a(11,8) = 1.0_wp/10.0_wp
   a(10,9) = -5183640.0_wp/26477681.0_wp
   a(11,9) = -1.0_wp/10.0_wp
   a(11,10) = -1403317093.0_wp/11371610250.0_wp
!
   r(1:11,1) = 0.0_wp; r(2,1:6) = 0.0_wp
   r(1,6) = -12134338393.0_wp/1050809760.0_wp
   r(1,5) = -1620741229.0_wp/50038560.0_wp
   r(1,4) = -2048058893.0_wp/59875200.0_wp
   r(1,3) = -87098480009.0_wp/5254048800.0_wp
   r(1,2) = -11513270273.0_wp/3502699200.0_wp
!
   r(3,6) = -33197340367.0_wp/1218433216.0_wp
   r(3,5) = -539868024987.0_wp/6092166080.0_wp
   r(3,4) = -39991188681.0_wp/374902528.0_wp
   r(3,3) = -69509738227.0_wp/1218433216.0_wp
   r(3,2) = -29327744613.0_wp/2436866432.0_wp
!
   r(4,6) = -284800997201.0_wp/19905339168.0_wp
   r(4,5) = -7896875450471.0_wp/165877826400.0_wp
   r(4,4) = -333945812879.0_wp/5671036800.0_wp
   r(4,3) = -16209923456237.0_wp/497633479200.0_wp
   r(4,2) = -2382590741699.0_wp/331755652800.0_wp
!
   r(5,6) = -540919.0_wp/741312.0_wp
   r(5,5) = -103626067.0_wp/43243200.0_wp
   r(5,4) = -633779.0_wp/211200.0_wp
   r(5,3) = -32406787.0_wp/18532800.0_wp
   r(5,2) = -36591193.0_wp/86486400.0_wp
!
   r(6,6) = 7157998304.0_wp/374350977.0_wp
   r(6,5) = 30405842464.0_wp/623918295.0_wp
   r(6,4) = 183022264.0_wp/5332635.0_wp
   r(6,3) = -3357024032.0_wp/1871754885.0_wp
   r(6,2) = -611586736.0_wp/89131185.0_wp
!
   r(7,6) = -138073.0_wp/9408.0_wp
   r(7,5) = -719433.0_wp/15680.0_wp
   r(7,4) = -1620541.0_wp/31360.0_wp
   r(7,3) = -385151.0_wp/15680.0_wp
   r(7,2) = -65403.0_wp/15680.0_wp
!
   r(8,6) = 1245.0_wp/64.0_wp
   r(8,5) = 3991.0_wp/64.0_wp
   r(8,4) = 4715.0_wp/64.0_wp
   r(8,3) = 2501.0_wp/64.0_wp
   r(8,2) = 149.0_wp/16.0_wp
   r(8,1) = 1.0_wp
!
   r(9,6) = 55.0_wp/3.0_wp
   r(9,5) = 71.0_wp
   r(9,4) = 103.0_wp
   r(9,3) = 199.0_wp/3.0_wp
   r(9,2) = 16.0d0
!
   r(10,6) = -1774004627.0_wp/75810735.0_wp
   r(10,5) = -1774004627.0_wp/25270245.0_wp
   r(10,4) = -26477681.0_wp/359975.0_wp
   r(10,3) = -11411880511.0_wp/379053675.0_wp
   r(10,2) = -423642896.0_wp/126351225.0_wp
!
   r(11,6) = 35.0_wp
   r(11,5) = 105.0_wp
   r(11,4) = 117.0_wp
   r(11,3) = 59.0_wp
   r(11,2) = 12.0_wp
!
case (3)
!
!  METHD = 3
!    This pair is from "High Order Embedded Runge-Kutta Formulae" by P.J.
!    Prince and J.R. Dormand, J. Comp. Appl. Math.,7, pp. 67-75, 1981.  The
!    authors are grateful to P. Prince and J. Dormand for their assistance in
!    implementing the pair.
!
   no_of_stages = 13; fsal = .false.; order = 7
   tan_angle = 11.0_wp; stability_radius = 5.2_wp
   safety = 0.8_wp; intrp_able = .false.
   intrp_needs_stages = .false.; intrp_degree = 0
   no_of_ge_steps = 2
!
   ptr(1:13) = (/ 0, 1, 2, 1, 3, 2, 4, 5, 6, 7, 8, 9, 1 /)
!
   a(2,1) = 5.55555555555555555555555555556e-2_wp
   a(3,1) = 2.08333333333333333333333333333e-2_wp
   a(3,2) = 6.25e-2_wp
   a(4,1) = 3.125e-2_wp
   a(4,2) = 0.0_wp
   a(4,3) = 9.375e-2_wp
   a(5,1) = 3.125e-1_wp
   a(5,2) = 0.0_wp
   a(5,3) = -1.171875_wp
   a(5,4) = 1.171875_wp
   a(6,1) = 3.75e-2_wp
   a(6,2) = 0.0_wp
   a(6,3) = 0.0_wp
   a(6,4) = 1.875e-1_wp
   a(6,5) = 1.5e-1_wp
   a(7,1) = 4.79101371111111111111111111111e-2_wp
   a(7,2) = 0.0_wp
   a(7,3) = 0.0_wp
   a(7,4) = 1.12248712777777777777777777778e-1_wp
   a(7,5) = -2.55056737777777777777777777778e-2_wp
   a(7,6) = 1.28468238888888888888888888889e-2_wp
   a(8,1) = 1.6917989787292281181431107136e-2_wp
   a(8,2) = 0.0_wp
   a(8,3) = 0.0_wp
   a(8,4) = 3.87848278486043169526545744159e-1_wp
   a(8,5) = 3.59773698515003278967008896348e-2_wp
   a(8,6) = 1.96970214215666060156715256072e-1_wp
   a(8,7) = -1.72713852340501838761392997002e-1_wp
   a(9,1) = 6.90957533591923006485645489846e-2_wp
   a(9,2) = 0.0_wp
   a(9,3) = 0.0_wp
   a(9,4) = -6.34247976728854151882807874972e-1_wp
   a(9,5) = -1.61197575224604080366876923982e-1_wp
   a(9,6) = 1.38650309458825255419866950133e-1_wp
   a(9,7) = 9.4092861403575626972423968413e-1_wp
   a(9,8) = 2.11636326481943981855372117132e-1_wp
   a(10,1) = 1.83556996839045385489806023537e-1_wp
   a(10,2) = 0.0_wp
   a(10,3) = 0.0_wp
   a(10,4) = -2.46876808431559245274431575997_wp
   a(10,5) = -2.91286887816300456388002572804e-1_wp
   a(10,6) = -2.6473020233117375688439799466e-2_wp
   a(10,7) = 2.84783876419280044916451825422_wp
   a(10,8) = 2.81387331469849792539403641827e-1_wp
   a(10,9) = 1.23744899863314657627030212664e-1_wp
   a(11,1) = -1.21542481739588805916051052503_wp
   a(11,2) = 0.0_wp
   a(11,3) = 0.0_wp
   a(11,4) = 1.66726086659457724322804132886e1_wp
   a(11,5) = 9.15741828416817960595718650451e-1_wp
   a(11,6) = -6.05660580435747094755450554309_wp
   a(11,7) = -1.60035735941561781118417064101e1_wp
   a(11,8) = 1.4849303086297662557545391898e1_wp
   a(11,9) = -1.33715757352898493182930413962e1_wp
   a(11,10) = 5.13418264817963793317325361166_wp
   a(12,1) = 2.58860916438264283815730932232e-1_wp
   a(12,2) = 0.0_wp
   a(12,3) = 0.0_wp
   a(12,4) = -4.77448578548920511231011750971_wp
   a(12,5) = -4.3509301377703250944070041181e-1_wp
   a(12,6) = -3.04948333207224150956051286631_wp
   a(12,7) = 5.57792003993609911742367663447_wp
   a(12,8) = 6.15583158986104009733868912669_wp
   a(12,9) = -5.06210458673693837007740643391_wp
   a(12,10) = 2.19392617318067906127491429047_wp
   a(12,11) = 1.34627998659334941535726237887e-1_wp
   a(13,1) = 8.22427599626507477963168204773e-1_wp
   a(13,2) = 0.0_wp
   a(13,3) = 0.0_wp
   a(13,4) = -1.16586732572776642839765530355e1_wp
   a(13,5) = -7.57622116690936195881116154088e-1_wp
   a(13,6) = 7.13973588159581527978269282765e-1_wp
   a(13,7) = 1.20757749868900567395661704486e1_wp
   a(13,8) = -2.12765911392040265639082085897_wp
   a(13,9) = 1.99016620704895541832807169835_wp
   a(13,10) = -2.34286471544040292660294691857e-1_wp
   a(13,11) = 1.7589857770794226507310510589e-1_wp
   a(13,12) = 0.0_wp
!
!  The coefficients BHAT refer to the formula used to advance the
!  integration, here the one of order 8.  The coefficients B refer
!  to the other formula, here the one of order 7.
!
   bhat(1) = 4.17474911415302462220859284685e-2_wp
   bhat(2) = 0.0_wp
   bhat(3) = 0.0_wp
   bhat(4) = 0.0_wp
   bhat(5) = 0.0_wp
   bhat(6) = -5.54523286112393089615218946547e-2_wp
   bhat(7) = 2.39312807201180097046747354249e-1_wp
   bhat(8) = 7.0351066940344302305804641089e-1_wp
   bhat(9) = -7.59759613814460929884487677085e-1_wp
   bhat(10) = 6.60563030922286341461378594838e-1_wp
   bhat(11) = 1.58187482510123335529614838601e-1_wp
   bhat(12) = -2.38109538752862804471863555306e-1_wp
   bhat(13) = 2.5e-1_wp
!
   b(1) = 2.9553213676353496981964883112e-2_wp
   b(2) = 0.0_wp
   b(3) = 0.0_wp
   b(4) = 0.0_wp
   b(5) = 0.0_wp
   b(6) = -8.28606276487797039766805612689e-1_wp
   b(7) = 3.11240900051118327929913751627e-1_wp
   b(8) = 2.46734519059988698196468570407_wp
   b(9) = -2.54694165184190873912738007542_wp
   b(10) = 1.44354858367677524030187495069_wp
   b(11) = 7.94155958811272872713019541622e-2_wp
   b(12) = 4.44444444444444444444444444445e-2_wp
   b(13) = 0.0_wp
!
   c(1) = 0.0_wp
   c(2) = 5.55555555555555555555555555556e-2_wp
   c(3) = 8.33333333333333333333333333334e-2_wp
   c(4) = 1.25e-1_wp
   c(5) = 3.125e-1_wp
   c(6) = 3.75e-1_wp
   c(7) = 1.475e-1_wp
   c(8) = 4.65e-1_wp
   c(9) = 5.64865451382259575398358501426e-1_wp
   c(10) = 6.5e-1_wp
   c(11) = 9.24656277640504446745013574318e-1_wp
   c(12) = 1.0_wp
   c(13) = c(12)
!
end select
!
!  The definitions of all pairs come here for the calculation of
!  LAST_STAGE - the position of the last evaluated stage in a method
!  RS1, RS2, RS3, RS4 - minimum and maximum rations used is step selection
!  COST - the cost of a step
!  MAX_STIFF_ITERS - the number of iterations permitted in stiffness detection
!     There are at most Q = 3 function calls per iteration. MAX_STIFF_ITERS
!     is determined so that  Q*MAX_STIFF_ITERS <= 5% of the cost of
!     50 steps and 1 <= MAX_STIFF_ITERS <= 8. This limits the number of
!     function calls in each diagnosis of stiffness to 24.
!  EXPON - an exponent for use in step slection
!  CDIFF - a coefficent used in determining the minimum permissible step
!
last_stage = ptr(no_of_stages)
if (fsal) then
   cost = real(no_of_stages-1,kind=wp) !zgs
else
   cost = real(no_of_stages,kind=wp) !zgs
end if
!
max_stiff_iters = min(8,max(1,int(fivepc*cost*fifty)))
!
expon = one/(order+one)
!
!     In calculating CDIFF it is assumed that there will be a non-zero
!     difference |C(I) - C(J)| less than one. If C(I) = C(J) for any I not
!     equal to J, they should be made precisely equal by assignment.
!
cdiff = one
do i = 1, no_of_stages - 1
   cdiff = min( cdiff, minval( &
  abs((c(i)-c(i+1:no_of_stages))), &
       mask=(c(i)-c(i+1:no_of_stages)/=0)) )
end do
!
rs = two; rs1 = one/rs; rs2 = rs**2
rs3 = rs*rs2; rs4 = one/rs3
!
end subroutine method_const
!endb!
!startc!
subroutine setup_r1(comm,t_start,y_start,t_end,tolerance,thresholds, &
                    method,task,error_assess,h_start,message)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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_end, t_start                          !indep!
real(kind=wp), intent(in) :: tolerance
real(kind=wp), dimension(:), intent(in) :: y_start                   !dep!
real(kind=wp), dimension(:), intent(in) :: thresholds                !shp-dep!
type(rk_comm_real_1d) :: comm
real(kind=wp), intent(in), optional :: h_start                       !indep!
logical, intent(in), optional :: error_assess, message
character(len=*), intent(in), optional :: task, method
!
character(len=*), parameter :: srname="SETUP"
!
real(kind=wp) :: hmin                                                !indep!
real(kind=wp) :: cdiff
integer :: ier, nrec, tr_dim_of_stages
logical :: legalt
character(len=1) :: task1, method1
!
integer, parameter :: not_ready=-1, fatal=911, just_fine=1
real(kind=wp), parameter :: zero=0.0_wp, pt01=0.01_wp, fivepc=0.05_wp, &
   third=1.0_wp/3.0_wp, one=1.0_wp, two=2.0_wp, ten=10.0_wp, fifty=50.0_wp
!
ier = just_fine; nrec = 0
!
!  Clear previous state of the suite.
!
call setup_global_stuff
nullify(comm%thresh, comm%err_estimates, comm%weights, comm%y_old, &
   comm%scratch, &
   comm%y, comm%yp, comm%y_new, comm%yp_old, comm%stages, comm%ge_y, &
   comm%ge_yp, comm%ge_err_estimates, comm%ge_assess, comm%ge_y_new, &
   comm%ge_stages, comm%v0, comm%v1, comm%v2, comm%v3, comm%vtemp, &
   comm%xstage, comm%ytemp, comm%p)
!
!  Fetch output channel and machine constants;
!
call machine_const(comm%round_off,comm%sqrrmc,comm%cubrmc,comm%sqtiny, &
                   comm%outch)
!
body: do
!
!  Check for valid shape
   if (size(shape(y_start))>0) then
      if (any(shape(y_start)==0)) then
         ier = fatal; nrec = 2; write (comm%rec,"(a)") &
" ** An extent of Y_START has zero length. This is not permitted."
         exit body
      end if
   end if
!
!  Check and process non-trivial optional arguments
   if (present(task)) then
      task1 = task(1:1); comm%use_range = task1 == "R" .or. task1 == "r"
      legalt = comm%use_range  .or.  task1 == "S" .or. task1 == "s"
      if (.not.legalt) then
         ier = fatal; nrec = 2; write (comm%rec,"(a,a,a/a)") &
" ** You have set the first character of TASK to be '",TASK1,"'.", &
" ** It must be one of 'R','r','S' or 's'."
         exit body
      end if
   end if
   if (present(method)) then
      method1 = method(1:1)
      select case (method1)
      case("L","l"); comm%rk_method = 1
      case("M","m"); comm%rk_method = 2
      case("H","h"); comm%rk_method = 3
      case default
         ier = fatal; nrec = 2; write (comm%rec,"(a,a,a/a)") &
" ** You have set the first character of METHOD to be '",METHOD1,"'.", &
" ** It must be one of 'L','l','M','m','H' or 'h'."
         exit body
      end select
   end if
   if (present(message)) comm%print_message = message
!
! Check consistency of array arguments
!
   if (any(shape(y_start)/=shape(thresholds))) then
      ier = fatal; nrec = 1; write (comm%rec,"(a)") &
" ** The shapes of Y_START and THRESHOLDS are not consistent."
      exit body
   end if
!
! Check and process compulsory arguments
   if (t_start == t_end) then
      ier = fatal; nrec = 1; write (comm%rec,"(a,e13.5,a)") &
" ** You have set T_START = T_END = ",T_START,"."
      exit body
   else
      comm%t_end = t_end; comm%t_start = t_start
      comm%t_old = t_start; comm%t = t_start
      comm%dir = sign(one,t_end-t_start)
   end if
   if ((tolerance > pt01) .or. (tolerance < comm%round_off)) then
      ier = fatal; nrec = 2; write (comm%rec,"(a,e13.5,a/a,e13.5,a)") &
" ** You have set TOLERANCE = ",tolerance," which is not permitted. The", &
" ** range of permitted values is (",comm%round_off,",0.01)."
      exit body
   else
      comm%tol = tolerance
   end if
   if (minval(thresholds) < comm%sqtiny) then                        !spec-ar!
      ier = fatal; nrec = 2; write (comm%rec,"(a/a,e13.5,a)") &
" ** You have set a component of THRESHOLDS to be less than the permitted", &
" ** minimum,",comm%sqtiny,"."
      exit body
   end if
!
!  Set formula definitions and characteristics
   call method_const(comm%rk_method, comm%a, comm%b, comm%c, comm%bhat, &
        comm%r, comm%e, comm%ptr, comm%no_of_stages, comm%intrp_degree, &
        comm%intrp_able, comm%intrp_needs_stages,  comm%cost, &
        comm%safety, comm%expon, comm%stability_radius, comm%tan_angle, &
        comm%rs, comm%rs1, comm%rs2, comm%rs3, comm%rs4, comm%order, &
        comm%last_stage, comm%max_stiff_iters, comm%no_of_ge_steps, comm%fsal,&
        cdiff)
!
   tr_dim_of_stages = maxval(comm%ptr(2:comm%no_of_stages))
   comm%toosml = comm%round_off/cdiff
!
!  In STEP_INTEGRATE the first step taken will have magnitude H.  If
!  H_START = ABS(H_START) is not equal to zero, H = H_START.  If H_START is
!  equal to zero, the code is to find an on-scale initial step size H.  To
!  start this process, H is set here to an upper bound on the first step
!  size that reflects the scale of the independent variable.
!  RANGE_INTEGRATE has some additional information, namely the first output
!  point, that is used to refine this bound in STEP_INTEGRATE when called
!  from RANGE_INTEGRATE.  If H_START is not zero, but it is either too big
!  or too small, the input H_START is ignored and H_START is set to zero to
!  activate the automatic determination of an on-scale initial step size.
!  
   hmin = max(comm%sqtiny,comm%toosml*max(abs(t_start),abs(t_end)))
   if (abs(t_end-t_start) < hmin) then
      ier = fatal; nrec = 4; write (comm%rec,"(a/a/a,e13.5,a/a,e13.5,a)") &
" ** You have set values for T_END and T_START that are not clearly", &
" ** distinguishable for the method and the precision of the computer", &
" ** being used. ABS(T_END-T_START) is ",ABS(T_END-T_START)," but should be", &
" **  at least ",hmin,"."
     exit body
   end if
   if (present(h_start)) comm%h_start = abs(h_start)
   if (comm%h_start > abs(t_end-t_start) .or. comm%h_start < hmin) &
      comm%h_start = zero
   if (comm%h_start == zero) then
      comm%h = max(abs(t_end-t_start)/comm%rs3,hmin)
   else
      comm%h = comm%h_start
   end if
!
!  Allocate a number of arrays using pointers.
!
   allocate(comm%thresh(size(y_start,1)), &                          !alloc!
            comm%err_estimates(size(y_start,1)), &                   !alloc!
            comm%weights(size(y_start,1)), &                         !alloc!
            comm%y_old(size(y_start,1)), &                           !alloc!
            comm%scratch(size(y_start,1)), &                         !alloc!
            comm%y(size(y_start,1)), &                               !alloc!
            comm%yp(size(y_start,1)), &                              !alloc!
            comm%stages(size(y_start,1),tr_dim_of_stages), &         !alloc!
            comm%ymax(size(y_start,1)),stat=ier)                     !alloc!
   if (ier /= 0) then
      ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
      exit body
   else
      comm%y = y_start; comm%ymax = abs(y_start) 
      comm%thresh = thresholds !zgs;
      comm%y_new => comm%scratch; comm%yp_old => comm%scratch
      comm%v0 => comm%err_estimates; comm%vtemp => comm%scratch
      comm%v1 => comm%stages(:,1); comm%v2 => comm%stages(:,2)
      comm%v3 => comm%stages(:,3)
   end if
!
!  Pre-allocate storage for interpolation if the TASK = `R' was specified. 
!
   if (comm%use_range) then
      if (comm%intrp_able) then
         if (comm%rk_method==1) then
            comm%p => comm%stages(:,1:2)
         else if (comm%rk_method==2) then
            allocate(comm%p(size(y_start,1),5), &                    !alloc!
                     comm%ytemp(size(y_start,1)), &                  !alloc!
                     comm%xstage(size(y_start,1)),stat=ier)          !alloc!
            if (ier /= 0) then
               ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
               exit body
            end if
         end if
      end if
   end if
!
!  Initialise state and allocate storage for global error assessment
!
   comm%t_ge_max_contrib = t_start
   if (present(error_assess)) comm%erason = error_assess
   if (comm%erason) then
!
!  Storage is required for the stages of a secondary integration. The
!  stages of the primary intergration can only be overwritten in the
!  cases where there is no interpolant or the interpolant does not
!  require information about the stages (e.g. METHOD 'H' and METHOD 'L',
!  respectively).
      if (.not.comm%intrp_needs_stages) then
         comm%ge_stages => comm%stages
      else
         allocate(comm%ge_stages(size(y_start,1),tr_dim_of_stages),stat=ier) !alloc!
         if (ier /= 0) then
            ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
            exit body
         end if
      end if
      allocate(comm%ge_y(size(y_start,1)), &                         !alloc!
               comm%ge_yp(size(y_start,1)), &                        !alloc!
               comm%ge_err_estimates(size(y_start,1)), &             !alloc!
               comm%ge_assess(size(y_start,1)), &                    !alloc!
               comm%ge_y_new(size(y_start,1)),stat=ier)              !alloc!
      if (ier /= 0) then
         ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
         exit body
      else
         comm%ge_assess = 0.0_wp; comm%ge_y = y_start
      end if
   end if
   exit body
end do body
!
call rkmsg_r1(ier,srname,nrec,comm)
!
contains

subroutine setup_global_stuff
!
comm%h_start=0.0_wp; comm%h_old=0.0_wp
comm%f_count=0; comm%full_f_count=0; comm%step_count=0; comm%bad_step_count=0
comm%at_t_start=.true.;  comm%at_t_end=.false.
comm%rk_method=2 !zgs; 
comm%ge_max_contrib=0.0_wp; comm%ge_f_count=0
comm%erason=.false.;  comm%erasfl=.false.
comm%print_message=.true.;  comm%use_range=.true.
comm%stiff_bad_step_count=0;  comm%hit_t_end_count=0
comm%errold=0.0_wp;  comm%h_average=0.0_wp
comm%chkeff=.false.;  comm%phase2=.true.
comm%save_states(1:7)=not_ready
comm%stop_on_fatal=.true.;  comm%saved_fatal_err=.false.
!
end subroutine setup_global_stuff

end subroutine setup_r1


subroutine collect_garbage_r1(comm)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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) :: comm
!
if (associated(comm%thresh)) then
   deallocate(comm%thresh); nullify(comm%thresh)
end if
if (associated(comm%y)) then
   deallocate(comm%y); nullify(comm%y)
end if
if (associated(comm%yp)) then
   deallocate(comm%yp); nullify(comm%yp)
end if
if (associated(comm%ymax)) then
   deallocate(comm%ymax); nullify(comm%ymax)
end if
if (associated(comm%scratch)) then
   deallocate(comm%scratch); nullify(comm%scratch)
   nullify(comm%y_new); nullify(comm%yp_old); nullify(comm%vtemp)
end if
if (associated(comm%weights)) then
   deallocate(comm%weights); nullify(comm%weights)
end if
if (associated(comm%ytemp)) then
   deallocate(comm%ytemp); nullify(comm%ytemp)
end if
if (associated(comm%y_old)) then
   deallocate(comm%y_old); nullify(comm%y_old)
end if
if (associated(comm%err_estimates)) then
   deallocate(comm%err_estimates); nullify(comm%err_estimates)
   nullify(comm%v0)
end if
if (associated(comm%p,comm%stages(:,1:2))) then
   nullify(comm%p)
end if
if (associated(comm%ge_stages,comm%stages)) then
   deallocate(comm%stages); nullify(comm%stages); nullify(comm%ge_stages)!zgs;
   nullify(comm%v1,comm%v2,comm%v3)
else if (associated(comm%ge_stages)) then
   deallocate(comm%ge_stages); nullify(comm%ge_stages)
end if
if (associated(comm%ge_y_new)) then
   deallocate(comm%ge_y_new); nullify(comm%ge_y_new)
end if
if (associated(comm%ge_assess)) then
   deallocate(comm%ge_assess); nullify(comm%ge_assess)
end if
if (associated(comm%ge_err_estimates)) then
   deallocate(comm%ge_err_estimates); nullify(comm%ge_err_estimates)
end if
if (associated(comm%ge_yp)) then
   deallocate(comm%ge_yp); nullify(comm%ge_yp)
end if
if (associated(comm%ge_y)) then
   deallocate(comm%ge_y); nullify(comm%ge_y)
end if
if (associated(comm%xstage)) then
   deallocate(comm%xstage); nullify(comm%xstage)
end if
if (associated(comm%p)) then
   deallocate(comm%p); nullify(comm%p)
end if
!
end subroutine collect_garbage_r1

recursive subroutine range_integrate_r1(comm,f,t_want,t_got,y_got,yderiv_got, &
                                        flag)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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), optional :: 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_want=y_got,yderiv_want=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_now=y_got,yderiv_now=yderiv_got, & 
                          flag=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

recursive subroutine step_integrate_r1(comm,f,t_now,y_now,yderiv_now,flag)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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), optional :: 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, phase2, 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
      phase2 = comm%phase2
      call step_r1(comm,f,comm%t,comm%y,comm%yp,comm%tol,htry, &
                   err,hmin,comm%phase2)
      comm%phase2 = 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 !zgs;
   if (present(y_now))y_now = comm%y
   if (present(yderiv_now)) 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


recursive subroutine truerr_r1(comm,f,ier)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
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 :: 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_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) !zgs
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) !zgs
   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_test1,hsec, &
         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 !zgs
      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 ! zgs; 
   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 !zgs; 
      ge_assess = ge_assess + (abs(ge_y-y)/weights)**2 
      ier = just_fine
   end if
   exit body
!
end do body
!
end subroutine truerr_r1

recursive subroutine step_r1(comm,f,tnow,y,yp,tol,htry,    &
                     err,hmin,phase_2)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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!
!
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 :: errest, y_new                !dep!
real(kind=wp), dimension(:,:), pointer :: stages                     !dep!
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.
!
main = present(hmin) .and. present(phase_2)
!
weights => comm%weights
thresh => comm%thresh
a => comm%a
b => comm%b
bhat => comm%bhat
c => comm%c
ptr => comm%ptr
!
if (main) then
   stages => comm%stages
   errest => comm%err_estimates
   y_new => comm%y_new
else
   stages => comm%ge_stages
   errest => comm%ge_err_estimates
   y_new => comm%ge_y_new
end if
!
attempt_step: do
!
   if (main) then
      if (phase_2) 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); phase_2 = .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 = argdif > comm%round_off*max(ynrm,ystgnm)
   if (cutbak) then
      cutbak = (abs(htry)*fdiff) > (comm%stability_radius*argdif)
      if (cutbak) then
         scl = max(comm%rs4,min(comm%rs1, &
                   (comm%stability_radius*argdif)/(abs(htry)*fdiff)))
         htry = scl*htry
      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


recursive subroutine stiff_r1(comm,f,toomch,sure_stiff)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
logical, intent(in) :: toomch
logical, intent(out) :: sure_stiff
!
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
!
logical :: maybe_stiff, lots_of_fails
!
real(kind=wp) :: alpha1, alpha2, beta1, beta2                        !dep!
real(kind=wp) :: rold, v1v0, v2v0, v2v1, v3v1, v3v2                  !dep!
real(kind=wp) :: dist, res2, scale, v0nrm, v3nrm, ynrm, rho, v0v0, v1v1, &
   v2v2, v3v3, yy, det1, det2
integer :: ntry
complex(kind=wp), dimension(2) :: root_pair, prv_root_pair
integer, parameter :: bigr=1, smlr=2
real(kind=wp), dimension(:), pointer :: v0, v1, v2, v3, y, y_old     !dep!
real(kind=wp), dimension(:), pointer :: weights, thresh              !shp-dep!
!
integer, parameter :: max_f_count=5000
real(kind=wp),parameter :: zero=0.0_wp, pt001=0.001_wp, pt9=0.9_wp, &
   fifth=0.2_wp, half=0.5_wp, one=1.0_wp, two=2.0_wp, five=5.0_wp, &
   large=1.0e+10_wp
!
v0 => comm%v0
v1 => comm%v1
v2 => comm%v2
v3 => comm%v3
weights => comm%weights
thresh => comm%thresh
y => comm%y
y_old => comm%y_old
!
sure_stiff = .false.
lots_of_fails = .false.
!
if (mod(comm%step_count-10,40)==0) then
   lots_of_fails = comm%stiff_bad_step_count >= 10
   comm%stiff_bad_step_count = 0
end if
!
!  If either too much work has been done or there are lots of failed steps,
!  test for stiffness.
!
maybe_stiff = toomch .or. lots_of_fails
if (maybe_stiff) then
!
!  Regenerate weight vector
!
   weights = max(half*(abs(y)+abs(y_old)),thresh)
   maybe_stiff = fifth <  abs(comm%h/comm%h_average) .and. &
                          abs(comm%h/comm%h_average) < five 
   if (maybe_stiff) then
!
!  The average step size is used to predict the cost in function evaluations
!  of finishing the integration to T_END.  If this cost is no more than 
!  MAX_F_COUNT, the problem is declared not stiff. If the step size is 
!  being restricted on grounds of stability, it will stay close to H_AVERAGE.
!  The prediction will then be good, but the cost is too low to consider
!  the problem stiff.  If the step size is not close to H_AVERAGE, the
!  problem is not stiff.  Either way there is no point to testing for a step
!  size restriction due to stability.
!
      maybe_stiff  = comm%cost*abs((comm%t_end-comm%t)/comm%h_average) > &
                     real(max_f_count,kind=wp) !zgs
      if (maybe_stiff) then
!
!  There have been many step failures or a lot of work has been done.  Now 
!  we must determine if this is due to the stability characteristics of the
!  formula.  This is done by calculating the dominant eigenvalues of the
!  local Jacobian and then testing whether H_AVERAGE corresponds to being
!  on the boundary of the stability region.
!  The size of Y provides scale information needed to approximate
!  the Jacobian by differences.
!
         v0v0 = wt_inner_prod(v0,v0)
         yy = wt_inner_prod(y,y)
         ynrm = sqrt(yy)
         scale = ynrm*comm%sqrrmc
         if (scale==zero) then
!
!  Degenerate case.  Y is (almost) the zero vector so the scale is not 
!  defined.  The input vector V0 is the difference between Y and a 
!  lower order approximation to the solution that is within the error 
!  tolerance.  When Y vanishes, V0 is itself an acceptable approximate
!  solution, so we take SCALE from it, if this is possible.
!
            scale = v0v0*comm%sqrrmc
            maybe_stiff = scale > zero
         end if
      end if
   end if
end if
!
if (.not. maybe_stiff) return
!
if (v0v0==zero) then
!
!  Degenerate case.  V0 is (almost) the zero vector so cannot
!  be used to define a direction for an increment to Y.  Try a
!  "random" direction.
!
   v0 = one;   v0v0 = wt_inner_prod(v0,v0)
end if
!
v0nrm = sqrt(v0v0)
v0 = v0/v0nrm; v0v0 = one
!
!  Use a nonlinear power method to estimate the two dominant eigenvalues.
!  V0 is often very rich in the two associated eigenvectors.  For this 
!  reason the computation is organized with the expectation that a minimal 
!  number of iterations will suffice.  Indeed, it is necessary to recognize 
!  a kind of degeneracy when there is a dominant eigenvalue.  The function
!  DOMINANT_EIGENVALUE does this.  In the first try, NTRY = 1, a Rayleigh 
!  quotient for such an eigenvalue is initialized as ROLD.  After each 
!  iteration, DOMINANT_EIGENVALUE computes a new Rayleigh quotient and
!  tests whether the two approximations agree to one tenth of one per cent
!  and the eigenvalue, eigenvector pair satisfy a stringent test on the 
!  residual.
!
ntry = 1
do
!
   v1 = approx_jacobian(f,v0,v0v0)
   v1v1 = wt_inner_prod(v1,v1)
!
!  The quantity SQRT(V1V1/V0V0) is a lower bound for the product of H_AVERAGE
!  and a Lipschitz constant.  If it should be LARGE, stiffness is not
!  restricting the step size to the stability region.  The principle is
!  clear enough, but the real reason for this test is to recognize an
!  extremely inaccurate computation of V1V1 due to finite precision
!  arithmetic in certain degenerate circumstances.
!
   if (sqrt(v1v1) > large*sqrt(v0v0)) return
!
   v1v0 = wt_inner_prod(v1,v0)
   if (ntry==1) then
      rold = v1v0/v0v0
!
!  This is the first Rayleigh quotient approximating the product of H_AVERAGE
!  and a dominant eigenvalue.  If it should be very small, the
!  problem is not stiff.  It is important to test for this possibility so
!  as to prevent underflow and degeneracies in the subsequent iteration.
!
      if (abs(rold) < comm%cubrmc) return
   else
!
      if (dominant_eigenvalue(v1v1,v1v0,v0v0)) exit
   end if
!
   v2 = approx_jacobian(f,v1,v1v1)
   v2v2 = wt_inner_prod(v2,v2)
   v2v0 = wt_inner_prod(v2,v0)
   v2v1 = wt_inner_prod(v2,v1)
   if (dominant_eigenvalue(v2v2,v2v1,v1v1)) exit
!
!  Fit a quadratic in the eigenvalue to the three successive iterates
!  V0,V1,V2 of the power method to get a first approximation to
!  a pair of eigenvalues.  A test made earlier in DOMINANT_EIGENVALUE
!  implies that the quantity DET1 here will not be too small.
!
   det1 = v0v0*v1v1 - v1v0*rev_wt_inner_prod(v1v0)
   alpha1 = (-v0v0*v2v1 + rev_wt_inner_prod(v1v0)*v2v0)/det1
   beta1 = (v1v0*v2v1 - v1v1*v2v0)/det1
!
!  Iterate again to get V3, test again for degeneracy, and then fit a
!  quadratic to V1,V2,V3 to get a second approximation to a pair
!  of eigenvalues.
!
   v3 = approx_jacobian(f,v2,v2v2)
   v3v3 = wt_inner_prod(v3,v3)
   v3v1 = wt_inner_prod(v3,v1)
   v3v2 = wt_inner_prod(v3,v2)
   if (dominant_eigenvalue(v3v3,v3v2,v2v2)) exit
!
   det2 = v1v1*v2v2 - v2v1*rev_wt_inner_prod(v2v1)
   alpha2 = (-v1v1*v3v2 + rev_wt_inner_prod(v2v1)*v3v1)/det2
   beta2 = (v2v1*v3v2 - v2v2*v3v1)/det2
!
!  First test the residual of the quadratic fit to see if we might
!  have determined a pair of eigenvalues.
!
   res2 = abs( v3v3 + rev_wt_inner_prod(alpha2)*v3v2 + &
               rev_wt_inner_prod(beta2)*v3v1 + &
               alpha2*rev_wt_inner_prod(v3v2) + &
               alpha2*rev_wt_inner_prod(alpha2)*v2v2 + &
               alpha2*rev_wt_inner_prod(beta2)*v2v1 + &
               beta2*rev_wt_inner_prod(v3v1) + &
               beta2*rev_wt_inner_prod(alpha2)*rev_wt_inner_prod(v2v1) + &
               beta2*rev_wt_inner_prod(beta2)*v1v1 )
   if (res2 <= abs(v3v3)*pt001**2) then
!
!  Calculate the two approximate pairs of eigenvalues.
!
      prv_root_pair(1:2) = quadratic_roots(alpha1,beta1)
      root_pair(1:2) = quadratic_roots(alpha2,beta2)
!
!  The test for convergence is done on the larger root of the second
!  approximation.  It is complicated by the fact that one pair of roots 
!  might be real and the other complex.  First calculate the spectral 
!  radius RHO of HAVG*J as the magnitude of ROOT1.  Then see if one of 
!  the roots R1,R2 is within one per cent of ROOT1.  A subdominant root 
!  may be very poorly approximated if its magnitude is much smaller than 
!  RHO -- this does not matter in our use of these eigenvalues.
!
      rho = abs(prv_root_pair(bigr))
      dist = min( abs(root_pair(bigr) - prv_root_pair(bigr)), &
                  abs(root_pair(bigr) - prv_root_pair(smlr)) )
      if (dist <= pt001*rho) exit
   end if
!
!  Do not have convergence yet.  Because the iterations are cheap, and
!  because the convergence criterion is stringent, we are willing to try
!  a few iterations.
!
   ntry = ntry + 1
   if (ntry > comm%max_stiff_iters) return
   v3nrm = sqrt(v3v3)
   v0 = v3/v3nrm
   v0v0 = one
!
end do
!
!  We now have the dominant eigenvalues.  Decide if the average step
!  size is being restricted on grounds of stability.  Check the real
!  parts of the eigenvalues.  First see if the dominant eigenvalue is
!  in the left half plane -- there won't be a stability restriction
!  unless it is. If there is another eigenvalue of comparable magnitude
!  with a positive real part, the problem is not stiff. If the dominant
!  eigenvalue is too close to the imaginary axis, we cannot diagnose
!  stiffness.
!
if ( real(root_pair(bigr)) < zero) then
   if ( .not. ( abs(root_pair(smlr)) >= pt9*rho .and. &
                real(root_pair(smlr)) > zero ) ) then
      if ( abs(aimag(root_pair(bigr))) <= &
           abs(real(root_pair(bigr)))*comm%tan_angle) then
!
!  If the average step size corresponds to being well within the
!  stability region, the step size is not being restricted because
!  of stability.
!
         sure_stiff = rho >= pt9*comm%stability_radius
      end if
   end if
end if
!
contains

! All non-procedure arguments to these contained procedures
! should have intent(in) but IBM RS/600 AIX xlf v3.2.1.x
! produces rubbish so intent attribute has been dropped

function approx_jacobian(f,v,vdotv)
!
real(kind=wp) :: vdotv
real(kind=wp), dimension(:) :: v                         !dep!
real(kind=wp), dimension(size(v,1)) :: approx_jacobian               !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) :: temp1
!
!  Scale V so that it can be used as an increment to Y
!  for an accurate difference approximation to the Jacobian.
!
temp1 = scale/sqrt(vdotv)
comm%vtemp = y + temp1*v
!
approx_jacobian = f(comm%t,comm%vtemp)
comm%f_count = comm%f_count + 1
!
!  Form the difference approximation.  At the same time undo
!  the scaling of V and introduce the factor of H_AVERAGE.
!
approx_jacobian = &
    (comm%h_average/temp1)*(approx_jacobian-comm%yp)
!
end function approx_jacobian


function quadratic_roots(alpha,beta)
!
real(kind=wp) :: alpha, beta                             !dep!
complex(kind=wp), dimension(2) :: quadratic_roots
!
complex(kind=wp) :: temp, sqdisc, r1, r2
!
!  For types other than real/complex, this procedure must be constructed
!  such that temp and sqdisc are evaluated as compelx quantities
!
temp = alpha/two; sqdisc = sqrt(temp**2 - beta)
!
! Do we have double root?
!
if (sqdisc==zero) then
   quadratic_roots = (/ -temp, -temp /)
!
! Distinct roots
!
else
   r1 = -temp + sqdisc; r2 = -temp + sqdisc
   if (abs(r1) > abs(r2)) then
      quadratic_roots = (/ r1, r2 /)
   else
      quadratic_roots = (/ r2, r1 /)
   end if
end if
!
end function quadratic_roots

function dominant_eigenvalue(v1v1,v1v0,v0v0)
!
real(kind=wp) :: v0v0, v1v1
real(kind=wp) :: v1v0                                    !dep!
logical :: dominant_eigenvalue
!
real(kind=wp) :: ratio                                               !dep!
real(kind=wp) :: res, det
logical :: big
!
ratio = v1v0/v0v0; rho = abs(ratio)
det = v0v0*v1v1 - v1v0*rev_wt_inner_prod(v1v0); res = abs(det/v0v0)
!
big = det == zero .or. &
                  (res<=abs(v1v1)*pt001**2 .and. abs(ratio-rold)<=pt001*rho)
!
if (big) then
   root_pair(bigr) = cmplx(ratio)
   root_pair(smlr) = cmplx(zero)
end if
!
rold = ratio
dominant_eigenvalue = big
!
end function dominant_eigenvalue

function wt_inner_prod(vec_1,vec_2)
!
real(kind=wp), dimension(:) :: vec_1, vec_2              !dep!
real(kind=wp) :: wt_inner_prod                                       !dep!
!
!
wt_inner_prod = sum ( (vec_1/weights) * (vec_2/weights) )            !spec-ar!
!
end function wt_inner_prod

function rev_wt_inner_prod(value)
!
real(kind=wp) :: value                                   !dep!
real(kind=wp) :: rev_wt_inner_prod                                   !dep!
!
! given result of inner product value = v1.v0
! must return the reverse, ie v0.v1
!
! for real variables the value is the same
! for complex need to conjugate
!
rev_wt_inner_prod = value                                            !spec-line!
!
end function rev_wt_inner_prod

end subroutine stiff_r1

subroutine statistics_r1(comm,total_f_calls,step_cost,waste,num_succ_steps,&
                         h_next,y_maxvals)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
real(kind=wp), optional, intent(out) :: h_next                       !indep!
real(kind=wp), optional, intent(out) :: waste
real(kind=wp), dimension(:), optional, intent(out) :: y_maxvals      !shp-dep!
integer, optional, intent(out) :: step_cost, num_succ_steps, total_f_calls
!
character(len=*), parameter :: srname="STATISTICS"
!
integer :: ier, nrec, state
!
integer, parameter :: not_ready=-1, not_reusable=-3, fatal=911, &
   catastrophe=912, just_fine=1
logical, parameter :: ask=.true.
real(kind=wp), parameter :: zero=0.0_wp
!
ier = just_fine; nrec = 0
!
body: do
!
!  Is it permissible to call STATISTICS?
!
   state = get_saved_state_r1(srname,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_reusable) then
      ier = fatal; nrec = 2; write (comm%rec,"(a/a)") &
" ** You have already made a call to STATISTICS after a hard failure was ", &
" ** reported from the integrator. You cannot call STATISTICS again."
      exit body
   end if
   state = get_saved_state_r1("STEP_INTEGRATE",comm%save_states)
   if (state==not_ready) then
      ier = fatal; nrec = 1
      if (comm%use_range) then
         write (comm%rec,"(a)") &
" ** You have not called RANGE_INTEGRATE, so you cannot use STATISTICS."
      else
         write (comm%rec,"(a)") &
" ** You have not called STEP_INTEGRATE, so you cannot use STATISTICS."
      end if
      exit body
   end if
   if (present(y_maxvals)) then
      if (any(shape(y_maxvals) /= shape(comm%y))) then
         ier = fatal; nrec = 2; write (comm%rec,"(a,i6,a/a,i6,a)") &
" ** The shape of Y_MAXVALS is not consistent with the shape of the", &
" ** dependent variables."
         exit body
      end if
   end if
!
!  Set flag so that the routine can only be called once after a hard 
!  failure from the integrator.
!
   if (state==5 .or. state==6) ier = not_reusable
!
   if (present(total_f_calls)) then
      total_f_calls = comm%full_f_count + comm%f_count
!      if (comm%erason) total_f_calls = total_f_calls + comm%ge_f_count
   end if
   if (present(step_cost)) step_cost = comm%cost
   if (present(num_succ_steps)) num_succ_steps = comm%step_count
   if (present(waste)) then
      if (comm%step_count<=1) then
         waste = zero
      else
         waste = real(comm%bad_step_count,kind=wp) / & !zgs
                 real(comm%bad_step_count+comm%step_count,kind=wp) !zgs
      end if
   end if
   if (present(h_next)) h_next = comm%h
   if (present(y_maxvals)) y_maxvals = comm%ymax
   exit body
end do body
!
call rkmsg_r1(ier,srname,nrec,comm)
!
end subroutine statistics_r1

subroutine global_error_r1(comm,rms_error,max_error,t_max_error)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
real(kind=wp), optional, intent(out) :: max_error
real(kind=wp), optional, intent(out) :: t_max_error                  !indep!
real(kind=wp), dimension(:), optional, intent(out) :: rms_error      !shp-dep!
!
character(len=*), parameter :: srname="GLOBAL_ERROR"
!
integer :: ier, nrec, state
!
intrinsic         sqrt
!
integer, parameter :: not_ready=-1, not_reusable=-3, fatal=911, &
   catastrophe=912, just_fine=1
logical, parameter :: ask=.true.
!
ier = just_fine; nrec = 0
!
body: do 
!
!  Is it permissible to call GLOBAL_ERROR?
!
   state = get_saved_state_r1(srname,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_reusable) then
      ier = fatal; nrec = 2; write (comm%rec,"(a/a)") &
" ** You have already made a call to GLOBAL_ERROR after a hard failure was", &
" ** reported from the integrator. You cannot call GLOBAL_ERROR again."
      exit body
   end if
   state = get_saved_state_r1("STEP_INTEGRATE",comm%save_states)
   if (state==not_ready) then
      ier = fatal; nrec = 1
      if (comm%use_range) then
         write (comm%rec,"(a)") &
" ** You have not yet called RANGE_INTEGRATE, so you cannot call GLOBAL_ERROR."
      else
         write (comm%rec,"(a)") &
" ** You have not yet called STEP_INTEGRATE, so you cannot call GLOBAL_ERROR."
      end if
      exit body
   end if
!
!  Set flag so that the routine can only be called once after a hard 
!  failure from the integrator.
!
   if (state==5 .or. state==6) ier = not_reusable
!
!  Check that ERROR_ASSESS was set properly for error assessment in SETUP.
!
   if (.not.comm%erason) then
      ier = fatal; nrec = 3; write (comm%rec,"(a/a/a)") &
" ** No error assessment is available since you did not ask for it in your",&
" ** call to the routine SETUP. Check your program carefully."
      exit body
   end if
!
! Check size of RMS_ERROR
!
   if (present(rms_error)) then
      if (any(shape(rms_error) /= shape(comm%y))) then
         ier = fatal; nrec = 2; write (comm%rec,"(a,a)") &
" ** The shape of RMS_ERROR is not consistent with the shape of the", &
" ** dependent variables."
         exit body
      end if
   end if
!
!  Check to see if the integrator has not actually taken a step.
!
   if (comm%step_count==0) then
      ier = fatal; nrec = 2; write (comm%rec,"(a/a)") &
" ** The integrator has not actually taken any successful steps. You cannot",&
" ** call GLOBAL_ERROR in this circumstance. Check your program carefully."
       exit body
   end if
!
!  Compute RMS error and set output variables.
!
   if (present(max_error)) max_error = comm%ge_max_contrib
   if (present(t_max_error)) t_max_error = comm%t_ge_max_contrib
   if (present(rms_error)) rms_error = &
      sqrt(comm%ge_assess/real(comm%step_count,kind=wp)) !zgs
!
   exit body
end do body
!
call rkmsg_r1(ier,srname,nrec,comm)
!
end subroutine global_error_r1

subroutine reset_t_end_r1(comm,t_end_new)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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_end_new                               !indep!
type(rk_comm_real_1d), intent(inout) :: comm
!
character(len=*), parameter :: srname="RESET_T_END"
!
real(kind=wp) :: hmin, tdiff                                         !indep!
integer ::           ier, nrec, state
!
integer, parameter :: not_ready=-1, usable=-2, fatal=911, catastrophe=912, &
   just_fine=1
logical, parameter :: ask=.true.
real(kind=wp), parameter :: zero=0.0_wp
!
ier = just_fine; nrec = 0
!
!  Is it permissible to call RESET_T_END?
!
body: do
!
   state = get_saved_state_r1("STEP_INTEGRATE",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 (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 RESET_T_END after you specified to SETUP that you were",&
" ** going to use RANGE_INTEGRATE. This is not permitted."
         exit body
      end if
   end if
   if (state==not_ready) then
      ier = fatal; nrec = 1; write (comm%rec,"(a)") &
" ** You have not called STEP_INTEGRATE, so you cannot use RESET_T_END."
      exit body
   end if
   if (state==5 .or. state==6) then
      ier = fatal; nrec = 2; write (comm%rec,"(a,i1,a/a)") &
" ** STEP_INTEGRATE has returned with FLAG =  ",STATE," You cannot call",&
" ** RESET_T_END inthis circumstance."
      exit body
   end if
!
!  Check value of T_END_NEW
!
   if (comm%dir>zero .and. t_end_new<=comm%t) then
      ier = fatal; nrec = 3; write (comm%rec,"(a/a,e13.5/a,e13.5,a)") &
" ** Integration is proceeding in the positive direction. The current value",&
" ** for the independent variable is ",comm%T," and you have set T_END_NEW =",&
" ** ",T_END_NEW,".  T_END_NEW must be greater than T."
      exit body
   else if (comm%dir<zero .and. t_end_new>=comm%t) then
      ier = fatal; nrec = 3; write (comm%rec,"(a/a,e13.5/a,e13.5,a)") &
" ** Integration is proceeding in the negative direction. The current value",&
" ** for the independent variable is ",comm%T," and you have set T_END_NEW =",&
" ** ",T_END_NEW,".  T_END_NEW must be less than T."
      exit body
   else
      hmin = max(comm%sqtiny,comm%toosml*max(abs(comm%t),abs(t_end_new)))
      tdiff = abs(t_end_new-comm%t)
      if (tdiff<hmin) then
         ier = fatal; nrec = 4 
         write (comm%rec,"(a,e13.5,a/a,e13.5,a/a/a,e13.5,a)")&
" ** The current value of the independent variable T is ",comm%T,". The",&
" ** T_END_NEW you supplied has ABS(T_END_NEW-T) = ",TDIFF,". For the METHOD",&
" ** and the precision of the computer being used, this difference must be",&
" ** at least ",HMIN,"."
         exit body
      end if
   end if
!
   comm%t_end = t_end_new; comm%at_t_end = .false.
!
   exit body
end do body
!
call rkmsg_r1(ier,srname,nrec,comm)
!
end subroutine reset_t_end_r1

recursive subroutine interpolate_r1(comm,f,t_want,y_want,yderiv_want)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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!
type(rk_comm_real_1d), intent(inout), target :: comm
real(kind=wp), dimension(:), intent(out), optional :: y_want         !dep!
real(kind=wp), dimension(:), intent(out), optional :: yderiv_want    !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="INTERPOLATE"
integer :: ier, jer, nrec, state, npcls
logical :: intrp_initialised
!
integer, parameter :: not_ready=-1, usable=-2, fatal=911, catastrophe=912, &
   just_fine=1
logical, parameter :: ask=.true.
!
ier = just_fine; nrec = 0
!
body: do
!
!  Is it permissible to call INTERPOLATE?
!
   state = get_saved_state_r1("STEP_INTEGRATE",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 (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 INTERPOLATE after you specified to SETUP that you were",&
" ** going to use RANGE_INTEGRATE. This is not permitted."
         exit body
      end if
   end if
   if (state==not_ready) then
      ier = fatal; nrec = 1; write (comm%rec,"(a)") &
" ** You have not called STEP_INTEGRATE, so you cannot use INTERPOLATE."
      exit body
   end if
   if (state > just_fine) then
      ier = fatal; nrec = 2; write (comm%rec,"(a/a)") &
" ** STEP_INTEGRATE has returned with a flag value greater than 1. You", &
" ** cannot call INTERPOLATE in this circumstance."
      exit body
   end if
!
!  Check sizes of arrays
!
   if (present(y_want)) then
      if (any(shape(y_want)/=shape(comm%y))) then
         ier = fatal; nrec = 3; write (comm%rec,"(a,i6,a/a,i6,a/a)") &
" ** The shape of the array Y_WANT is not consistent with the shape of the ", &
" ** dependent variables."
         exit body
      end if
   end if
   if (present(yderiv_want)) then
      if (any(shape(yderiv_want)/=shape(comm%y))) then
         ier = fatal; nrec = 3; write (comm%rec,"(a,i6,a/a,i6,a/a)") &
" ** The shape of the array YDERIV_WANT is not consistent with the shape of", &
" ** the dependent variables."
         exit body
      end if
   end if
!
!  Check METHOD is ok to interpolate with
!
   if (comm%rk_method==3) then
      ier = fatal; nrec = 5; write (comm%rec,"(a/a/a/a/a)") &
" ** You have been using STEP_INTEGRATE with METHOD = 'H' to integrate your",&
" ** equations. You have just called INTERPOLATE, but interpolation is not",&
" ** available for this METHOD. Either use METHOD = 'M', for which",&
" ** interpolation is available, or use RESET_T_END to make STEP_INTEGRATE",&
" ** step exactly to the points where you want output."
      exit body
   end if
!
!  Get some workspace -
!     can overwrite STAGES in METHOD 'L' since they're not requird for the
!     interpolant
!
   select case(comm%rk_method)
   case(1)
      if (.not. associated(comm%p)) comm%p => comm%stages(:,1:2)
      npcls = 2
      if (.not. associated(comm%ytemp)) comm%p => comm%stages(:,1:3)
   case(2)
      jer = 0
      if (.not.associated(comm%xstage)) then
         allocate(comm%xstage(size(comm%y,1)),stat=jer)              !alloc!
      end if
      if (.not.associated(comm%ytemp)) then
         allocate(comm%ytemp(size(comm%y,1)),stat=jer)               !alloc!
      end if
      if (.not.associated(comm%p)) then
         allocate(comm%p(size(comm%y,1),5),stat=jer)                 !alloc!
      end if
      npcls = 5
      if (jer /= 0) then
         ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
         exit body
      end if
   end select
!
!  Check data to see if interpolant has already been calculated for this
!  step
!
   intrp_initialised = get_saved_state_r1(srname,comm%save_states) /= usable
!
!  Some initialization must be done before interpolation is possible.
!
   if (.not.intrp_initialised) call form_intrp(f,comm%p)
!
!  The actual evaluation of the interpolating polynomial and/or its first
!  derivative is done in EVALUATE_INTRP.
!
   call evaluate_intrp(comm%p,y_want,yderiv_want)
   exit body
!
end do body
!
call rkmsg_r1(ier,srname,nrec,comm)
!
contains
!
recursive subroutine form_intrp(f,p)
!
real(kind=wp), intent(out), dimension(:,:) :: p                      !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), dimension(:,:), pointer :: r                          !real!
real(kind=wp), dimension(:,:), pointer :: stages                     !dep!
real(kind=wp), dimension(:), pointer :: y, yp, y_old, yp_old         !dep!
real(kind=wp), dimension(:), pointer :: xstage                       !dep!
!
stages => comm%stages
r => comm%r
y => comm%y
yp => comm%yp
y_old => comm%y_old
yp_old => comm%yp_old
xstage => comm%xstage
!
select case(comm%rk_method)
case(1)
!
!  METHOD = 'L'.  Use the cubic Hermite interpolant that is fully
!  specified by the values and slopes at the two ends of the step.
!
   p(:,2) = y - y_old
   p(:,1) = comm%h_old*yp - p(:,2)
   p(:,2) = p(:,1) - (p(:,2)-comm%h_old*yp_old)
   p(:,1) = p(:,1) + p(:,2)
!
case(2)
!
!  METHOD = 'M'.
!
   if (.not.intrp_initialised) call extra_stages(f,comm%ytemp,xstage)
!
!  Form the coefficients of the interpolating polynomial in its shifted
!  and scaled form.  The transformation from the form in which the
!  polynomial is derived can be somewhat ill-conditioned.  The terms 
!  are grouped so as to minimize the errors of the transformation.
!
!  Coefficient of SIGMA**6
   p(:,5) = r(5,6)*stages(:,4) + &
                     ((r(10,6)*xstage+r(8,6)*yp)+ &
                     (r(7,6)*stages(:,6)+r(6,6)*stages(:,5))) + &
                     ((r(4,6)*stages(:,3)+r(9,6)*stages(:,7))+ &
                     (r(3,6)*stages(:,2)+r(11,6)*stages(:,1))+ &
                     r(1,6)*yp_old)
!
!  Coefficient of SIGMA**5
   p(:,4) = (r(10,5)*xstage+r(9,5)*stages(:,7)) + &
                     ((r(7,5)*stages(:,6)+r(6,5)*stages(:,5))+ &
                     r(5,5)*stages(:,4)) + ((r(4,5)*stages(:,3)+ &
                     r(8,5)*yp)+(r(3,5)*stages(:,2)+r(11,5)* &
                     stages(:,1))+r(1,5)*yp_old)
!
!  Coefficient of SIGMA**4
   p(:,3) = ((r(4,4)*stages(:,3)+r(8,4)*yp)+ &
                     (r(7,4)*stages(:,6)+r(6,4)*stages(:,5))+ &
                     r(5,4)*stages(:,4)) + ((r(10,4)*xstage+ &
                     r(9,4)*stages(:,7))+(r(3,4)*stages(:,2)+ &
                     r(11,4)*stages(:,1))+r(1,4)*yp_old)
!
!  Coefficient of SIGMA**3
   p(:,2) = r(5,3)*stages(:,4) + r(6,3)*stages(:,5) + &
                     ((r(3,3)*stages(:,2)+r(9,3)*stages(:,7))+ &
                     (r(10,3)*xstage+r(8,3)*yp)+r(1,3)* &
                     yp_old)+((r(4,3)*stages(:,3)+r(11,3)* &
                     stages(:,1))+r(7,3)*stages(:,6))
!
!  Coefficient of SIGMA**2
   p(:,1) = r(5,2)*stages(:,4) + ((r(6,2)*stages(:,5)+ &
                   r(8,2)*yp)+r(1,2)*yp_old) + &
                   ((r(3,2)*stages(:,2)+r(9,2)*stages(:,7))+ &
                   r(10,2)*xstage) + ((r(4,2)*stages(:,3)+ &
                   r(11,2)*stages(:,1))+r(7,2)*stages(:,6))
!
!  Scale all the coefficients by the step size.
   p(:,:) = comm%h_old*p(:,:)
!
end select
!
end subroutine form_intrp


subroutine evaluate_intrp(p,y_want,yderiv_want)
!
real(kind=wp), dimension(:), optional, intent(out) :: y_want         !dep!
real(kind=wp), dimension(:), optional, intent(out) :: yderiv_want    !dep!
real(kind=wp), dimension(:,:), intent(in) :: p                       !dep!
!
real :: sigma
integer :: i
!
sigma = (t_want-comm%t)/comm%h_old
!
if (present(y_want)) then
   y_want = p(:,comm%intrp_degree-1)*sigma
   do i = comm%intrp_degree - 2, 1, -1
      y_want = (y_want+p(:,i))*sigma
   end do
   y_want = (y_want+comm%h_old*comm%yp)*sigma + comm%y
end if
!
if (present(yderiv_want)) then
   yderiv_want = comm%intrp_degree*p(:,comm%intrp_degree-1)*sigma
   do i = comm%intrp_degree - 1, 2, -1
      yderiv_want = (yderiv_want+i*p(:,i-1))*sigma
   end do
   yderiv_want = (yderiv_want+comm%h_old*comm%yp)/comm%h_old
end if
!
end subroutine evaluate_intrp


recursive subroutine extra_stages(f,ytemp,xstage)
!
real(kind=wp), dimension(:), intent(out) :: ytemp, xstage            !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), dimension(:,:), pointer :: stages                     !dep!
real(kind=wp), dimension(:), pointer :: yp, y_old, yp_old            !dep!
!
real(kind=wp), dimension(:,:), pointer :: a                          !real!
real(kind=wp), dimension(:), pointer :: c                            !real!
real(kind=wp), pointer :: h_old, t_old                               !indep!
!
integer :: i, j
!
a => comm%a
stages => comm%stages
c => comm%c
yp => comm%yp
y_old => comm%y_old
yp_old => comm%yp_old
h_old => comm%h_old
t_old => comm%t_old
!
!  Compute the extra stages needed for interpolation using the facts that
!       1. Stage 1 is YP_OLD.
!       2. Stage i (i>1) is stored in STAGES(...,i-1).
!       3. This pair is FSAL, i.e. STAGES(...,7)=YP, which frees
!          up STAGES(...,7) for use by stage 9.
!       4. XSTAGE is used for stage 10.
!       5. The coefficient of stage 2 in the interpolant is always 0, so
!          STAGES(...,1) is used for stage 11.
!
do i = 9, 11
   do j = 1, i-1
      select case (j)
         case(1); ytemp = a(i,1)*yp_old
! could have used matmul here but that prevents increasing rank of dep-var
         case(2:7);ytemp = ytemp + a(i,j)*stages(:,j-1)
         case(8); ytemp = ytemp + a(i,j)*yp
         case(9); ytemp = ytemp + a(i,j)*stages(:,7)
         case(10); ytemp = ytemp + a(i,j)*xstage
      end select
   end do
   ytemp = y_old + h_old*ytemp
   select case(i)
      case(9)
         stages(:,7) = f(t_old+c(i)*h_old,ytemp)
      case(10)
         xstage = f(t_old+c(i)*h_old,ytemp)
      case(11)
         stages(:,1) = f(t_old+c(i)*h_old,ytemp)
   end select
   comm%f_count = comm%f_count + 1
end do
!
end subroutine extra_stages

end subroutine interpolate_r1

subroutine rkmsg_r1(ier,srname,nrec,comm,flag)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
!
integer, intent(in) :: ier, nrec
integer, intent(out), optional :: flag
character(len=*), intent(in) :: srname
type(rk_comm_real_1d), intent(inout) :: comm
!
logical :: ok, on, range_call
!
integer, parameter :: fatal=911, catastrophe=912, just_fine=1
logical, parameter :: tell=.false.
!
!  Check where the call came from - if it is an indirect call from 
!  RANGE_INTEGRATE the run is not STOPped.
!
range_call = (srname=="RESET_T_END" .or. srname=="STEP_INTEGRATE" .or. &
          srname=="INTERPOLATE") .and. comm%use_range
!
!  Check if can continue with integrator.
!
ok = (srname=="STEP_INTEGRATE" .or. srname=="RANGE_INTEGRATE") .and. &
     (ier==2 .or. ier==3 .or. ier==4)
!
!  Check if program termination has been overridden.
!
on = get_stop_on_fatal_r1(comm)
!
if ((comm%print_message.and.ier>just_fine) .or. ier>=fatal) then
   write (comm%outch,"(/a)") " **"
   write (comm%outch,"(a)") comm%rec(1:nrec)
   if (ier>=fatal) then
      write (comm%outch,"(a/a,a,a/a/)") &
" **",&
" ** Catastrophic error detected in ", srname, ".",&
" **"
      if ((.not.range_call.and.on.and.ier==fatal) .or. ier==catastrophe) then
         write (comm%outch,"(a/a/a)") &
" **",&
" ** Execution of your program is being terminated.",&
" **"
         stop
      end if
   else 
      if (ok) then
         write (comm%outch,"(a/a,a,a,i2,a/a/a)")  &
" **", &
" ** Warning from routine ", srname, " with flag set ",ier, ".",&
" ** You can continue integrating this problem.",&
" **"
      else
         write (comm%outch,"(a/a,a,a,i2,a/a/a)")  &
" **", &
" ** Warning from routine ", srname, " with flag set ",ier, ".", &
" ** You cannot continue integrating this problem.", &
" **"
      end if
      if (.not.present(flag)) then
         write (comm%outch,"(a/a/a)") &
" **",&
" ** Execution of your program is being terminated.",&
" **"
         stop
      end if
   end if
end if
!
if (present(flag)) flag = ier
comm%rec(nrec+1:10) = " "
!
!  Save the status of the routine associated with SRNAME
!
call set_saved_state_r1(srname,ier,comm)
!
!  Indicate that a catastrophic error has been detected
!
!call set_saved_fatal_r1(comm,ier >= catastrophe)
!
end subroutine rkmsg_r1

subroutine set_saved_state_r1(srname,state,comm)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
!
integer, intent(in) :: state
type(rk_comm_real_1d), intent(inout) :: comm
character(len=*), intent(in) :: srname
!
integer :: name
!
integer, parameter :: fatal=911
!
select case (srname)
   case("SETUP"); name = 1
   case("RANGE_INTEGRATE"); name = 2
   case("STATISTICS"); name = 3
   case("GLOBAL_ERROR"); name = 4
   case("STEP_INTEGRATE"); name = 5
   case("INTERPOLATE"); name= 6
   case("RESET_T_END"); name = 7
   case default; name = 0
end select
!
comm%save_states(name) = state
comm%saved_fatal_err = state >= fatal
!
end subroutine set_saved_state_r1

function get_saved_state_r1(srname,save_states)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
!
integer, dimension(7), intent(inout) :: save_states
character(len=*), intent(in) :: srname
integer :: get_saved_state_r1
!
integer :: name
!
integer, parameter :: fatal=911
!
select case (srname)
   case("SETUP"); name = 1
   case("RANGE_INTEGRATE"); name = 2
   case("STATISTICS"); name = 3
   case("GLOBAL_ERROR"); name = 4
   case("STEP_INTEGRATE"); name = 5
   case("INTERPOLATE"); name= 6
   case("RESET_T_END"); name = 7
   case default; name = 0
end select
!
!  Check for status of given routine but check for any fatal errors first
!
if (any(save_states(1:7)==fatal)) then
   get_saved_state_r1 = fatal
else
   get_saved_state_r1 = save_states(name)
end if
!
end function get_saved_state_r1

function get_saved_fatal_r1(comm)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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(in) :: comm
logical :: get_saved_fatal_r1
!
get_saved_fatal_r1 = comm%saved_fatal_err
!
end function get_saved_fatal_r1

subroutine set_stop_on_fatal_r1(comm,action)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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
logical, intent(in) :: action
!
comm%stop_on_fatal = action
!
end subroutine set_stop_on_fatal_r1

function get_stop_on_fatal_r1(comm)
!
! Part of rksuite_90 v1.1 (May 1995)
!         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(in) :: comm
logical get_stop_on_fatal_r1
!
get_stop_on_fatal_r1 = comm%stop_on_fatal
!
end function get_stop_on_fatal_r1

!endc!
end module rksuite_90

! The following are examples

module invar_imbed_fw

use rksuite_90_prec


contains

function g(t,y)
real(kind=wp), intent(in) :: t
real(kind=wp), dimension(:), intent(in) :: y
real(kind=wp), dimension(size(y)) :: g

g(:) = (/ 1.0_wp - y(1)**2, -y(1)*y(2) /)

end function g

end module invar_imbed_fw

module invar_imbed_bw

use rksuite_90
use invar_imbed_fw

real(kind=wp) :: u, v

real(kind=wp) :: t_fw_start=0.0_wp, t_fw_end, t_fw_got, tol=1.0e-5_wp
real(kind=wp), dimension(2) :: thres(2)=1.0e-20_wp, &
    y_fw_start, y_fw_got
type(rk_comm_real_1d) :: comm_fw

contains

function f(t,y)

real(kind=wp), intent(in) :: t
real(kind=wp), intent(in) :: y
real(kind=wp) :: f

t_fw_end = t
u = 0.0_wp; v = 1.0_wp; y_fw_start = (/u, v/)

if (t_fw_end /= t_fw_start) then
   call setup(comm_fw,t_fw_start,y_fw_start,t_fw_end,tol,thres)
   call range_integrate(comm_fw,g,t_fw_end,t_fw_got,y_got=y_fw_got)
   call collect_garbage(comm_fw)
   u = y_fw_got(1); v = y_fw_got(2)
end if

f = u*y + v

end function f

end module invar_imbed_bw

module two_body_f

! PROBLEM
! 
! Integrate a two body problem.  The equations for the coordinates
! (x(t),y(t)) of one body as functions of time t in a suitable frame of
! reference are
! 
!     x'' = - x/r**3,
!     y'' = - y/r**3,   where   r = SQRT(x**2 + y**2).
! 
! The initial conditions lead to elliptic motion with eccentricity ECC.
! This parameter will be taken to be 0.9.
! 
!     x(0) = 1-ECC,     x'(0) = 0,
!     y(0) = 0,         y'(0) = SQRT((1+ECC)/(1-ECC)).
! 
! An accurate solution that shows the general behavior of the orbit is
! desired.  The coordinates will be returned at every time step in [0,20].
! This is a standard test problem for which there is a semi-analytical
! solution.  It will be compared to the computed solution at the end of
! the interval.
! 
! SOLUTION
! 
! For illustrative purposes we shall pose the system using a two dimensional
! array of unkowns
!        
!      ( x   x' )
!      ( y   y' )
! 
! to obtain
! 
!      ( x   x' )'  =  ( x'  -x/r**3 )
!      ( y   y' )      ( y'  -y/r**3 )
! 
! Since it is the general behavior of the solution that is desired, it is
! best to let the integrator choose where to provide answers.  It will
! produce answers more frequently where the solution changes more rapidly.
! Because the solution is inspected at every step, the task is to use
! STEP_INTEGRATE.
! 
! On physical grounds the solution is expected to be somewhat unstable
! when one body approaches the other. To obtain an accurate solution, a
! stringent relative error tolerance should be imposed -- TOL = 1.0D-10
! will be used.  At a tolerance this stringent the highest order pair,
! METHOD = 'H', is likely to be the most efficient choice. This method is
! inefficient when it is to produce answers at a great many specific
! points. It is most effective when used as in this template.  The
! solution components are expected to be of order 1, so threshold values
! of THRESHOLDS = 1.0D-13 are reasonable.  When a solution component is smaller
! in magnitude than this threshold, the code will control the local error
! to be no more than TOL*THRESHOLDS = 1.0D-23.  The reasonableness of this
! choice will be monitored by printing out the maximum value seen for each
! solution component in the course of the integration.  Error and warning
! messages, if any, will be printed out.
! 
! This is the standard test problem D5 of T.E. Hull, W.H. Enright,
! B.M. Fellen, and A.E. Sedgwick, "Comparing Numerical Methods
! for Ordinary Differential Equations," SIAM Journal on Numerical
! Analysis, Vol. 9, pp. 603-637, 1972.  The analytical solution
! in terms of the numerical solution of Kepler's equation can be
! found there as well as in most discussions of the two body
! problem.  The results for the particular choice of eccentricity,
! initial conditions, and interval  are providedin TRUE_Y_AT_TEND.
! 
! Global error assessment has been selected as a check on the reliability
! of the results. From the results generated it will be seen that the
! global error at the end of the run is about 2.3E-9, rather bigger than
! the local error tolerance TOL.  This illustrates the point that at best
! one can anticipate global errors comparable to the tolerance.  In point
! of fact, this problem is unstable at some points of the integration and
! the global error assessment reveals that the worst global error is
! considerably worse than the error at the end -- an example of the value
! of the global error assessment capability.

use rksuite_90_prec

contains

function f(t,y)
real(kind=wp), intent(in) :: t
real(kind=wp), dimension(:,:), intent(in) :: y
real(kind=wp), dimension(size(y,1),size(y,2)) :: f

real(kind=wp) :: r_cubed

r_cubed = sqrt( sum ( y(:,1)**2 ) )**3
f(:,1) =  y(:,2) !zgs;
f(:,2) = -y(:,1)/r_cubed

end function f

end module two_body_f

program two_body

use two_body_f
use rksuite_90

implicit none

integer :: flag, i, totf
real(kind=wp) :: t_start=0.0_wp, t_end=20.0_wp, tol=1.0e-10_wp, t_now, &
     ecc=0.9_wp, true_error, waste, max_error, t_max_error
real(kind=wp), dimension(2,2) :: y_start, thres = 1.0e-13_wp, y_now, &
     y_maxvals, true_y_at_t_end, assessed_error

type(rk_comm_real_1d) :: comm

true_y_at_t_end(:,1) = (/-1.29526625098758_wp,  0.400393896379232_wp /)
true_y_at_t_end(:,2) = (/-0.67753909247075_wp, -0.127083815427869_wp /)

y_start(:,1) = (/ 1.0_wp - ecc, 0.0_wp /)
y_start(:,2) = (/ 0.0_wp,  sqrt((1.0_wp+ecc)/(1.0_wp-ecc)) /)

! call setup(comm,t_start,y_start,t_end,tol,thres,method='high',task='step', &
!          error_assess=.true.)
! IBM F90 does not like it, lxm 5/9/96
do
!  call step_integrate(comm,f,t_now,y_now=y_now,flag=flag)
! IBM F90 does not like it, lxm 5/9/96
!  write (*,'(5e14.5') t_now, y_now
   if (t_now==t_end .or. flag > 3) exit
end do

! call statistics(comm,y_maxvals=y_maxvals,total_f_calls=totf,waste=waste)
! IBM F90 does not like it, lxm 5/9/96

write (*,'(/a,1p2e9.2/21x,1p2e9.2)') &
      '             ymax    ', (y_maxvals(i,1:2),i=1,2)
write (*,'(/a,1pe10.2,0p/a,i10/a,f10.2)') &
      ' The integration reached                      ', t_now, &
      ' The cost of the integration in calls to F was', totf, &
      ' The fraction of failed steps was             ', waste

if (t_now==t_end) then
   true_error = maxval(abs(true_y_at_t_end - y_now)/max(abs(y_now),thres))
!  call global_error(comm,rms_error=assessed_error,max_error=max_error, &
!                    t_max_error=t_max_error)
! IBM F90 does not like it, lxm 5/9/96
   write (*,'(/(a,1p,e9.2/a,e9.2/a,e9.2,a,e9.2))') &
      ' At t = 20, the true error is',true_error, &
      '            assessed error is',maxval(assessed_error), &
      '            maximum error was',max_error,' at ',t_max_error
end if

call collect_garbage(comm)

end program two_body
