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.0 (Aug 1994)
! software for initial value problems in ODEs
!
! Authors: R.W. Brankin (NAG Ltd., Oxford, England)
! I. Gladwell (Math Dept., SMU, Dallas, TX, USA)
! see main doc for contact details
!
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)
else
cost = real(no_of_stages,kind=wp)
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