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