1 SUBROUTINE STOD(NEQ,Y,YH,NYH,YH1,EWT,SAVF,ACOR,WM,IWM,F,JAC,RPAR, 2 1 IPAR) 3C***BEGIN PROLOGUE STOD 4C***REFER TO DEBDF 5C 6C STOD integrates a system of first order odes over one step in the 7C integrator package DEBDF. 8C ---------------------------------------------------------------------- 9C STOD performs one step of the integration of an initial value 10C problem for a system of ordinary differential equations. 11C Note.. STOD is independent of the value of the iteration method 12C indicator MITER, when this is .NE. 0, and hence is independent 13C of the type of chord method used, or the Jacobian structure. 14C Communication with STOD is done with the following variables.. 15C 16C Y = An array of length .GE. n used as the Y argument in 17C all calls to F and JAC. 18C NEQ = Integer array containing problem size in NEQ(1), and 19C passed as the NEQ argument in all calls to F and JAC. 20C YH = An NYH by LMAX array containing the dependent variables 21C and their approximate scaled derivatives, where 22C LMAX = MAXORD + 1. YH(I,J+1) contains the approximate 23C J-th derivative of Y(I), scaled by H**J/Factorial(j) 24C (J = 0,1,...,NQ). On entry for the first step, the first 25C two columns of YH must be set from the initial values. 26C NYH = A constant integer .GE. N, the first dimension of YH. 27C YH1 = A one-dimensional array occupying the same space as YH. 28C EWT = An array of N elements with which the estimated local 29C errors in YH are compared. 30C SAVF = An array of working storage, of length N. 31C ACOR = A work array of length N, used for the accumulated 32C corrections. On a successful return, ACOR(I) contains 33C the estimated one-step local error in Y(I). 34C WM,IWM = Real and integer work arrays associated with matrix 35C operations in chord iteration (MITER .NE. 0). 36C PJAC = Name of routine to evaluate and preprocess Jacobian matrix 37C if a chord method is being used. 38C SLVS = Name of routine to solve linear system in chord iteration. 39C H = The step size to be attempted on the next step. 40C H is altered by the error control algorithm during the 41C problem. H can be either positive or negative, but its 42C sign must remain constant throughout the problem. 43C HMIN = The minimum absolute value of the step size H to be used. 44C HMXI = Inverse of the maximum absolute value of H to be used. 45C HMXI = 0.0 is allowed and corresponds to an infinite HMAX. 46C HMIN and HMXI may be changed at any time, but will not 47C take effect until the next change of H is considered. 48C TN = The independent variable. TN is updated on each step taken. 49C JSTART = An integer used for input only, with the following 50C values and meanings.. 51C 0 Perform the first step. 52C .GT.0 Take a new step continuing from the last. 53C -1 Take the next step with a new value of H, MAXORD, 54C N, METH, MITER, and/or matrix parameters. 55C -2 Take the next step with a new value of H, 56C but with other inputs unchanged. 57C On return, JSTART is set to 1 to facilitate continuation. 58C KFLAG = a completion code with the following meanings.. 59C 0 The step was succesful. 60C -1 The requested error could not be achieved. 61C -2 Corrector convergence could not be achieved. 62C A return with KFLAG = -1 or -2 means either 63C ABS(H) = HMIN or 10 consecutive failures occurred. 64C On a return with KFLAG negative, the values of TN and 65C the YH array are as of the beginning of the last 66C step, and H is the last step size attempted. 67C MAXORD = The maximum order of integration method to be allowed. 68C METH/MITER = The method flags. See description in driver. 69C N = The number of first-order differential equations. 70C ---------------------------------------------------------------------- 71C***ROUTINES CALLED CFOD,PJAC,SLVS,VNWRMS 72C***COMMON BLOCKS DEBDF1 73C***REVISION HISTORY (YYMMDD) 74C 000330 Modified array declarations. (JEC) 75C 76C***END PROLOGUE STOD 77 EXTERNAL F, JAC 78C 79CLLL. OPTIMIZE 80 INTEGER NEQ, NYH, IWM, I, I1, IALTH, IER, IOWND, IREDO, IRET, 81 1 IPUP, J, JB, JSTART, KFLAG, L, LMAX, M, MAXORD, MEO, METH, 82 2 MITER, N, NCF, NEWQ, NFE, NJE, NQ, NQNYH, NQU, NST, NSTEPJ 83 REAL Y, YH, YH1, EWT, SAVF, ACOR, WM, 84 1 ROWND, CONIT, CRATE, EL, ELCO, HOLD, RC, RMAX, TESCO, 85 2 EL0, H, HMIN, HMXI, HU, TN, UROUND, 86 3 DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, 87 4 R, RH, RHDN, RHSM, RHUP, TOLD, VNWRMS 88 DIMENSION Y(*), YH(NYH,MAXORD+1), YH1(*), EWT(*), SAVF(*), 89 1 ACOR(*), WM(*), IWM(*), RPAR(*), IPAR(*) 90 COMMON /DEBDF1/ ROWND, CONIT, CRATE, EL(13), ELCO(13,12), 91 1 HOLD, RC, RMAX, TESCO(3,12), 92 2 EL0, H, HMIN, HMXI, HU, TN, UROUND, IOWND(7), KSTEPS, IOD(6), 93 3 IALTH, IPUP, LMAX, MEO, NQNYH, NSTEPJ, 94 4 IER, JSTART, KFLAG, L, METH, MITER, MAXORD, N, NQ, NST, NFE, 95 5 NJE, NQU 96C 97C 98C***FIRST EXECUTABLE STATEMENT STOD 99 KFLAG = 0 100 TOLD = TN 101 NCF = 0 102 IF (JSTART .GT. 0) GO TO 200 103 IF (JSTART .EQ. -1) GO TO 100 104 IF (JSTART .EQ. -2) GO TO 160 105C----------------------------------------------------------------------- 106C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE 107C INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED 108C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL 109C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE 110C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2 111C FOR THE NEXT INCREASE. 112C----------------------------------------------------------------------- 113 LMAX = MAXORD + 1 114 NQ = 1 115 L = 2 116 IALTH = 2 117 RMAX = 10000.0E0 118 RC = 0.0E0 119 EL0 = 1.0E0 120 CRATE = 0.7E0 121 DELP = 0.0E0 122 HOLD = H 123 MEO = METH 124 NSTEPJ = 0 125 IRET = 3 126 GO TO 140 127C----------------------------------------------------------------------- 128C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1. 129C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE. 130C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1), 131C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP. 132C IF THE CALLER HAS CHANGED METH, CFOD IS CALLED TO RESET 133C THE COEFFICIENTS OF THE METHOD. 134C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT 135C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY. 136C IF H IS TO BE CHANGED, YH MUST BE RESCALED. 137C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1 138C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS. 139C----------------------------------------------------------------------- 140 100 IPUP = MITER 141 LMAX = MAXORD + 1 142 IF (IALTH .EQ. 1) IALTH = 2 143 IF (METH .EQ. MEO) GO TO 110 144 CALL CFOD (METH, ELCO, TESCO) 145 MEO = METH 146 IF (NQ .GT. MAXORD) GO TO 120 147 IALTH = L 148 IRET = 1 149 GO TO 150 150 110 IF (NQ .LE. MAXORD) GO TO 160 151 120 NQ = MAXORD 152 L = LMAX 153 DO 125 I = 1,L 154 125 EL(I) = ELCO(I,NQ) 155 NQNYH = NQ*NYH 156 RC = RC*EL(1)/EL0 157 EL0 = EL(1) 158 CONIT = 0.5E0/FLOAT(NQ+2) 159 DDN = VNWRMS (N, SAVF, EWT)/TESCO(1,L) 160 EXDN = 1.0E0/FLOAT(L) 161 RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) 162 RH = AMIN1(RHDN,1.0E0) 163 IREDO = 3 164 IF (H .EQ. HOLD) GO TO 170 165 RH = AMIN1(RH,ABS(H/HOLD)) 166 H = HOLD 167 GO TO 175 168C----------------------------------------------------------------------- 169C CFOD IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE 170C CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET 171C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM. 172C----------------------------------------------------------------------- 173 140 CALL CFOD (METH, ELCO, TESCO) 174 150 DO 155 I = 1,L 175 155 EL(I) = ELCO(I,NQ) 176 NQNYH = NQ*NYH 177 RC = RC*EL(1)/EL0 178 EL0 = EL(1) 179 CONIT = 0.5E0/FLOAT(NQ+2) 180 GO TO (160, 170, 200), IRET 181C----------------------------------------------------------------------- 182C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST 183C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO 184C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS 185C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE. 186C----------------------------------------------------------------------- 187 160 IF (H .EQ. HOLD) GO TO 200 188 RH = H/HOLD 189 H = HOLD 190 IREDO = 3 191 GO TO 175 192 170 RH = AMAX1(RH,HMIN/ABS(H)) 193 175 RH = AMIN1(RH,RMAX) 194 RH = RH/AMAX1(1.0E0,ABS(H)*HMXI*RH) 195 R = 1.0E0 196 DO 180 J = 2,L 197 R = R*RH 198 DO 180 I = 1,N 199 180 YH(I,J) = YH(I,J)*R 200 H = H*RH 201 RC = RC*RH 202 IALTH = L 203 IF (IREDO .EQ. 0) GO TO 680 204C----------------------------------------------------------------------- 205C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY 206C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX. 207C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1). 208C WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER 209C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED. 210C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY 20-TH STEP. 211C----------------------------------------------------------------------- 212 200 IF (ABS(RC-1.0E0) .GT. 0.3E0) IPUP = MITER 213 IF (NST .GE. NSTEPJ+20) IPUP = MITER 214 TN = TN + H 215 I1 = NQNYH + 1 216 DO 215 JB = 1,NQ 217 I1 = I1 - NYH 218 DO 210 I = I1,NQNYH 219 210 YH1(I) = YH1(I) + YH1(I+NYH) 220 215 CONTINUE 221 KSTEPS = KSTEPS + 1 222C----------------------------------------------------------------------- 223C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS 224C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR 225C WEIGHT VECTOR EWT. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE 226C VECTOR ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP. 227C----------------------------------------------------------------------- 228 220 M = 0 229 DO 230 I = 1,N 230 230 Y(I) = YH(I,1) 231 CALL F (TN, Y, SAVF, RPAR, IPAR) 232 NFE = NFE + 1 233 IF (IPUP .LE. 0) GO TO 250 234C----------------------------------------------------------------------- 235C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND 236C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET 237C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE. 238C----------------------------------------------------------------------- 239 IPUP = 0 240 RC = 1.0E0 241 NSTEPJ = NST 242 CRATE = 0.7E0 243 CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC, 244 1 RPAR, IPAR) 245 IF (IER .NE. 0) GO TO 430 246 250 DO 260 I = 1,N 247 260 ACOR(I) = 0.0E0 248 270 IF (MITER .NE. 0) GO TO 350 249C----------------------------------------------------------------------- 250C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM 251C THE RESULT OF THE LAST FUNCTION EVALUATION. 252C----------------------------------------------------------------------- 253 DO 290 I = 1,N 254 SAVF(I) = H*SAVF(I) - YH(I,2) 255 290 Y(I) = SAVF(I) - ACOR(I) 256 DEL = VNWRMS (N, Y, EWT) 257 DO 300 I = 1,N 258 Y(I) = YH(I,1) + EL(1)*SAVF(I) 259 300 ACOR(I) = SAVF(I) 260 GO TO 400 261C----------------------------------------------------------------------- 262C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR, 263C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND 264C P AS COEFFICIENT MATRIX. 265C----------------------------------------------------------------------- 266 350 DO 360 I = 1,N 267 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I)) 268 CALL SLVS (WM, IWM, Y, SAVF) 269 IF (IER .NE. 0) GO TO 410 270 DEL = VNWRMS (N, Y, EWT) 271 DO 380 I = 1,N 272 ACOR(I) = ACOR(I) + Y(I) 273 380 Y(I) = YH(I,1) + EL(1)*ACOR(I) 274C----------------------------------------------------------------------- 275C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE 276C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. 277C----------------------------------------------------------------------- 278 400 IF (M .NE. 0) CRATE = AMAX1(0.2E0*CRATE,DEL/DELP) 279 DCON = DEL*AMIN1(1.0E0,1.5E0*CRATE)/(TESCO(2,NQ)*CONIT) 280 IF (DCON .LE. 1.0E0) GO TO 450 281 M = M + 1 282 IF (M .EQ. 3) GO TO 410 283 IF (M .GE. 2 .AND. DEL .GT. 2.0E0*DELP) GO TO 410 284 DELP = DEL 285 CALL F (TN, Y, SAVF, RPAR, IPAR) 286 NFE = NFE + 1 287 GO TO 270 288C----------------------------------------------------------------------- 289C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. 290C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR 291C THE NEXT TRY. OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES 292C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF H CANNOT BE 293C REDUCED OR 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2. 294C----------------------------------------------------------------------- 295 410 IF (IPUP .EQ. 0) GO TO 430 296 IPUP = MITER 297 GO TO 220 298 430 TN = TOLD 299 NCF = NCF + 1 300 RMAX = 2.0E0 301 I1 = NQNYH + 1 302 DO 445 JB = 1,NQ 303 I1 = I1 - NYH 304 DO 440 I = I1,NQNYH 305 440 YH1(I) = YH1(I) - YH1(I+NYH) 306 445 CONTINUE 307 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 670 308 IF (NCF .EQ. 10) GO TO 670 309 RH = 0.25E0 310 IPUP = MITER 311 IREDO = 1 312 GO TO 170 313C----------------------------------------------------------------------- 314C THE CORRECTOR HAS CONVERGED. IPUP IS SET TO -1 IF MITER .NE. 0, 315C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER. 316C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500 317C IF IT FAILS. 318C----------------------------------------------------------------------- 319 450 IF (MITER .NE. 0) IPUP = -1 320 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) 321 IF (M .GT. 0) DSM = VNWRMS (N, ACOR, EWT)/TESCO(2,NQ) 322 IF (DSM .GT. 1.0E0) GO TO 500 323C----------------------------------------------------------------------- 324C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY. 325C CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1. 326C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR 327C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP. 328C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER 329C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A 330C FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT 331C TESTING FOR THAT MANY STEPS. 332C----------------------------------------------------------------------- 333 KFLAG = 0 334 IREDO = 0 335 NST = NST + 1 336 HU = H 337 NQU = NQ 338 DO 470 J = 1,L 339 DO 470 I = 1,N 340 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I) 341 IALTH = IALTH - 1 342 IF (IALTH .EQ. 0) GO TO 520 343 IF (IALTH .GT. 1) GO TO 690 344 IF (L .EQ. LMAX) GO TO 690 345 DO 490 I = 1,N 346 490 YH(I,LMAX) = ACOR(I) 347 GO TO 690 348C----------------------------------------------------------------------- 349C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES. 350C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE 351C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR 352C ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE 353C BY A FACTOR OF 0.2 OR LESS. 354C----------------------------------------------------------------------- 355 500 KFLAG = KFLAG - 1 356 TN = TOLD 357 I1 = NQNYH + 1 358 DO 515 JB = 1,NQ 359 I1 = I1 - NYH 360 DO 510 I = I1,NQNYH 361 510 YH1(I) = YH1(I) - YH1(I+NYH) 362 515 CONTINUE 363 RMAX = 2.0E0 364 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 660 365 IF (KFLAG .LE. -3) GO TO 640 366 IREDO = 2 367 RHUP = 0.0E0 368 GO TO 540 369C----------------------------------------------------------------------- 370C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS 371C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED 372C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY. 373C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE. 374C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN 375C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE 376C ADDITIONAL SCALED DERIVATIVE. 377C----------------------------------------------------------------------- 378 520 RHUP = 0.0E0 379 IF (L .EQ. LMAX) GO TO 540 380 DO 530 I = 1,N 381 530 SAVF(I) = ACOR(I) - YH(I,LMAX) 382 DUP = VNWRMS (N, SAVF, EWT)/TESCO(3,NQ) 383 EXUP = 1.0E0/FLOAT(L+1) 384 RHUP = 1.0E0/(1.4E0*DUP**EXUP + 0.0000014E0) 385 540 EXSM = 1.0E0/FLOAT(L) 386 RHSM = 1.0E0/(1.2E0*DSM**EXSM + 0.0000012E0) 387 RHDN = 0.0E0 388 IF (NQ .EQ. 1) GO TO 560 389 DDN = VNWRMS (N, YH(1,L), EWT)/TESCO(1,NQ) 390 EXDN = 1.0E0/FLOAT(NQ) 391 RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) 392 560 IF (RHSM .GE. RHUP) GO TO 570 393 IF (RHUP .GT. RHDN) GO TO 590 394 GO TO 580 395 570 IF (RHSM .LT. RHDN) GO TO 580 396 NEWQ = NQ 397 RH = RHSM 398 GO TO 620 399 580 NEWQ = NQ - 1 400 RH = RHDN 401 IF (KFLAG .LT. 0 .AND. RH .GT. 1.0E0) RH = 1.0E0 402 GO TO 620 403 590 NEWQ = L 404 RH = RHUP 405 IF (RH .LT. 1.1E0) GO TO 610 406 R = EL(L)/FLOAT(L) 407 DO 600 I = 1,N 408 600 YH(I,NEWQ+1) = ACOR(I)*R 409 GO TO 630 410 610 IALTH = 3 411 GO TO 690 412 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1E0)) GO TO 610 413 IF (KFLAG .LE. -2) RH = AMIN1(RH,0.2E0) 414C----------------------------------------------------------------------- 415C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. 416C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED. 417C THEN EXIT FROM 680 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE. 418C----------------------------------------------------------------------- 419 IF (NEWQ .EQ. NQ) GO TO 170 420 630 NQ = NEWQ 421 L = NQ + 1 422 IRET = 2 423 GO TO 150 424C----------------------------------------------------------------------- 425C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED. 426C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1. 427C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE 428C YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST 429C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN 430C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED, 431C UNTIL IT SUCCEEDS OR H REACHES HMIN. 432C----------------------------------------------------------------------- 433 640 IF (KFLAG .EQ. -10) GO TO 660 434 RH = 0.1E0 435 RH = AMAX1(HMIN/ABS(H),RH) 436 H = H*RH 437 DO 645 I = 1,N 438 645 Y(I) = YH(I,1) 439 CALL F (TN, Y, SAVF, RPAR, IPAR) 440 NFE = NFE + 1 441 DO 650 I = 1,N 442 650 YH(I,2) = H*SAVF(I) 443 IPUP = MITER 444 IALTH = 5 445 IF (NQ .EQ. 1) GO TO 200 446 NQ = 1 447 L = 2 448 IRET = 3 449 GO TO 150 450C----------------------------------------------------------------------- 451C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD 452C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. 453C----------------------------------------------------------------------- 454 660 KFLAG = -1 455 GO TO 700 456 670 KFLAG = -2 457 GO TO 700 458 680 RMAX = 10.0E0 459 690 R = 1.0E0/TESCO(2,NQU) 460 DO 695 I = 1,N 461 695 ACOR(I) = ACOR(I)*R 462 700 HOLD = H 463 JSTART = 1 464 RETURN 465C----------------------- END OF SUBROUTINE STOD ----------------------- 466 END 467