1      SUBROUTINE STODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
2     1   WM, IWM, F, JAC, PJAC, SLVS, IERR)
3CLLL. OPTIMIZE
4      EXTERNAL F, JAC, PJAC, SLVS
5      INTEGER NEQ, NYH, IWM
6      INTEGER ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
7     1   MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, CNYH,
8     2   IALTH, IPUP, LMAX, MEO, NQNYH, NSLP
9      INTEGER ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
10     1   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
11      INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
12      DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM
13      DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
14     2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
15      DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
16     1   R, RH, RHDN, RHSM, RHUP, TOLD, VNORM
17      DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),
18     1   ACOR(*), WM(*), IWM(*)
19      COMMON /LS0001/ CONIT, CRATE, EL(13), ELCO(13,12),
20     1   HOLD, RMAX, TESCO(3,12),
21     2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
22     2   ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
23     3   MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, CNYH,
24     3   IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
25     4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
26     5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
27C-----------------------------------------------------------------------
28C STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
29C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
30C NOTE.. STODE IS INDEPENDENT OF THE VALUE OF THE ITERATION METHOD
31C INDICATOR MITER, WHEN THIS IS .NE. 0, AND HENCE IS INDEPENDENT
32C OF THE TYPE OF CHORD METHOD USED, OR THE JACOBIAN STRUCTURE.
33C COMMUNICATION WITH STODE IS DONE WITH THE FOLLOWING VARIABLES..
34C
35C NEQ    = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND
36C          PASSED AS THE NEQ ARGUMENT IN ALL CALLS TO F AND JAC.
37C Y      = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN
38C          ALL CALLS TO F AND JAC.
39C YH     = AN NYH BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
40C          AND THEIR APPROXIMATE SCALED DERIVATIVES, WHERE
41C          LMAX = MAXORD + 1.  YH(I,J+1) CONTAINS THE APPROXIMATE
42C          J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J)
43C          (J = 0,1,...,NQ).  ON ENTRY FOR THE FIRST STEP, THE FIRST
44C          TWO COLUMNS OF YH MUST BE SET FROM THE INITIAL VALUES.
45C NYH    = A CONSTANT INTEGER .GE. N, THE FIRST DIMENSION OF YH.
46C YH1    = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH.
47C EWT    = AN ARRAY OF LENGTH N CONTAINING MULTIPLICATIVE WEIGHTS
48C          FOR LOCAL ERROR MEASUREMENTS.  LOCAL ERRORS IN Y(I) ARE
49C          COMPARED TO 1.0/EWT(I) IN VARIOUS ERROR TESTS.
50C SAVF   = AN ARRAY OF WORKING STORAGE, OF LENGTH N.
51C          ALSO USED FOR INPUT OF YH(*,MAXORD+2) WHEN JSTART = -1
52C          AND MAXORD .LT. THE CURRENT ORDER NQ.
53C ACOR   = A WORK ARRAY OF LENGTH N, USED FOR THE ACCUMULATED
54C          CORRECTIONS.  ON A SUCCESSFUL RETURN, ACOR(I) CONTAINS
55C          THE ESTIMATED ONE-STEP LOCAL ERROR IN Y(I).
56C WM,IWM = REAL AND INTEGER WORK ARRAYS ASSOCIATED WITH MATRIX
57C          OPERATIONS IN CHORD ITERATION (MITER .NE. 0).
58C PJAC   = NAME OF ROUTINE TO EVALUATE AND PREPROCESS JACOBIAN MATRIX
59C          AND P = I - H*EL0*JAC, IF A CHORD METHOD IS BEING USED.
60C SLVS   = NAME OF ROUTINE TO SOLVE LINEAR SYSTEM IN CHORD ITERATION.
61C CCMAX  = MAXIMUM RELATIVE CHANGE IN H*EL0 BEFORE PJAC IS CALLED.
62C H      = THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
63C          H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
64C          PROBLEM.  H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
65C          SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
66C HMIN   = THE MINIMUM ABSOLUTE VALUE OF THE STEP SIZE H TO BE USED.
67C HMXI   = INVERSE OF THE MAXIMUM ABSOLUTE VALUE OF H TO BE USED.
68C          HMXI = 0.0 IS ALLOWED AND CORRESPONDS TO AN INFINITE HMAX.
69C          HMIN AND HMXI MAY BE CHANGED AT ANY TIME, BUT WILL NOT
70C          TAKE EFFECT UNTIL THE NEXT CHANGE OF H IS CONSIDERED.
71C TN     = THE INDEPENDENT VARIABLE. TN IS UPDATED ON EACH STEP TAKEN.
72C JSTART = AN INTEGER USED FOR INPUT ONLY, WITH THE FOLLOWING
73C          VALUES AND MEANINGS..
74C               0  PERFORM THE FIRST STEP.
75C           .GT.0  TAKE A NEW STEP CONTINUING FROM THE LAST.
76C              -1  TAKE THE NEXT STEP WITH A NEW VALUE OF H, MAXORD,
77C                    N, METH, MITER, AND/OR MATRIX PARAMETERS.
78C              -2  TAKE THE NEXT STEP WITH A NEW VALUE OF H,
79C                    BUT WITH OTHER INPUTS UNCHANGED.
80C          ON RETURN, JSTART IS SET TO 1 TO FACILITATE CONTINUATION.
81C KFLAG  = A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
82C               0  THE STEP WAS SUCCESFUL.
83C              -1  THE REQUESTED ERROR COULD NOT BE ACHIEVED.
84C              -2  CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED.
85C              -3  FATAL ERROR IN PJAC OR SLVS.
86C          A RETURN WITH KFLAG = -1 OR -2 MEANS EITHER
87C          ABS(H) = HMIN OR 10 CONSECUTIVE FAILURES OCCURRED.
88C          ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF TN AND
89C          THE YH ARRAY ARE AS OF THE BEGINNING OF THE LAST
90C          STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
91C MAXORD = THE MAXIMUM ORDER OF INTEGRATION METHOD TO BE ALLOWED.
92C MAXCOR = THE MAXIMUM NUMBER OF CORRECTOR ITERATIONS ALLOWED.
93C MSBP   = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0).
94C MXNCF  = MAXIMUM NUMBER OF CONVERGENCE FAILURES ALLOWED.
95C METH/MITER = THE METHOD FLAGS.  SEE DESCRIPTION IN DRIVER.
96C N      = THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
97C IERR   = ERROR FLAG FROM USER-SUPPLIED FUNCTION
98C-----------------------------------------------------------------------
99      KFLAG = 0
100      TOLD = TN
101      NCF = 0
102      IERPJ = 0
103      IERSL = 0
104      JCUR = 0
105      ICF = 0
106      DELP = 0.0D0
107      IF (JSTART .GT. 0) GO TO 200
108      IF (JSTART .EQ. -1) GO TO 100
109      IF (JSTART .EQ. -2) GO TO 160
110C-----------------------------------------------------------------------
111C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE
112C INITIALIZED.  RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
113C IN A SINGLE STEP.  IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
114C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10.  IF A FAILURE
115C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
116C FOR THE NEXT INCREASE.
117C-----------------------------------------------------------------------
118      LMAX = MAXORD + 1
119      NQ = 1
120      L = 2
121      IALTH = 2
122      RMAX = 10000.0D0
123      RC = 0.0D0
124      EL0 = 1.0D0
125      CRATE = 0.7D0
126      HOLD = H
127      MEO = METH
128      NSLP = 0
129      IPUP = MITER
130      IRET = 3
131      GO TO 140
132C-----------------------------------------------------------------------
133C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1.
134C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE.
135C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1),
136C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP.
137C IF THE CALLER HAS CHANGED METH, CFODE IS CALLED TO RESET
138C THE COEFFICIENTS OF THE METHOD.
139C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
140C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY.
141C IF H IS TO BE CHANGED, YH MUST BE RESCALED.
142C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1
143C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS.
144C-----------------------------------------------------------------------
145 100  IPUP = MITER
146      LMAX = MAXORD + 1
147      IF (IALTH .EQ. 1) IALTH = 2
148      IF (METH .EQ. MEO) GO TO 110
149      CALL CFODE (METH, ELCO, TESCO)
150      MEO = METH
151      IF (NQ .GT. MAXORD) GO TO 120
152      IALTH = L
153      IRET = 1
154      GO TO 150
155 110  IF (NQ .LE. MAXORD) GO TO 160
156 120  NQ = MAXORD
157      L = LMAX
158      DO 125 I = 1,L
159 125    EL(I) = ELCO(I,NQ)
160      NQNYH = NQ*NYH
161      RC = RC*EL(1)/EL0
162      EL0 = EL(1)
163      CONIT = 0.5D0/DBLE(NQ+2)
164      DDN = VNORM (N, SAVF, EWT)/TESCO(1,L)
165      EXDN = 1.0D0/DBLE(L)
166      RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
167      RH = DMIN1(RHDN,1.0D0)
168      IREDO = 3
169      IF (H .EQ. HOLD) GO TO 170
170      RH = DMIN1(RH,DABS(H/HOLD))
171      H = HOLD
172      GO TO 175
173C-----------------------------------------------------------------------
174C CFODE IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE
175C CURRENT METH.  THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET
176C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM.
177C-----------------------------------------------------------------------
178 140  CALL CFODE (METH, ELCO, TESCO)
179 150  DO 155 I = 1,L
180 155    EL(I) = ELCO(I,NQ)
181      NQNYH = NQ*NYH
182      RC = RC*EL(1)/EL0
183      EL0 = EL(1)
184      CONIT = 0.5D0/DBLE(NQ+2)
185      GO TO (160, 170, 200), IRET
186C-----------------------------------------------------------------------
187C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
188C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED.  IALTH IS SET TO
189C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS
190C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE.
191C-----------------------------------------------------------------------
192 160  IF (H .EQ. HOLD) GO TO 200
193      RH = H/HOLD
194      H = HOLD
195      IREDO = 3
196      GO TO 175
197 170  RH = DMAX1(RH,HMIN/DABS(H))
198 175  RH = DMIN1(RH,RMAX)
199      RH = RH/DMAX1(1.0D0,DABS(H)*HMXI*RH)
200      R = 1.0D0
201      DO 180 J = 2,L
202        R = R*RH
203        DO 180 I = 1,N
204 180      YH(I,J) = YH(I,J)*R
205      H = H*RH
206      RC = RC*RH
207      IALTH = L
208      IF (IREDO .EQ. 0) GO TO 690
209C-----------------------------------------------------------------------
210C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
211C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX.
212C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT  H*EL(1).
213C WHEN RC DIFFERS FROM 1 BY MORE THAN CCMAX, IPUP IS SET TO MITER
214C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED.
215C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS.
216C-----------------------------------------------------------------------
217 200  IF (DABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
218      IF (NST .GE. NSLP+MSBP) IPUP = MITER
219      TN = TN + H
220      I1 = NQNYH + 1
221      DO 215 JB = 1,NQ
222        I1 = I1 - NYH
223CDIR$ IVDEP
224        DO 210 I = I1,NQNYH
225 210      YH1(I) = YH1(I) + YH1(I+NYH)
226 215    CONTINUE
227C-----------------------------------------------------------------------
228C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN.  A CONVERGENCE TEST IS
229C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR
230C WEIGHT VECTOR EWT.  THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE
231C VECTOR ACOR(I).  THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP.
232C-----------------------------------------------------------------------
233 220  M = 0
234      DO 230 I = 1,N
235 230    Y(I) = YH(I,1)
236      IERR = 0
237      CALL F (NEQ, TN, Y, SAVF, IERR)
238      IF (IERR .LT. 0) RETURN
239      NFE = NFE + 1
240      IF (IPUP .LE. 0) GO TO 250
241C-----------------------------------------------------------------------
242C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND
243C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION.  IPUP IS SET
244C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE.
245C-----------------------------------------------------------------------
246      IERR = 0
247      CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC,
248     1   IERR)
249      IF (IERR .LT. 0) RETURN
250      IPUP = 0
251      RC = 1.0D0
252      NSLP = NST
253      CRATE = 0.7D0
254      IF (IERPJ .NE. 0) GO TO 430
255 250  DO 260 I = 1,N
256 260    ACOR(I) = 0.0D0
257 270  IF (MITER .NE. 0) GO TO 350
258C-----------------------------------------------------------------------
259C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM
260C THE RESULT OF THE LAST FUNCTION EVALUATION.
261C-----------------------------------------------------------------------
262      DO 290 I = 1,N
263        SAVF(I) = H*SAVF(I) - YH(I,2)
264 290    Y(I) = SAVF(I) - ACOR(I)
265      DEL = VNORM (N, Y, EWT)
266      DO 300 I = 1,N
267        Y(I) = YH(I,1) + EL(1)*SAVF(I)
268 300    ACOR(I) = SAVF(I)
269      GO TO 400
270C-----------------------------------------------------------------------
271C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
272C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND
273C P AS COEFFICIENT MATRIX.
274C-----------------------------------------------------------------------
275 350  DO 360 I = 1,N
276 360    Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
277      CALL SLVS (WM, IWM, Y, SAVF)
278      IF (IERSL .LT. 0) GO TO 430
279      IF (IERSL .GT. 0) GO TO 410
280      DEL = VNORM (N, Y, EWT)
281      DO 380 I = 1,N
282        ACOR(I) = ACOR(I) + Y(I)
283 380    Y(I) = YH(I,1) + EL(1)*ACOR(I)
284C-----------------------------------------------------------------------
285C TEST FOR CONVERGENCE.  IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
286C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
287C-----------------------------------------------------------------------
288 400  IF (M .NE. 0) CRATE = DMAX1(0.2D0*CRATE,DEL/DELP)
289      DCON = DEL*DMIN1(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
290      IF (DCON .LE. 1.0D0) GO TO 450
291      M = M + 1
292      IF (M .EQ. MAXCOR) GO TO 410
293      IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
294      DELP = DEL
295      IERR = 0
296      CALL F (NEQ, TN, Y, SAVF, IERR)
297      IF (IERR .LT. 0) RETURN
298      NFE = NFE + 1
299      GO TO 270
300C-----------------------------------------------------------------------
301C THE CORRECTOR ITERATION FAILED TO CONVERGE.
302C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR
303C THE NEXT TRY.  OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES
304C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE.  IF H CANNOT BE
305C REDUCED OR MXNCF FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2.
306C-----------------------------------------------------------------------
307 410  IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
308      ICF = 1
309      IPUP = MITER
310      GO TO 220
311 430  ICF = 2
312      NCF = NCF + 1
313      RMAX = 2.0D0
314      TN = TOLD
315      I1 = NQNYH + 1
316      DO 445 JB = 1,NQ
317        I1 = I1 - NYH
318CDIR$ IVDEP
319        DO 440 I = I1,NQNYH
320 440      YH1(I) = YH1(I) - YH1(I+NYH)
321 445    CONTINUE
322      IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
323      IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 670
324      IF (NCF .EQ. MXNCF) GO TO 670
325      RH = 0.25D0
326      IPUP = MITER
327      IREDO = 1
328      GO TO 170
329C-----------------------------------------------------------------------
330C THE CORRECTOR HAS CONVERGED.  JCUR IS SET TO 0
331C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER.
332C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500
333C IF IT FAILS.
334C-----------------------------------------------------------------------
335 450  JCUR = 0
336      IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
337      IF (M .GT. 0) DSM = VNORM (N, ACOR, EWT)/TESCO(2,NQ)
338      IF (DSM .GT. 1.0D0) GO TO 500
339C-----------------------------------------------------------------------
340C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY.
341C CONSIDER CHANGING H IF IALTH = 1.  OTHERWISE DECREASE IALTH BY 1.
342C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR
343C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
344C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
345C BY ONE IS CONSIDERED ALSO.  A CHANGE IN H IS MADE ONLY IF IT IS BY A
346C FACTOR OF AT LEAST 1.1.  IF NOT, IALTH IS SET TO 3 TO PREVENT
347C TESTING FOR THAT MANY STEPS.
348C-----------------------------------------------------------------------
349      KFLAG = 0
350      IREDO = 0
351      NST = NST + 1
352      HU = H
353      NQU = NQ
354      DO 470 J = 1,L
355        DO 470 I = 1,N
356 470      YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
357      IALTH = IALTH - 1
358      IF (IALTH .EQ. 0) GO TO 520
359      IF (IALTH .GT. 1) GO TO 700
360      IF (L .EQ. LMAX) GO TO 700
361      DO 490 I = 1,N
362 490    YH(I,LMAX) = ACOR(I)
363      GO TO 700
364C-----------------------------------------------------------------------
365C THE ERROR TEST FAILED.  KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
366C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
367C TO TRY THE STEP AGAIN.  COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
368C ONE LOWER ORDER.  AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE
369C BY A FACTOR OF 0.2 OR LESS.
370C-----------------------------------------------------------------------
371 500  KFLAG = KFLAG - 1
372      TN = TOLD
373      I1 = NQNYH + 1
374      DO 515 JB = 1,NQ
375        I1 = I1 - NYH
376CDIR$ IVDEP
377        DO 510 I = I1,NQNYH
378 510      YH1(I) = YH1(I) - YH1(I+NYH)
379 515    CONTINUE
380      RMAX = 2.0D0
381      IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 660
382      IF (KFLAG .LE. -3) GO TO 640
383      IREDO = 2
384      RHUP = 0.0D0
385      GO TO 540
386C-----------------------------------------------------------------------
387C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
388C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED
389C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
390C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE.
391C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
392C ACCORDINGLY.  IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
393C ADDITIONAL SCALED DERIVATIVE.
394C-----------------------------------------------------------------------
395 520  RHUP = 0.0D0
396      IF (L .EQ. LMAX) GO TO 540
397      DO 530 I = 1,N
398 530    SAVF(I) = ACOR(I) - YH(I,LMAX)
399      DUP = VNORM (N, SAVF, EWT)/TESCO(3,NQ)
400      EXUP = 1.0D0/DBLE(L+1)
401      RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0)
402 540  EXSM = 1.0D0/DBLE(L)
403      RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
404      RHDN = 0.0D0
405      IF (NQ .EQ. 1) GO TO 560
406      DDN = VNORM (N, YH(1,L), EWT)/TESCO(1,NQ)
407      EXDN = 1.0D0/DBLE(NQ)
408      RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
409 560  IF (RHSM .GE. RHUP) GO TO 570
410      IF (RHUP .GT. RHDN) GO TO 590
411      GO TO 580
412 570  IF (RHSM .LT. RHDN) GO TO 580
413      NEWQ = NQ
414      RH = RHSM
415      GO TO 620
416 580  NEWQ = NQ - 1
417      RH = RHDN
418      IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
419      GO TO 620
420 590  NEWQ = L
421      RH = RHUP
422      IF (RH .LT. 1.1D0) GO TO 610
423      R = EL(L)/DBLE(L)
424      DO 600 I = 1,N
425 600    YH(I,NEWQ+1) = ACOR(I)*R
426      GO TO 630
427 610  IALTH = 3
428      GO TO 700
429 620  IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610
430      IF (KFLAG .LE. -2) RH = DMIN1(RH,0.2D0)
431C-----------------------------------------------------------------------
432C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
433C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED.
434C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
435C-----------------------------------------------------------------------
436      IF (NEWQ .EQ. NQ) GO TO 170
437 630  NQ = NEWQ
438      L = NQ + 1
439      IRET = 2
440      GO TO 150
441C-----------------------------------------------------------------------
442C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURRED.
443C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1.
444C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
445C YH ARRAY HAVE ERRORS OF THE WRONG ORDER.  HENCE THE FIRST
446C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1.  THEN
447C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED,
448C UNTIL IT SUCCEEDS OR H REACHES HMIN.
449C-----------------------------------------------------------------------
450 640  IF (KFLAG .EQ. -10) GO TO 660
451      RH = 0.1D0
452      RH = DMAX1(HMIN/DABS(H),RH)
453      H = H*RH
454      DO 645 I = 1,N
455 645    Y(I) = YH(I,1)
456      IERR = 0
457      CALL F (NEQ, TN, Y, SAVF, IERR)
458      IF (IERR .LT. 0) RETURN
459      NFE = NFE + 1
460      DO 650 I = 1,N
461 650    YH(I,2) = H*SAVF(I)
462      IPUP = MITER
463      IALTH = 5
464      IF (NQ .EQ. 1) GO TO 200
465      NQ = 1
466      L = 2
467      IRET = 3
468      GO TO 150
469C-----------------------------------------------------------------------
470C ALL RETURNS ARE MADE THROUGH THIS SECTION.  H IS SAVED IN HOLD
471C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
472C-----------------------------------------------------------------------
473 660  KFLAG = -1
474      GO TO 720
475 670  KFLAG = -2
476      GO TO 720
477 680  KFLAG = -3
478      GO TO 720
479 690  RMAX = 10.0D0
480 700  R = 1.0D0/TESCO(2,NQU)
481      DO 710 I = 1,N
482 710    ACOR(I) = ACOR(I)*R
483 720  HOLD = H
484      JSTART = 1
485      RETURN
486C----------------------- END OF SUBROUTINE STODE -----------------------
487      END
488