1 SUBROUTINE stode(neq, y, yh, nyh, yh1, ewt, savf, acor, wm, iwm,
9 INTEGER,
DIMENSION(*) :: neq, iwm
10 REAL(rprec),
DIMENSION(*) :: y
11 REAL(rprec),
DIMENSION(nyh,*) :: yh
12 REAL(rprec),
DIMENSION(*) :: yh1, ewt, savf, acor, wm
17 common /ls0001/ conit, crate, el(13), elco(13,12), hold, rmax,
18 1 tesco(3,12), ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
19 2 iownd(14), ialth, ipup, lmax, meo, nqnyh, nslp, icf, ierpj,
20 3 iersl, jcur, jstart, kflag, l, meth, miter, maxord, maxcor,
21 4 msbp, mxncf, n, nq, nst, nfe, nje, nqu
22 INTEGER iownd, ialth, ipup, lmax, meo, nqnyh, nslp, icf, ierpj,
23 1 iersl, jcur, jstart, kflag, l, meth, miter, maxord, maxcor,
24 2 msbp, mxncf, n, nq, nst, nfe, nje, nqu
25 REAL(rprec) ::conit, crate, el, elco, hold, rmax, tesco,
26 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
30 REAL(rprec),
PARAMETER :: one = 1
31 INTEGER :: i1, iredo, iret, j, jb, m, ncf, newq
32 REAL(rprec) :: dcon, ddn, del, delp, dsm, dup, exdn, exsm,
33 1 exup, r, rh, rhdn, rhsm, rhup, told
37 REAL(rprec) ,
EXTERNAL :: vnorm
38 EXTERNAL f, jac, pjac, slvs
120 IF (jstart > 0)
GOTO 200
121 IF (jstart /= (-1))
THEN
122 IF (jstart == (-2))
GOTO 160
160 IF (ialth == 1) ialth = 2
161 IF (meth /= meo)
THEN
162 CALL cfode (meth, elco, tesco)
164 IF (nq > maxord)
GOTO 120
169 IF (nq <= maxord)
GOTO 160
177 conit = 0.5_dp/(nq + 2)
178 ddn = vnorm(n,savf,ewt)/tesco(1,l)
180 rhdn = one/(1.3_dp*ddn**exdn + 0.0000013_dp)
183 IF (h == hold)
GOTO 170
184 rh = min(rh,abs(h/hold))
193 CALL cfode (meth, elco, tesco)
199 conit = 0.5_dp/(nq + 2)
200 GOTO (160,170,200) iret
208 IF (h == hold)
GOTO 200
214 rh = max(rh,hmin/abs(h))
217 rh = rh/max(one,abs(h)*hmxi*rh)
221 yh(:n,j) = yh(:n,j)*r
226 IF (iredo == 0)
GOTO 690
236 IF (abs(rc - 1) > ccmax) ipup = miter
237 IF (nst >= nslp + msbp) ipup = miter
243 yh1(i1:nqnyh) = yh1(i1:nqnyh) + yh1(i1+nyh:nqnyh+nyh)
254 CALL f (neq, tn, y, savf)
256 IF (ipup <= 0)
GOTO 250
262 CALL pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac)
267 IF (ierpj /= 0)
GOTO 430
276 savf(:n) = h*savf(:n) - yh(:n,2)
277 y(:n) = savf(:n) - acor(:n)
279 y(:n) = yh(:n,1) + el(1)*savf(:n)
287 y(:n) = h*savf(:n) - (yh(:n,2)+acor(:n))
288 CALL slvs (wm, iwm, y, savf)
289 IF (iersl < 0)
GOTO 430
290 IF (iersl > 0)
GOTO 410
292 acor(:n) = acor(:n) + y(:n)
293 y(:n) = yh(:n,1) + el(1)*acor(:n)
299 IF (m /= 0) crate = max(0.2_dp*crate,del/delp)
300 dcon = del*min(one,1.5_dp*crate)/(tesco(2,nq)*conit)
301 IF (dcon <= one)
GOTO 450
303 IF (m == maxcor)
GOTO 410
304 IF (m>=2 .and. del>2.0_dp*delp)
GOTO 410
306 CALL f (neq, tn, y, savf)
317 IF (miter==0 .or. jcur==1)
GOTO 430
330 yh1(i1:nqnyh) = yh1(i1:nqnyh) - yh1(i1+nyh:nqnyh+nyh)
332 IF (ierpj<0 .or. iersl<0)
GOTO 680
333 IF (abs(h) <= hmin*1.00001_dp)
GOTO 670
334 IF (ncf == mxncf)
GOTO 670
347 IF (m == 0) dsm = del/tesco(2,nq)
348 IF (m > 0) dsm = vnorm(n,acor,ewt)/tesco(2,nq)
366 yh(:n,j) = yh(:n,j) + el(j)*acor(:n)
369 IF (ialth == 0)
GOTO 520
370 IF (ialth > 1)
GOTO 700
371 IF (l == lmax)
GOTO 700
372 yh(:n,lmax) = acor(:n)
388 yh1(i1:nqnyh) = yh1(i1:nqnyh) - yh1(i1+nyh:nqnyh+nyh)
391 IF (abs(h) <= hmin*1.00001_dp)
GOTO 660
392 IF (kflag <= (-3))
GOTO 640
408 savf(:n) = acor(:n) - yh(:n,lmax)
409 dup = vnorm(n,savf,ewt)/tesco(3,nq)
411 rhup = one/(1.4_dp*dup**exup + 0.0000014_dp)
415 rhsm = one/(1.2_dp*dsm**exsm + 0.0000012_dp)
418 ddn = vnorm(n,yh(1,l),ewt)/tesco(1,nq)
420 rhdn = one/(1.3_dp*ddn**exdn + 0.0000013_dp)
422 IF (rhsm < rhup)
THEN
423 IF (rhup > rhdn)
GOTO 590
425 IF (rhsm < rhdn)
GOTO 580
433 IF (kflag<0 .and. rh>one) rh = 1
438 IF (rh < 1.1_dp)
GOTO 610
440 yh(:n,newq+1) = acor(:n)*r
446 IF (kflag==0 .and. rh<1.1_dp)
GOTO 610
447 IF (kflag <= (-2)) rh = min(rh,0.2_dp)
453 IF (newq == nq)
GOTO 170
469 IF (kflag == (-10))
GOTO 660
471 rh = max(hmin/abs(h),rh)
474 CALL f (neq, tn, y, savf)
476 yh(:n,2) = h*savf(:n)
479 IF (nq == 1)
GOTO 200
501 acor(:n) = acor(:n)*r