*DECK DRKFS SUBROUTINE DRKFS (DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID, H, + TOLFAC, YP, F1, F2, F3, F4, F5, YS, TOLD, DTSIGN, U26, RER, + INIT, KSTEPS, KOP, IQUIT, STIFF, NONSTF, NTSTEP, NSTIFS, RPAR, + IPAR) C***BEGIN PROLOGUE DRKFS C***SUBSIDIARY C***PURPOSE Subsidiary to DDERKF C***LIBRARY SLATEC C***TYPE DOUBLE PRECISION (DERKFS-S, DRKFS-D) C***AUTHOR Watts, H. A., (SNLA) C***DESCRIPTION C C Fehlberg Fourth-Fifth Order Runge-Kutta Method C ********************************************************************** C C DRKFS integrates a system of first order ordinary differential C equations as described in the comments for DDERKF . C C The arrays YP,F1,F2,F3,F4,F5,and YS (of length at least NEQ) C appear in the call list for variable dimensioning purposes. C C The variables H,TOLFAC,TOLD,DTSIGN,U26,RER,INIT,KSTEPS,KOP,IQUIT, C STIFF,NONSTF,NTSTEP, and NSTIFS are used internally by the code C and appear in the call list to eliminate local retention of C variables between calls. Accordingly, these variables and the C array YP should not be altered. C Items of possible interest are C H - An appropriate step size to be used for the next step C TOLFAC - Factor of change in the tolerances C YP - Derivative of solution vector at T C KSTEPS - Counter on the number of steps attempted C C ********************************************************************** C C***SEE ALSO DDERKF C***ROUTINES CALLED D1MACH, DFEHL, DHSTRT, DHVNRM, XERMSG C***REVISION HISTORY (YYMMDD) C 820301 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 891024 Changed references from DVNORM to DHVNRM. (WRB) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900328 Added TYPE section. (WRB) C 900510 Convert XERRWV calls to XERMSG calls, change GOTOs to C IF-THEN-ELSEs. (RWC) C 910722 Updated AUTHOR section. (ALS) C***END PROLOGUE DRKFS C INTEGER IDID, INFO, INIT, IPAR, IQUIT, K, KOP, KSTEPS, KTOL, 1 MXKOP, MXSTEP, NATOLP, NEQ, NRTOLP, NSTIFS, NTSTEP DOUBLE PRECISION A, ATOL, BIG, D1MACH, 1 DT, DTSIGN, DHVNRM, DY, EE, EEOET, ES, ESTIFF, 2 ESTTOL, ET, F1, F2, F3, F4, F5, H, HMIN, REMIN, RER, RPAR, 3 RTOL, S, T, TOL, TOLD, TOLFAC, TOUT, U, U26, UTE, Y, YAVG, 4 YP, YS LOGICAL HFAILD,OUTPUT,STIFF,NONSTF CHARACTER*8 XERN1 CHARACTER*16 XERN3, XERN4 C DIMENSION Y(*),YP(*),F1(*),F2(*),F3(*),F4(*),F5(*), 1 YS(*),INFO(15),RTOL(*),ATOL(*),RPAR(*),IPAR(*) C EXTERNAL DF C C .................................................................. C C A FIFTH ORDER METHOD WILL GENERALLY NOT BE CAPABLE OF DELIVERING C ACCURACIES NEAR LIMITING PRECISION ON COMPUTERS WITH LONG C WORDLENGTHS. TO PROTECT AGAINST LIMITING PRECISION DIFFICULTIES C ARISING FROM UNREASONABLE ACCURACY REQUESTS, AN APPROPRIATE C TOLERANCE THRESHOLD REMIN IS ASSIGNED FOR THIS METHOD. THIS C VALUE SHOULD NOT BE CHANGED ACROSS DIFFERENT MACHINES. C SAVE REMIN, MXSTEP, MXKOP DATA REMIN /1.0D-12/ C C .................................................................. C C THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE C NUMBER OF STEPS ATTEMPTED. WHEN THIS EXCEEDS MXSTEP, THE C COUNTER IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE C EXCESSIVE WORK. C DATA MXSTEP /500/ C C .................................................................. C C INEFFICIENCY CAUSED BY TOO FREQUENT OUTPUT IS MONITORED BY C COUNTING THE NUMBER OF STEP SIZES WHICH ARE SEVERELY SHORTENED C DUE SOLELY TO THE CHOICE OF OUTPUT POINTS. WHEN THE NUMBER OF C ABUSES EXCEED MXKOP, THE COUNTER IS RESET TO ZERO AND THE USER C IS INFORMED ABOUT POSSIBLE MISUSE OF THE CODE. C DATA MXKOP /100/ C C .................................................................. C C***FIRST EXECUTABLE STATEMENT DRKFS IF (INFO(1) .EQ. 0) THEN C C ON THE FIRST CALL , PERFORM INITIALIZATION -- C DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY U BY CALLING THE C FUNCTION ROUTINE D1MACH. THE USER MUST MAKE SURE THAT THE C VALUES SET IN D1MACH ARE RELEVANT TO THE COMPUTER BEING USED. C U = D1MACH(4) C -- SET ASSOCIATED MACHINE DEPENDENT PARAMETERS U26 = 26.0D0*U RER = 2.0D0*U + REMIN C -- SET TERMINATION FLAG IQUIT = 0 C -- SET INITIALIZATION INDICATOR INIT = 0 C -- SET COUNTER FOR IMPACT OF OUTPUT POINTS KOP = 0 C -- SET COUNTER FOR ATTEMPTED STEPS KSTEPS = 0 C -- SET INDICATORS FOR STIFFNESS DETECTION STIFF = .FALSE. NONSTF = .FALSE. C -- SET STEP COUNTERS FOR STIFFNESS DETECTION NTSTEP = 0 NSTIFS = 0 C -- RESET INFO(1) FOR SUBSEQUENT CALLS INFO(1) = 1 ENDIF C C....................................................................... C C CHECK VALIDITY OF INPUT PARAMETERS ON EACH ENTRY C IF (INFO(1) .NE. 0 .AND. INFO(1) .NE. 1) THEN WRITE (XERN1, '(I8)') INFO(1) CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, INFO(1) MUST BE SET TO 0 ' // * 'FOR THE START OF A NEW PROBLEM, AND MUST BE SET TO 1 ' // * 'FOLLOWING AN INTERRUPTED TASK. YOU ARE ATTEMPTING TO ' // * 'CONTINUE THE INTEGRATION ILLEGALLY BY CALLING THE CODE ' // * 'WITH INFO(1) = ' // XERN1, 3, 1) IDID = -33 ENDIF C IF (INFO(2) .NE. 0 .AND. INFO(2) .NE. 1) THEN WRITE (XERN1, '(I8)') INFO(2) CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, INFO(2) MUST BE 0 OR 1 ' // * 'INDICATING SCALAR AND VECTOR ERROR TOLERANCES, ' // * 'RESPECTIVELY. YOU HAVE CALLED THE CODE WITH INFO(2) = ' // * XERN1, 4, 1) IDID = -33 ENDIF C IF (INFO(3) .NE. 0 .AND. INFO(3) .NE. 1) THEN WRITE (XERN1, '(I8)') INFO(3) CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, INFO(3) MUST BE 0 OR 1 ' // * 'INDICATING THE INTERVAL OR INTERMEDIATE-OUTPUT MODE OF ' // * 'INTEGRATION, RESPECTIVELY. YOU HAVE CALLED THE CODE ' // * 'WITH INFO(3) = ' // XERN1, 5, 1) IDID = -33 ENDIF C IF (NEQ .LT. 1) THEN WRITE (XERN1, '(I8)') NEQ CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, THE NUMBER OF EQUATIONS ' // * 'NEQ MUST BE A POSITIVE INTEGER. YOU HAVE CALLED THE ' // * 'CODE WITH NEQ = ' // XERN1, 6, 1) IDID = -33 ENDIF C NRTOLP = 0 NATOLP = 0 DO 10 K=1,NEQ IF (NRTOLP .EQ. 0 .AND. RTOL(K) .LT. 0.D0) THEN WRITE (XERN1, '(I8)') K WRITE (XERN3, '(1PE15.6)') RTOL(K) CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, THE RELATIVE ERROR ' // * 'TOLERANCES RTOL MUST BE NON-NEGATIVE. YOU HAVE ' // * 'CALLED THE CODE WITH RTOL(' // XERN1 // ') = ' // * XERN3 // '. IN THE CASE OF VECTOR ERROR TOLERANCES, ' // * 'NO FURTHER CHECKING OF RTOL COMPONENTS IS DONE.', 7, 1) IDID = -33 NRTOLP = 1 ENDIF C IF (NATOLP .EQ. 0 .AND. ATOL(K) .LT. 0.D0) THEN WRITE (XERN1, '(I8)') K WRITE (XERN3, '(1PE15.6)') ATOL(K) CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, THE ABSOLUTE ERROR ' // * 'TOLERANCES ATOL MUST BE NON-NEGATIVE. YOU HAVE ' // * 'CALLED THE CODE WITH ATOL(' // XERN1 // ') = ' // * XERN3 // '. IN THE CASE OF VECTOR ERROR TOLERANCES, ' // * 'NO FURTHER CHECKING OF ATOL COMPONENTS IS DONE.', 8, 1) IDID = -33 NATOLP = 1 ENDIF C IF (INFO(2) .EQ. 0) GO TO 20 IF (NATOLP.GT.0 .AND. NRTOLP.GT.0) GO TO 20 10 CONTINUE C C C CHECK SOME CONTINUATION POSSIBILITIES C 20 IF (INIT .NE. 0) THEN IF (T .EQ. TOUT) THEN WRITE (XERN3, '(1PE15.6)') T CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, YOU HAVE CALLED THE ' // * 'CODE WITH T = TOUT = ' // XERN3 // '$$THIS IS NOT ' // * 'ALLOWED ON CONTINUATION CALLS.', 9, 1) IDID=-33 ENDIF C IF (T .NE. TOLD) THEN WRITE (XERN3, '(1PE15.6)') TOLD WRITE (XERN4, '(1PE15.6)') T CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, YOU HAVE CHANGED THE ' // * 'VALUE OF T FROM ' // XERN3 // ' TO ' // XERN4 // * '$$THIS IS NOT ALLOWED ON CONTINUATION CALLS.', 10, 1) IDID=-33 ENDIF C IF (INIT .NE. 1) THEN IF (DTSIGN*(TOUT-T) .LT. 0.D0) THEN WRITE (XERN3, '(1PE15.6)') TOUT CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, BY CALLING THE CODE WITH TOUT = ' // * XERN3 // ' YOU ARE ATTEMPTING TO CHANGE THE ' // * 'DIRECTION OF INTEGRATION.$$THIS IS NOT ALLOWED ' // * 'WITHOUT RESTARTING.', 11, 1) IDID=-33 ENDIF ENDIF ENDIF C C INVALID INPUT DETECTED C IF (IDID .EQ. (-33)) THEN IF (IQUIT .NE. (-33)) THEN IQUIT = -33 GOTO 540 ELSE CALL XERMSG ('SLATEC', 'DRKFS', * 'IN DDERKF, INVALID INPUT WAS ' // * 'DETECTED ON SUCCESSIVE ENTRIES. IT IS IMPOSSIBLE ' // * 'TO PROCEED BECAUSE YOU HAVE NOT CORRECTED THE ' // * 'PROBLEM, SO EXECUTION IS BEING TERMINATED.', 12, 2) RETURN ENDIF ENDIF C C ............................................................ C C RTOL = ATOL = 0. IS ALLOWED AS VALID INPUT AND C INTERPRETED AS ASKING FOR THE MOST ACCURATE SOLUTION C POSSIBLE. IN THIS CASE, THE RELATIVE ERROR TOLERANCE C RTOL IS RESET TO THE SMALLEST VALUE RER WHICH IS LIKELY C TO BE REASONABLE FOR THIS METHOD AND MACHINE. C DO 190 K = 1, NEQ IF (RTOL(K) + ATOL(K) .GT. 0.0D0) GO TO 180 RTOL(K) = RER IDID = -2 180 CONTINUE C ...EXIT IF (INFO(2) .EQ. 0) GO TO 200 190 CONTINUE 200 CONTINUE C IF (IDID .NE. (-2)) GO TO 210 C C RTOL=ATOL=0 ON INPUT, SO RTOL WAS CHANGED TO A C SMALL POSITIVE VALUE TOLFAC = 1.0D0 GO TO 530 210 CONTINUE C C BRANCH ON STATUS OF INITIALIZATION INDICATOR C INIT=0 MEANS INITIAL DERIVATIVES AND C STARTING STEP SIZE C NOT YET COMPUTED C INIT=1 MEANS STARTING STEP SIZE NOT YET C COMPUTED INIT=2 MEANS NO FURTHER C INITIALIZATION REQUIRED C IF (INIT .EQ. 0) GO TO 220 C ......EXIT IF (INIT .EQ. 1) GO TO 240 C .........EXIT GO TO 260 220 CONTINUE C C ................................................ C C MORE INITIALIZATION -- C -- EVALUATE INITIAL C DERIVATIVES C INIT = 1 A = T CALL DF(A,Y,YP,RPAR,IPAR) IF (T .NE. TOUT) GO TO 230 C C INTERVAL MODE IDID = 2 T = TOUT TOLD = T C .....................EXIT GO TO 560 230 CONTINUE 240 CONTINUE C C -- SET SIGN OF INTEGRATION DIRECTION AND C -- ESTIMATE STARTING STEP SIZE C INIT = 2 DTSIGN = SIGN(1.0D0,TOUT-T) U = D1MACH(4) BIG = SQRT(D1MACH(2)) UTE = U**0.375D0 DY = UTE*DHVNRM(Y,NEQ) IF (DY .EQ. 0.0D0) DY = UTE KTOL = 1 DO 250 K = 1, NEQ IF (INFO(2) .EQ. 1) KTOL = K TOL = RTOL(KTOL)*ABS(Y(K)) + ATOL(KTOL) IF (TOL .EQ. 0.0D0) TOL = DY*RTOL(KTOL) F1(K) = TOL 250 CONTINUE C CALL DHSTRT(DF,NEQ,T,TOUT,Y,YP,F1,4,U,BIG,F2,F3,F4, 1 F5,RPAR,IPAR,H) 260 CONTINUE C C ...................................................... C C SET STEP SIZE FOR INTEGRATION IN THE DIRECTION C FROM T TO TOUT AND SET OUTPUT POINT INDICATOR C DT = TOUT - T H = SIGN(H,DT) OUTPUT = .FALSE. C C TEST TO SEE IF DDERKF IS BEING SEVERELY IMPACTED BY C TOO MANY OUTPUT POINTS C IF (ABS(H) .GE. 2.0D0*ABS(DT)) KOP = KOP + 1 IF (KOP .LE. MXKOP) GO TO 270 C C UNNECESSARY FREQUENCY OF OUTPUT IS RESTRICTING C THE STEP SIZE CHOICE IDID = -5 KOP = 0 GO TO 510 270 CONTINUE C IF (ABS(DT) .GT. U26*ABS(T)) GO TO 290 C C IF TOO CLOSE TO OUTPUT POINT,EXTRAPOLATE AND C RETURN C DO 280 K = 1, NEQ Y(K) = Y(K) + DT*YP(K) 280 CONTINUE A = TOUT CALL DF(A,Y,YP,RPAR,IPAR) KSTEPS = KSTEPS + 1 GO TO 500 290 CONTINUE C BEGIN BLOCK PERMITTING ...EXITS TO 490 C C ********************************************* C ********************************************* C STEP BY STEP INTEGRATION C 300 CONTINUE C BEGIN BLOCK PERMITTING ...EXITS TO 480 HFAILD = .FALSE. C C TO PROTECT AGAINST IMPOSSIBLE ACCURACY C REQUESTS, COMPUTE A TOLERANCE FACTOR C BASED ON THE REQUESTED ERROR TOLERANCE C AND A LEVEL OF ACCURACY ACHIEVABLE AT C LIMITING PRECISION C TOLFAC = 0.0D0 KTOL = 1 DO 330 K = 1, NEQ IF (INFO(2) .EQ. 1) KTOL = K ET = RTOL(KTOL)*ABS(Y(K)) 1 + ATOL(KTOL) IF (ET .GT. 0.0D0) GO TO 310 TOLFAC = MAX(TOLFAC, 1 RER/RTOL(KTOL)) GO TO 320 310 CONTINUE TOLFAC = MAX(TOLFAC, 1 ABS(Y(K)) 2 *(RER/ET)) 320 CONTINUE 330 CONTINUE IF (TOLFAC .LE. 1.0D0) GO TO 340 C C REQUESTED ERROR UNATTAINABLE DUE TO LIMITED C PRECISION AVAILABLE TOLFAC = 2.0D0*TOLFAC IDID = -2 C .....................EXIT GO TO 520 340 CONTINUE C C SET SMALLEST ALLOWABLE STEP SIZE C HMIN = U26*ABS(T) C C ADJUST STEP SIZE IF NECESSARY TO HIT C THE OUTPUT POINT -- LOOK AHEAD TWO C STEPS TO AVOID DRASTIC CHANGES IN THE C STEP SIZE AND THUS LESSEN THE IMPACT OF C OUTPUT POINTS ON THE CODE. STRETCH THE C STEP SIZE BY, AT MOST, AN AMOUNT EQUAL C TO THE SAFETY FACTOR OF 9/10. C DT = TOUT - T IF (ABS(DT) .GE. 2.0D0*ABS(H)) 1 GO TO 370 IF (ABS(DT) .GT. ABS(H)/0.9D0) 1 GO TO 350 C C THE NEXT STEP, IF SUCCESSFUL, C WILL COMPLETE THE INTEGRATION TO C THE OUTPUT POINT C OUTPUT = .TRUE. H = DT GO TO 360 350 CONTINUE C H = 0.5D0*DT 360 CONTINUE 370 CONTINUE C C C *************************************** C CORE INTEGRATOR FOR TAKING A C SINGLE STEP C *************************************** C TO AVOID PROBLEMS WITH ZERO C CROSSINGS, RELATIVE ERROR IS C MEASURED USING THE AVERAGE OF THE C MAGNITUDES OF THE SOLUTION AT THE C BEGINNING AND END OF A STEP. C THE ERROR ESTIMATE FORMULA HAS C BEEN GROUPED TO CONTROL LOSS OF C SIGNIFICANCE. C LOCAL ERROR ESTIMATES FOR A FIRST C ORDER METHOD USING THE SAME C STEP SIZE AS THE FEHLBERG METHOD C ARE CALCULATED AS PART OF THE C TEST FOR STIFFNESS. C TO DISTINGUISH THE VARIOUS C ARGUMENTS, H IS NOT PERMITTED C TO BECOME SMALLER THAN 26 UNITS OF C ROUNDOFF IN T. PRACTICAL LIMITS C ON THE CHANGE IN THE STEP SIZE ARE C ENFORCED TO SMOOTH THE STEP SIZE C SELECTION PROCESS AND TO AVOID C EXCESSIVE CHATTERING ON PROBLEMS C HAVING DISCONTINUITIES. TO C PREVENT UNNECESSARY FAILURES, THE C CODE USES 9/10 THE STEP SIZE C IT ESTIMATES WILL SUCCEED. C AFTER A STEP FAILURE, THE STEP C SIZE IS NOT ALLOWED TO INCREASE C FOR THE NEXT ATTEMPTED STEP. THIS C MAKES THE CODE MORE EFFICIENT ON C PROBLEMS HAVING DISCONTINUITIES C AND MORE EFFECTIVE IN GENERAL C SINCE LOCAL EXTRAPOLATION IS BEING C USED AND EXTRA CAUTION SEEMS C WARRANTED. C ....................................... C C MONITOR NUMBER OF STEPS ATTEMPTED C 380 CONTINUE IF (KSTEPS .LE. MXSTEP) GO TO 390 C C A SIGNIFICANT AMOUNT OF WORK HAS C BEEN EXPENDED IDID = -1 KSTEPS = 0 C ........................EXIT IF (.NOT.STIFF) GO TO 520 C C PROBLEM APPEARS TO BE STIFF IDID = -4 STIFF = .FALSE. NONSTF = .FALSE. NTSTEP = 0 NSTIFS = 0 C ........................EXIT GO TO 520 390 CONTINUE C C ADVANCE AN APPROXIMATE SOLUTION OVER C ONE STEP OF LENGTH H C CALL DFEHL(DF,NEQ,T,Y,H,YP,F1,F2,F3, 1 F4,F5,YS,RPAR,IPAR) KSTEPS = KSTEPS + 1 C C .................................... C C COMPUTE AND TEST ALLOWABLE C TOLERANCES VERSUS LOCAL ERROR C ESTIMATES. NOTE THAT RELATIVE C ERROR IS MEASURED WITH RESPECT C TO THE AVERAGE OF THE C MAGNITUDES OF THE SOLUTION AT C THE BEGINNING AND END OF THE C STEP. LOCAL ERROR ESTIMATES C FOR A SPECIAL FIRST ORDER C METHOD ARE CALCULATED ONLY WHEN C THE STIFFNESS DETECTION IS C TURNED ON. C EEOET = 0.0D0 ESTIFF = 0.0D0 KTOL = 1 DO 420 K = 1, NEQ YAVG = 0.5D0 1 *(ABS(Y(K)) 2 + ABS(YS(K))) IF (INFO(2) .EQ. 1) KTOL = K ET = RTOL(KTOL)*YAVG + ATOL(KTOL) IF (ET .GT. 0.0D0) GO TO 400 C C PURE RELATIVE ERROR INAPPROPRIATE WHEN SOLUTION C VANISHES IDID = -3 C ...........................EXIT GO TO 520 400 CONTINUE C EE = ABS((-2090.0D0*YP(K) 1 +(21970.0D0*F3(K) 2 -15048.0D0*F4(K))) 3 +(22528.0D0*F2(K) 4 -27360.0D0*F5(K))) IF (STIFF .OR. NONSTF) GO TO 410 ES = ABS(H 1 *(0.055455D0*YP(K) 2 -0.035493D0*F1(K) 3 -0.036571D0*F2(K) 4 +0.023107D0*F3(K) 5 -0.009515D0*F4(K) 6 +0.003017D0*F5(K)) 7 ) ESTIFF = MAX(ESTIFF,ES/ET) 410 CONTINUE EEOET = MAX(EEOET,EE/ET) 420 CONTINUE C ESTTOL = ABS(H)*EEOET/752400.0D0 C C ...EXIT IF (ESTTOL .LE. 1.0D0) GO TO 440 C C .................................... C C UNSUCCESSFUL STEP C IF (ABS(H) .GT. HMIN) GO TO 430 C C REQUESTED ERROR UNATTAINABLE AT SMALLEST C ALLOWABLE STEP SIZE TOLFAC = 1.69D0*ESTTOL IDID = -2 C ........................EXIT GO TO 520 430 CONTINUE C C REDUCE THE STEP SIZE , TRY AGAIN C THE DECREASE IS LIMITED TO A FACTOR C OF 1/10 C HFAILD = .TRUE. OUTPUT = .FALSE. S = 0.1D0 IF (ESTTOL .LT. 59049.0D0) 1 S = 0.9D0/ESTTOL**0.2D0 H = SIGN(MAX(S*ABS(H),HMIN),H) GO TO 380 440 CONTINUE C C ....................................... C C SUCCESSFUL STEP C STORE SOLUTION AT T+H C AND EVALUATE C DERIVATIVES THERE C T = T + H DO 450 K = 1, NEQ Y(K) = YS(K) 450 CONTINUE A = T CALL DF(A,Y,YP,RPAR,IPAR) C C CHOOSE NEXT STEP SIZE C THE INCREASE IS LIMITED TO A FACTOR OF C 5 IF STEP FAILURE HAS JUST OCCURRED, C NEXT C STEP SIZE IS NOT ALLOWED TO INCREASE C S = 5.0D0 IF (ESTTOL .GT. 1.889568D-4) 1 S = 0.9D0/ESTTOL**0.2D0 IF (HFAILD) S = MIN(S,1.0D0) H = SIGN(MAX(S*ABS(H),HMIN),H) C C ....................................... C C CHECK FOR STIFFNESS (IF NOT C ALREADY DETECTED) C C IN A SEQUENCE OF 50 SUCCESSFUL C STEPS BY THE FEHLBERG METHOD, 25 C SUCCESSFUL STEPS BY THE FIRST C ORDER METHOD INDICATES STIFFNESS C AND TURNS THE TEST OFF. IF 26 C FAILURES BY THE FIRST ORDER METHOD C OCCUR, THE TEST IS TURNED OFF C UNTIL THIS SEQUENCE OF 50 STEPS BY C THE FEHLBERG METHOD IS COMPLETED. C C ...EXIT IF (STIFF) GO TO 480 NTSTEP = MOD(NTSTEP+1,50) IF (NTSTEP .EQ. 1) NONSTF = .FALSE. C ...EXIT IF (NONSTF) GO TO 480 IF (ESTIFF .GT. 1.0D0) GO TO 460 C C SUCCESSFUL STEP WITH FIRST ORDER C METHOD NSTIFS = NSTIFS + 1 C TURN TEST OFF AFTER 25 INDICATIONS C OF STIFFNESS IF (NSTIFS .EQ. 25) STIFF = .TRUE. GO TO 470 460 CONTINUE C C UNSUCCESSFUL STEP WITH FIRST ORDER C METHOD IF (NTSTEP - NSTIFS .LE. 25) GO TO 470 C TURN STIFFNESS DETECTION OFF FOR THIS BLOCK OF C FIFTY STEPS NONSTF = .TRUE. C RESET STIFF STEP COUNTER NSTIFS = 0 470 CONTINUE 480 CONTINUE C C ****************************************** C END OF CORE INTEGRATOR C ****************************************** C C C SHOULD WE TAKE ANOTHER STEP C C ......EXIT IF (OUTPUT) GO TO 490 IF (INFO(3) .EQ. 0) GO TO 300 C C ********************************************* C ********************************************* C C INTEGRATION SUCCESSFULLY COMPLETED C C ONE-STEP MODE IDID = 1 TOLD = T C .....................EXIT GO TO 560 490 CONTINUE 500 CONTINUE C C INTERVAL MODE IDID = 2 T = TOUT TOLD = T C ...............EXIT GO TO 560 510 CONTINUE 520 CONTINUE 530 CONTINUE 540 CONTINUE C C INTEGRATION TASK INTERRUPTED C INFO(1) = -1 TOLD = T C ...EXIT IF (IDID .NE. (-2)) GO TO 560 C C THE ERROR TOLERANCES ARE INCREASED TO VALUES C WHICH ARE APPROPRIATE FOR CONTINUING RTOL(1) = TOLFAC*RTOL(1) ATOL(1) = TOLFAC*ATOL(1) C ...EXIT IF (INFO(2) .EQ. 0) GO TO 560 DO 550 K = 2, NEQ RTOL(K) = TOLFAC*RTOL(K) ATOL(K) = TOLFAC*ATOL(K) 550 CONTINUE 560 CONTINUE RETURN END