1 SUBROUTINE lsode(f, neq, y, t, tout, itol, rtol, atol, itask,
2 1 istate, iopt, rwork, lrw, iwork, liw, jac, mf)
11 INTEGER :: itol, itask, istate, iopt, lrw, liw, jac, mf
12 REAL(rprec) :: t, tout
13 INTEGER,
DIMENSION(*) :: neq
14 INTEGER,
DIMENSION(liw) :: iwork
15 REAL(rprec),
DIMENSION(*) :: y, rtol, atol
16 REAL(rprec),
DIMENSION(lrw) :: rwork
21 COMMON /ls0001/ rowns(209),
22 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
23 2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
24 3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6),
25 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
26 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
27 INTEGER :: illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
28 1 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns,
29 2 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
30 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
31 REAL(rprec) ::rowns, ccmax, el0, h, hmin, hmxi, hu, rc, tn,
36 INTEGER,
PARAMETER :: mord(2) = (/12, 5 /)
37 INTEGER,
PARAMETER :: mxhnl0 = 10, mxstp0 = 500
38 REAL(rprec),
PARAMETER :: zero = 0, one = 1
39 INTEGER :: i,i1,i2,iflag,imxer,kgo,lf0,leniw,lenrw,lenwm,ml
41 REAL(rprec) :: atoli, ayi, big, ewti, h0, hmax, hmx, rh,
42 1 rtoli, tcrit, tdist, tnext, tol, tolsf, tp, size0, sum0, w0
47 REAL(rprec),
EXTERNAL :: prepj, solsy
48 REAL(rprec) ,
EXTERNAL :: vnorm
1008 IF (istate>=1 .and. istate<=3)
THEN
1009 IF (itask<1 .or. itask>5)
GOTO 602
1010 IF (istate /= 1)
THEN
1011 IF (init == 0)
GOTO 603
1012 IF (istate == 2)
GOTO 200
1015 IF (tout == t)
GOTO 430
1027 IF (neq(1) <= 0)
GOTO 604
1028 IF (istate == 1)
GOTO 25
1029 IF (neq(1) > n)
GOTO 605
1032 IF (itol<1 .or. itol>4)
GOTO 606
1033 IF (iopt<0 .or. iopt>1)
GOTO 607
1035 miter = mf - 10*meth
1036 IF (meth<1 .or. meth>2)
GOTO 608
1037 IF (miter<0 .or. miter>5)
GOTO 608
1038 IF (miter <= 3)
GOTO 30
1041 IF (ml<0 .or. ml>=n)
GOTO 609
1042 IF (mu<0 .or. mu>=n)
GOTO 610
1049 IF (istate == 1) h0 = 0
1054 IF (maxord < 0)
GOTO 611
1055 IF (maxord == 0) maxord = 100
1056 maxord = min(maxord,mord(meth))
1058 IF (mxstep < 0)
GOTO 612
1059 IF (mxstep == 0) mxstep = mxstp0
1061 IF (mxhnil < 0)
GOTO 613
1062 IF (mxhnil == 0) mxhnil = mxhnl0
1063 IF (istate /= 1)
GOTO 50
1065 IF ((tout - t)*h0 < zero)
GOTO 614
1068 IF (hmax < zero)
GOTO 615
1070 IF (hmax > zero) hmxi = 1/hmax
1072 IF (hmin < zero)
GOTO 616
1081 IF (istate == 1) nyh = n
1082 lwm = lyh + (maxord + 1)*nyh
1083 IF (miter == 0) lenwm = 0
1084 IF (miter==1 .or. miter==2) lenwm = n*n + 2
1085 IF (miter == 3) lenwm = n + 2
1086 IF (miter >= 4) lenwm = (2*ml + mu + 1)*n + 2
1090 lenrw = lacor + n - 1
1094 IF (miter==0 .or. miter==3) leniw = 20
1096 IF (lenrw > lrw)
GOTO 617
1097 IF (leniw > liw)
GOTO 618
1102 IF (itol >= 3) rtoli = rtol(i)
1103 IF (itol==2 .or. itol==4) atoli = atol(i)
1104 IF (rtoli < zero)
GOTO 619
1105 IF (atoli < zero)
GOTO 620
1107 IF (istate /= 1)
THEN
1110 IF (nq > maxord)
THEN
1113 rwork(i+lsavf-1) = rwork(i+lwm-1)
1117 IF (miter > 0) rwork(lwm) = sqrt(uround)
1118 IF (n == nyh)
GOTO 200
1121 i2 = lyh + (maxord + 1)*nyh - 1
1122 IF (i1 > i2)
GOTO 200
1134 uround = epsilon(uround)
1136 IF (itask==4 .or. itask==5)
THEN
1138 IF ((tcrit - tout)*(tout - t) < zero)
GOTO 625
1139 IF (h0/=zero .and. (t+h0-tcrit)*h0>zero) h0 = tcrit - t
1142 IF (miter > 0) rwork(lwm) = sqrt(uround)
1155 CALL f (neq, t, y, rwork(lf0))
1158 rwork(lyh:n-1+lyh) = y(:n)
1162 CALL ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1164 IF (rwork(i+lewt-1) <= zero)
GOTO 621
1165 rwork(i+lewt-1) = 1/rwork(i+lewt-1)
1183 IF (h0 == zero)
THEN
1184 tdist = abs(tout - t)
1185 w0 = max(abs(t),abs(tout))
1186 IF (tdist < 2*uround*w0)
GOTO 622
1189 tol = max(tol,maxval(rtol(:n)))
1191 IF (tol <= zero)
THEN
1194 IF (itol==2 .or. itol==4) atoli = atol(i)
1196 IF (ayi /= zero) tol = max(tol,atoli/ayi)
1199 tol = max(tol,100*uround)
1200 tol = min(tol,0.001_dp)
1201 sum0 = vnorm(n,rwork(lf0),rwork(lewt))
1202 sum0 = 1/(tol*w0*w0) + tol*sum0**2
1205 h0 = sign(h0,tout - t)
1209 IF (rh > one) h0 = h0/rh
1212 rwork(lf0:n-1+lf0) = h0*rwork(lf0:n-1+lf0)
1221 GOTO (210,250,220,230,240) itask
1223 IF ((tn - tout)*h < zero)
GOTO 250
1224 CALL intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1225 IF (iflag /= 0)
GOTO 627
1229 tp = tn - hu*(1 + 100*uround)
1230 IF ((tp - tout)*h > zero)
GOTO 623
1231 IF ((tn - tout)*h < zero)
GOTO 250
1235 IF ((tn - tcrit)*h > zero)
GOTO 624
1236 IF ((tcrit - tout)*h < zero)
GOTO 625
1237 IF ((tn - tout)*h < zero)
GOTO 245
1238 CALL intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1239 IF (iflag /= 0)
GOTO 627
1244 IF ((tn - tcrit)*h > zero)
GOTO 624
1246 hmx = abs(tn) + abs(h)
1247 ihit = abs(tn - tcrit) <= 100*uround*hmx
1249 tnext = tn + h*(1 + 4*uround)
1250 IF ((tnext - tcrit)*h > zero)
THEN
1251 h = (tcrit - tn)*(1 - 4*uround)
1252 IF (istate == 2) jstart = -2
1266 IF (nst - nslast >= mxstep)
GOTO 500
1267 CALL ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1269 IF (rwork(i+lewt-1) <= zero)
GOTO 510
1270 rwork(i+lewt-1) = 1/rwork(i+lewt-1)
1273 tolsf = uround*vnorm(n,rwork(lyh),rwork(lewt))
1274 IF (tolsf > one)
THEN
1276 IF (nst == 0)
GOTO 626
1279 IF (tn + h == tn)
THEN
1304 CALL stode (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1305 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm), f, jac
1308 GOTO (300,530,540) kgo
1316 GOTO (310,400,330,340,350) itask
1319 IF ((tn - tout)*h < zero)
GOTO 250
1320 CALL intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1325 IF ((tn - tout)*h >= zero)
GOTO 400
1329 IF ((tn - tout)*h >= zero)
THEN
1330 CALL intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1334 hmx = abs(tn) + abs(h)
1335 ihit = abs(tn - tcrit) <= 100*uround*hmx
1337 tnext = tn + h*(1 + 4*uround)
1338 IF ((tnext - tcrit)*h <= zero)
GOTO 250
1339 h = (tcrit - tn)*(1 - 4*uround)
1344 hmx = abs(tn) + abs(h)
1345 ihit = abs(tn - tcrit) <= 100*uround*hmx
1356 y(:n) = rwork(lyh:n-1+lyh)
1358 IF (itask==4 .or. itask==5)
THEN
1376 IF (ntrep < 5)
RETURN
1378 1
'lsode-- repeated calls with istate = 1 and tout = t (=r1) ',
1379 2 60, 301, 0, 0, 0, 0, 1, t, zero)
1393 1
'lsode-- at current t (=r1), mxstep (=i1) steps ', 50,
1394 2 201, 0, 0, 0, 0, 0, zero, zero)
1396 1
' taken on this call before reaching tout ', 50,
1397 2 201, 0, 1, mxstep, 0, 1, tn, zero)
1402 ewti = rwork(lewt+i-1)
1404 1 .le.
'lsode-- at t (=r1), ewt(i1) has become r2 0.', 50,
1405 2 202, 0, 1, i, 0, 2, tn, ewti)
1411 1
'lsode-- at t (=r1), too much accuracy requested ', 50,
1412 2 203, 0, 0, 0, 0, 0, zero, zero)
1414 1
' for precision of machine.. see tolsf (=r2) ', 50,
1415 2 203, 0, 0, 0, 0, 2, tn, tolsf)
1422 1
'lsode-- at t(=r1) and step size h(=r2), the error', 50,
1423 2 204, 0, 0, 0, 0, 0, zero, zero)
1425 1
' test failed repeatedly or with abs(h) = hmin', 50,
1426 2 204, 0, 0, 0, 0, 2, tn, h)
1432 1
'lsode-- at t (=r1) and step size h (=r2), the ', 50,
1433 2 205, 0, 0, 0, 0, 0, zero, zero)
1435 1
' corrector convergence failed repeatedly ', 50,
1436 2 205, 0, 0, 0, 0, 0, zero, zero)
1437 CALL xerrwv (
' or with abs(h) = hmin ', 30, 205, 0, 0,
1445 size0 = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
1446 IF (big < size0)
THEN
1454 y(:n) = rwork(lyh:n-1+lyh)
1475 CALL xerrwv (
'lsode-- istate (=i1) illegal ', 30, 1, 0, 1,
1476 1 istate, 0, 0, zero, zero)
1479 CALL xerrwv (
'lsode-- itask (=i1) illegal ', 30, 2, 0, 1,
1480 1 itask, 0, 0, zero, zero)
1483 CALL xerrwv (.gt.
'lsode-- istate 1 but lsode not initialized '
1484 1 , 50, 3, 0, 0, 0, 0, 0, zero, zero)
1487 CALL xerrwv (.lt.
'lsode-- neq (=i1) 1 ', 30, 4, 0, 1,
1488 1 neq(1), 0, 0, zero, zero)
1491 CALL xerrwv (
'lsode-- istate = 3 and neq increased (i1 to i2) '
1492 1 , 50, 5, 0, 2, n, neq(1), 0, zero, zero)
1495 CALL xerrwv (
'lsode-- itol (=i1) illegal ', 30, 6, 0, 1, itol
1496 1 , 0, 0, zero, zero)
1499 CALL xerrwv (
'lsode-- iopt (=i1) illegal ', 30, 7, 0, 1, iopt
1500 1 , 0, 0, zero, zero)
1503 CALL xerrwv (
'lsode-- mf (=i1) illegal ', 30, 8, 0, 1, mf, 0
1507 CALL xerrwv (.lt..ge.
'lsode-- ml (=i1) illegal.. 0 or neq (=i2)'
1508 1 , 50, 9, 0, 2, ml, neq(1), 0, zero, zero)
1511 CALL xerrwv (.lt..ge.
'lsode-- mu (=i1) illegal.. 0 or neq (=i2)'
1512 1 , 50, 10, 0, 2, mu, neq(1), 0, zero, zero)
1515 CALL xerrwv (.lt.
'lsode-- maxord (=i1) 0 ', 30, 11, 0, 1,
1516 1 maxord, 0, 0, zero, zero)
1519 CALL xerrwv (.lt.
'lsode-- mxstep (=i1) 0 ', 30, 12, 0, 1,
1520 1 mxstep, 0, 0, zero, zero)
1523 CALL xerrwv (.lt.
'lsode-- mxhnil (=i1) 0 ', 30, 13, 0, 1,
1524 1 mxhnil, 0, 0, zero, zero)
1527 CALL xerrwv (
'lsode-- tout (=r1) behind t (=r2) ', 40, 14,
1528 1 0, 0, 0, 0, 2, tout, t)
1529 CALL xerrwv (
' integration direction is given by h0 (=r1) '
1530 1 , 50, 14, 0, 0, 0, 0, 1, h0, zero)
1533 CALL xerrwv (.lt.
'lsode-- hmax (=r1) 0.0 ', 30, 15, 0, 0, 0, 0
1537 CALL xerrwv (.lt.
'lsode-- hmin (=r1) 0.0 ', 30, 16, 0, 0, 0, 0
1542 1
'lsode-- rwork length needed, lenrw (=i1), exceeds lrw (=i2)',
1543 2 60, 17, 0, 2, lenrw, lrw, 0, zero, zero)
1547 1
'lsode-- iwork length needed, leniw (=i1), exceeds liw (=i2)',
1548 2 60, 18, 0, 2, leniw, liw, 0, zero, zero)
1551 CALL xerrwv (.lt.
'lsode-- rtol(i1) is r1 0.0 ', 40, 19,
1552 1 0, 1, i, 0, 1, rtoli, zero)
1555 CALL xerrwv (.lt.
'lsode-- atol(i1) is r1 0.0 ', 40, 20,
1556 1 0, 1, i, 0, 1, atoli, zero)
1559 ewti = rwork(lewt+i-1)
1560 CALL xerrwv (.le.
'lsode-- ewt(i1) is r1 0.0 ', 40, 21,
1561 1 0, 1, i, 0, 1, ewti, zero)
1565 1
'lsode-- tout (=r1) too close to t(=r2) to start integration',
1566 2 60, 22, 0, 0, 0, 0, 2, tout, t)
1570 1
'lsode-- itask = i1 and tout (=r1) behind tcur - hu (= r2) ',
1571 2 60, 23, 0, 1, itask, 0, 2, tout, tp)
1575 1
'lsode-- itask = 4 or 5 and tcrit (=r1) behind tcur (=r2) ',
1576 2 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1580 1
'lsode-- itask = 4 or 5 and tcrit (=r1) behind tout (=r2) ',
1581 2 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1584 CALL xerrwv (
'lsode-- at start of problem, too much accuracy '
1585 1 , 50, 26, 0, 0, 0, 0, 0, zero, zero)
1587 1
' requested for precision of machine.. see tolsf (=r1) ',
1588 2 60, 26, 0, 0, 0, 0, 1, tolsf, zero)
1592 CALL xerrwv (
'lsode-- trouble from intdy. itask = i1, tout = r1'
1593 1 , 50, 27, 0, 1, itask, 0, 1, tout, zero)
1596 IF (illin /= 5)
THEN
1601 CALL xerrwv (
'lsode-- repeated occurrences of illegal input '
1602 1 , 50, 302, 0, 0, 0, 0, 0, zero, zero)
1605 CALL xerrwv (
'lsode-- run aborted.. apparent infinite loop '
1606 1 , 50, 303, 2, 0, 0, 0, 0, zero, zero)
1608 END SUBROUTINE lsode